Mercurial > octave-nkf

annotate libcruft/lapack/zlals0.f @ 7072:b48d486f641d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-26 15:52:57 by jwe]
author jwe
date Fri, 26 Oct 2007 15:52:58 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7072
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
1 SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
2 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
3 $ POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
4 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
5 * -- LAPACK routine (version 3.1) --
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
7 * November 2006
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
8 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
9 * .. Scalar Arguments ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
10 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
11 $ LDGNUM, NL, NR, NRHS, SQRE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
12 DOUBLE PRECISION C, S
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
13 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
14 * .. Array Arguments ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
15 INTEGER GIVCOL( LDGCOL, * ), PERM( * )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
16 DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
17 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
18 $ RWORK( * ), Z( * )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
19 COMPLEX*16 B( LDB, * ), BX( LDBX, * )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
20 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
21 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
22 * Purpose
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
23 * =======
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
24 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
25 * ZLALS0 applies back the multiplying factors of either the left or the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
26 * right singular vector matrix of a diagonal matrix appended by a row
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
27 * to the right hand side matrix B in solving the least squares problem
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
28 * using the divide-and-conquer SVD approach.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
29 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
30 * For the left singular vector matrix, three types of orthogonal
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
31 * matrices are involved:
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
32 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
33 * (1L) Givens rotations: the number of such rotations is GIVPTR; the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
34 * pairs of columns/rows they were applied to are stored in GIVCOL;
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
35 * and the C- and S-values of these rotations are stored in GIVNUM.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
36 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
37 * (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
38 * row, and for J=2:N, PERM(J)-th row of B is to be moved to the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
39 * J-th row.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
40 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
41 * (3L) The left singular vector matrix of the remaining matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
42 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
43 * For the right singular vector matrix, four types of orthogonal
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
44 * matrices are involved:
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
45 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
46 * (1R) The right singular vector matrix of the remaining matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
47 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
48 * (2R) If SQRE = 1, one extra Givens rotation to generate the right
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
49 * null space.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
50 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
51 * (3R) The inverse transformation of (2L).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
52 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
53 * (4R) The inverse transformation of (1L).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
54 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
55 * Arguments
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
56 * =========
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
57 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
58 * ICOMPQ (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
59 * Specifies whether singular vectors are to be computed in
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
60 * factored form:
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
61 * = 0: Left singular vector matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
62 * = 1: Right singular vector matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
63 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
64 * NL (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
65 * The row dimension of the upper block. NL >= 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
66 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
67 * NR (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
68 * The row dimension of the lower block. NR >= 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
69 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
70 * SQRE (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
71 * = 0: the lower block is an NR-by-NR square matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
72 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
73 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
74 * The bidiagonal matrix has row dimension N = NL + NR + 1,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
75 * and column dimension M = N + SQRE.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
76 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
77 * NRHS (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
78 * The number of columns of B and BX. NRHS must be at least 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
79 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
80 * B (input/output) COMPLEX*16 array, dimension ( LDB, NRHS )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
81 * On input, B contains the right hand sides of the least
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
82 * squares problem in rows 1 through M. On output, B contains
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
83 * the solution X in rows 1 through N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
84 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
85 * LDB (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
86 * The leading dimension of B. LDB must be at least
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
87 * max(1,MAX( M, N ) ).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
88 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
89 * BX (workspace) COMPLEX*16 array, dimension ( LDBX, NRHS )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
90 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
91 * LDBX (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
92 * The leading dimension of BX.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
93 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
94 * PERM (input) INTEGER array, dimension ( N )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
95 * The permutations (from deflation and sorting) applied
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
96 * to the two blocks.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
97 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
98 * GIVPTR (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
99 * The number of Givens rotations which took place in this
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
100 * subproblem.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
101 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
102 * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
103 * Each pair of numbers indicates a pair of rows/columns
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
104 * involved in a Givens rotation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
105 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
106 * LDGCOL (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
107 * The leading dimension of GIVCOL, must be at least N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
108 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
109 * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
110 * Each number indicates the C or S value used in the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
111 * corresponding Givens rotation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
112 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
113 * LDGNUM (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
114 * The leading dimension of arrays DIFR, POLES and
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
115 * GIVNUM, must be at least K.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
116 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
117 * POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
118 * On entry, POLES(1:K, 1) contains the new singular
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
119 * values obtained from solving the secular equation, and
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
120 * POLES(1:K, 2) is an array containing the poles in the secular
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
121 * equation.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
122 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
123 * DIFL (input) DOUBLE PRECISION array, dimension ( K ).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
124 * On entry, DIFL(I) is the distance between I-th updated
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
125 * (undeflated) singular value and the I-th (undeflated) old
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
126 * singular value.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
127 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
128 * DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
129 * On entry, DIFR(I, 1) contains the distances between I-th
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
130 * updated (undeflated) singular value and the I+1-th
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
131 * (undeflated) old singular value. And DIFR(I, 2) is the
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
132 * normalizing factor for the I-th right singular vector.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
133 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
134 * Z (input) DOUBLE PRECISION array, dimension ( K )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
135 * Contain the components of the deflation-adjusted updating row
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
136 * vector.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
137 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
138 * K (input) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
139 * Contains the dimension of the non-deflated matrix,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
140 * This is the order of the related secular equation. 1 <= K <=N.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
141 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
142 * C (input) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
143 * C contains garbage if SQRE =0 and the C-value of a Givens
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
144 * rotation related to the right null space if SQRE = 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
145 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
146 * S (input) DOUBLE PRECISION
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
147 * S contains garbage if SQRE =0 and the S-value of a Givens
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
148 * rotation related to the right null space if SQRE = 1.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
149 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
150 * RWORK (workspace) DOUBLE PRECISION array, dimension
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
151 * ( K*(1+NRHS) + 2*NRHS )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
152 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
153 * INFO (output) INTEGER
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
154 * = 0: successful exit.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
155 * < 0: if INFO = -i, the i-th argument had an illegal value.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
156 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
157 * Further Details
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
158 * ===============
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
159 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
160 * Based on contributions by
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
161 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
162 * California at Berkeley, USA
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
163 * Osni Marques, LBNL/NERSC, USA
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
164 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
165 * =====================================================================
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
166 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
167 * .. Parameters ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
168 DOUBLE PRECISION ONE, ZERO, NEGONE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
169 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
170 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
171 * .. Local Scalars ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
172 INTEGER I, J, JCOL, JROW, M, N, NLP1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
173 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
174 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
175 * .. External Subroutines ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
176 EXTERNAL DGEMV, XERBLA, ZCOPY, ZDROT, ZDSCAL, ZLACPY,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
177 $ ZLASCL
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
178 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
179 * .. External Functions ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
180 DOUBLE PRECISION DLAMC3, DNRM2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
181 EXTERNAL DLAMC3, DNRM2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
182 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
183 * .. Intrinsic Functions ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
184 INTRINSIC DBLE, DCMPLX, DIMAG, MAX
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
185 * ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
186 * .. Executable Statements ..
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
187 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
188 * Test the input parameters.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
189 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
190 INFO = 0
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
191 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
192 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
193 INFO = -1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
194 ELSE IF( NL.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
195 INFO = -2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
196 ELSE IF( NR.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
197 INFO = -3
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
198 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
199 INFO = -4
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
200 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
201 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
202 N = NL + NR + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
203 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
204 IF( NRHS.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
205 INFO = -5
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
206 ELSE IF( LDB.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
207 INFO = -7
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
208 ELSE IF( LDBX.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
209 INFO = -9
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
210 ELSE IF( GIVPTR.LT.0 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
211 INFO = -11
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
212 ELSE IF( LDGCOL.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
213 INFO = -13
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
214 ELSE IF( LDGNUM.LT.N ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
215 INFO = -15
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
216 ELSE IF( K.LT.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
217 INFO = -20
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
218 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
219 IF( INFO.NE.0 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
220 CALL XERBLA( 'ZLALS0', -INFO )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
221 RETURN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
222 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
223 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
224 M = N + SQRE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
225 NLP1 = NL + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
226 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
227 IF( ICOMPQ.EQ.0 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
228 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
229 * Apply back orthogonal transformations from the left.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
230 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
231 * Step (1L): apply back the Givens rotations performed.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
232 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
233 DO 10 I = 1, GIVPTR
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
234 CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
235 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
236 $ GIVNUM( I, 1 ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
237 10 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
238 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
239 * Step (2L): permute rows of B.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
240 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
241 CALL ZCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
242 DO 20 I = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
243 CALL ZCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
244 20 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
245 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
246 * Step (3L): apply the inverse of the left singular vector
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
247 * matrix to BX.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
248 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
249 IF( K.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
250 CALL ZCOPY( NRHS, BX, LDBX, B, LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
251 IF( Z( 1 ).LT.ZERO ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
252 CALL ZDSCAL( NRHS, NEGONE, B, LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
253 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
254 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
255 DO 100 J = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
256 DIFLJ = DIFL( J )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
257 DJ = POLES( J, 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
258 DSIGJ = -POLES( J, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
259 IF( J.LT.K ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
260 DIFRJ = -DIFR( J, 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
261 DSIGJP = -POLES( J+1, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
262 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
263 IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
264 $ THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
265 RWORK( J ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
266 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
267 RWORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
268 $ ( POLES( J, 2 )+DJ )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
269 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
270 DO 30 I = 1, J - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
271 IF( ( Z( I ).EQ.ZERO ) .OR.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
272 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
273 RWORK( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
274 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
275 RWORK( I ) = POLES( I, 2 )*Z( I ) /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
276 $ ( DLAMC3( POLES( I, 2 ), DSIGJ )-
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
277 $ DIFLJ ) / ( POLES( I, 2 )+DJ )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
278 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
279 30 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
280 DO 40 I = J + 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
281 IF( ( Z( I ).EQ.ZERO ) .OR.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
282 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
283 RWORK( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
284 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
285 RWORK( I ) = POLES( I, 2 )*Z( I ) /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
286 $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
287 $ DIFRJ ) / ( POLES( I, 2 )+DJ )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
288 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
289 40 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
290 RWORK( 1 ) = NEGONE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
291 TEMP = DNRM2( K, RWORK, 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
292 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
293 * Since B and BX are complex, the following call to DGEMV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
294 * is performed in two steps (real and imaginary parts).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
295 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
296 * CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
297 * $ B( J, 1 ), LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
298 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
299 I = K + NRHS*2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
300 DO 60 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
301 DO 50 JROW = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
302 I = I + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
303 RWORK( I ) = DBLE( BX( JROW, JCOL ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
304 50 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
305 60 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
306 CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
307 $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
308 I = K + NRHS*2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
309 DO 80 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
310 DO 70 JROW = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
311 I = I + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
312 RWORK( I ) = DIMAG( BX( JROW, JCOL ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
313 70 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
314 80 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
315 CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
316 $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
317 DO 90 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
318 B( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
319 $ RWORK( JCOL+K+NRHS ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
320 90 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
321 CALL ZLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
322 $ LDB, INFO )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
323 100 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
324 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
325 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
326 * Move the deflated rows of BX to B also.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
327 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
328 IF( K.LT.MAX( M, N ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
329 $ CALL ZLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
330 $ B( K+1, 1 ), LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
331 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
332 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
333 * Apply back the right orthogonal transformations.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
334 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
335 * Step (1R): apply back the new right singular vector matrix
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
336 * to B.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
337 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
338 IF( K.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
339 CALL ZCOPY( NRHS, B, LDB, BX, LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
340 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
341 DO 180 J = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
342 DSIGJ = POLES( J, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
343 IF( Z( J ).EQ.ZERO ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
344 RWORK( J ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
345 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
346 RWORK( J ) = -Z( J ) / DIFL( J ) /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
347 $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
348 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
349 DO 110 I = 1, J - 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
350 IF( Z( J ).EQ.ZERO ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
351 RWORK( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
352 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
353 RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
354 $ 2 ) )-DIFR( I, 1 ) ) /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
355 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
356 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
357 110 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
358 DO 120 I = J + 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
359 IF( Z( J ).EQ.ZERO ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
360 RWORK( I ) = ZERO
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
361 ELSE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
362 RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
363 $ 2 ) )-DIFL( I ) ) /
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
364 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
365 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
366 120 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
367 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
368 * Since B and BX are complex, the following call to DGEMV
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
369 * is performed in two steps (real and imaginary parts).
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
370 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
371 * CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
372 * $ BX( J, 1 ), LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
373 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
374 I = K + NRHS*2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
375 DO 140 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
376 DO 130 JROW = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
377 I = I + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
378 RWORK( I ) = DBLE( B( JROW, JCOL ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
379 130 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
380 140 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
381 CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
382 $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
383 I = K + NRHS*2
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
384 DO 160 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
385 DO 150 JROW = 1, K
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
386 I = I + 1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
387 RWORK( I ) = DIMAG( B( JROW, JCOL ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
388 150 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
389 160 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
390 CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
391 $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
392 DO 170 JCOL = 1, NRHS
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
393 BX( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
394 $ RWORK( JCOL+K+NRHS ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
395 170 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
396 180 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
397 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
398 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
399 * Step (2R): if SQRE = 1, apply back the rotation that is
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
400 * related to the right null space of the subproblem.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
401 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
402 IF( SQRE.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
403 CALL ZCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
404 CALL ZDROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
405 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
406 IF( K.LT.MAX( M, N ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
407 $ CALL ZLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
408 $ LDBX )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
409 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
410 * Step (3R): permute rows of B.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
411 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
412 CALL ZCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
413 IF( SQRE.EQ.1 ) THEN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
414 CALL ZCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
415 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
416 DO 190 I = 2, N
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
417 CALL ZCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
418 190 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
419 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
420 * Step (4R): apply back the Givens rotations performed.
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
421 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
422 DO 200 I = GIVPTR, 1, -1
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
423 CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
424 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
425 $ -GIVNUM( I, 1 ) )
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
426 200 CONTINUE
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
427 END IF
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
428 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
429 RETURN
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
430 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
431 * End of ZLALS0
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
432 *
b48d486f641d [project @ 2007-10-26 15:52:57 by jwe]
jwe
parents:
diff changeset
433 END