|
2329
|
1 SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, |
|
|
2 $ LDU, C, LDC, RWORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK routine (version 3.0) -- |
|
2329
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3596
|
7 * October 31, 1999 |
|
2329
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER UPLO |
|
|
11 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION D( * ), E( * ), RWORK( * ) |
|
|
15 COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * ZBDSQR computes the singular value decomposition (SVD) of a real |
|
|
22 * N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P' |
|
|
23 * denotes the transpose of P), where S is a diagonal matrix with |
|
|
24 * non-negative diagonal elements (the singular values of B), and Q |
|
|
25 * and P are orthogonal matrices. |
|
|
26 * |
|
|
27 * The routine computes S, and optionally computes U * Q, P' * VT, |
|
|
28 * or Q' * C, for given complex input matrices U, VT, and C. |
|
|
29 * |
|
|
30 * See "Computing Small Singular Values of Bidiagonal Matrices With |
|
|
31 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, |
|
|
32 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, |
|
|
33 * no. 5, pp. 873-912, Sept 1990) and |
|
|
34 * "Accurate singular values and differential qd algorithms," by |
|
|
35 * B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics |
|
|
36 * Department, University of California at Berkeley, July 1992 |
|
|
37 * for a detailed description of the algorithm. |
|
|
38 * |
|
|
39 * Arguments |
|
|
40 * ========= |
|
|
41 * |
|
|
42 * UPLO (input) CHARACTER*1 |
|
|
43 * = 'U': B is upper bidiagonal; |
|
|
44 * = 'L': B is lower bidiagonal. |
|
|
45 * |
|
|
46 * N (input) INTEGER |
|
|
47 * The order of the matrix B. N >= 0. |
|
|
48 * |
|
|
49 * NCVT (input) INTEGER |
|
|
50 * The number of columns of the matrix VT. NCVT >= 0. |
|
|
51 * |
|
|
52 * NRU (input) INTEGER |
|
|
53 * The number of rows of the matrix U. NRU >= 0. |
|
|
54 * |
|
|
55 * NCC (input) INTEGER |
|
|
56 * The number of columns of the matrix C. NCC >= 0. |
|
|
57 * |
|
|
58 * D (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
59 * On entry, the n diagonal elements of the bidiagonal matrix B. |
|
|
60 * On exit, if INFO=0, the singular values of B in decreasing |
|
|
61 * order. |
|
|
62 * |
|
|
63 * E (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
64 * On entry, the elements of E contain the |
|
|
65 * offdiagonal elements of of the bidiagonal matrix whose SVD |
|
|
66 * is desired. On normal exit (INFO = 0), E is destroyed. |
|
|
67 * If the algorithm does not converge (INFO > 0), D and E |
|
|
68 * will contain the diagonal and superdiagonal elements of a |
|
|
69 * bidiagonal matrix orthogonally equivalent to the one given |
|
|
70 * as input. E(N) is used for workspace. |
|
|
71 * |
|
|
72 * VT (input/output) COMPLEX*16 array, dimension (LDVT, NCVT) |
|
|
73 * On entry, an N-by-NCVT matrix VT. |
|
|
74 * On exit, VT is overwritten by P' * VT. |
|
|
75 * VT is not referenced if NCVT = 0. |
|
|
76 * |
|
|
77 * LDVT (input) INTEGER |
|
|
78 * The leading dimension of the array VT. |
|
|
79 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. |
|
|
80 * |
|
|
81 * U (input/output) COMPLEX*16 array, dimension (LDU, N) |
|
|
82 * On entry, an NRU-by-N matrix U. |
|
|
83 * On exit, U is overwritten by U * Q. |
|
|
84 * U is not referenced if NRU = 0. |
|
|
85 * |
|
|
86 * LDU (input) INTEGER |
|
|
87 * The leading dimension of the array U. LDU >= max(1,NRU). |
|
|
88 * |
|
|
89 * C (input/output) COMPLEX*16 array, dimension (LDC, NCC) |
|
|
90 * On entry, an N-by-NCC matrix C. |
|
|
91 * On exit, C is overwritten by Q' * C. |
|
|
92 * C is not referenced if NCC = 0. |
|
|
93 * |
|
|
94 * LDC (input) INTEGER |
|
|
95 * The leading dimension of the array C. |
|
|
96 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. |
|
|
97 * |
|
3596
|
98 * RWORK (workspace) DOUBLE PRECISION array, dimension (4*N) |
|
2329
|
99 * |
|
|
100 * INFO (output) INTEGER |
|
|
101 * = 0: successful exit |
|
|
102 * < 0: If INFO = -i, the i-th argument had an illegal value |
|
|
103 * > 0: the algorithm did not converge; D and E contain the |
|
|
104 * elements of a bidiagonal matrix which is orthogonally |
|
|
105 * similar to the input matrix B; if INFO = i, i |
|
|
106 * elements of E have not converged to zero. |
|
|
107 * |
|
|
108 * Internal Parameters |
|
|
109 * =================== |
|
|
110 * |
|
|
111 * TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) |
|
|
112 * TOLMUL controls the convergence criterion of the QR loop. |
|
|
113 * If it is positive, TOLMUL*EPS is the desired relative |
|
|
114 * precision in the computed singular values. |
|
|
115 * If it is negative, abs(TOLMUL*EPS*sigma_max) is the |
|
|
116 * desired absolute accuracy in the computed singular |
|
|
117 * values (corresponds to relative accuracy |
|
|
118 * abs(TOLMUL*EPS) in the largest singular value. |
|
|
119 * abs(TOLMUL) should be between 1 and 1/EPS, and preferably |
|
|
120 * between 10 (for fast convergence) and .1/EPS |
|
|
121 * (for there to be some accuracy in the results). |
|
|
122 * Default is to lose at either one eighth or 2 of the |
|
|
123 * available decimal digits in each computed singular value |
|
|
124 * (whichever is smaller). |
|
|
125 * |
|
|
126 * MAXITR INTEGER, default = 6 |
|
|
127 * MAXITR controls the maximum number of passes of the |
|
|
128 * algorithm through its inner loop. The algorithms stops |
|
|
129 * (and so fails to converge) if the number of passes |
|
|
130 * through the inner loop exceeds MAXITR*N**2. |
|
|
131 * |
|
|
132 * ===================================================================== |
|
|
133 * |
|
|
134 * .. Parameters .. |
|
|
135 DOUBLE PRECISION ZERO |
|
|
136 PARAMETER ( ZERO = 0.0D0 ) |
|
|
137 DOUBLE PRECISION ONE |
|
|
138 PARAMETER ( ONE = 1.0D0 ) |
|
|
139 DOUBLE PRECISION NEGONE |
|
|
140 PARAMETER ( NEGONE = -1.0D0 ) |
|
|
141 DOUBLE PRECISION HNDRTH |
|
|
142 PARAMETER ( HNDRTH = 0.01D0 ) |
|
|
143 DOUBLE PRECISION TEN |
|
|
144 PARAMETER ( TEN = 10.0D0 ) |
|
|
145 DOUBLE PRECISION HNDRD |
|
|
146 PARAMETER ( HNDRD = 100.0D0 ) |
|
|
147 DOUBLE PRECISION MEIGTH |
|
|
148 PARAMETER ( MEIGTH = -0.125D0 ) |
|
|
149 INTEGER MAXITR |
|
|
150 PARAMETER ( MAXITR = 6 ) |
|
|
151 * .. |
|
|
152 * .. Local Scalars .. |
|
3333
|
153 LOGICAL LOWER, ROTATE |
|
|
154 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1, |
|
|
155 $ NM12, NM13, OLDLL, OLDM |
|
2329
|
156 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU, |
|
|
157 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL, |
|
|
158 $ SINR, SLL, SMAX, SMIN, SMINL, SMINLO, SMINOA, |
|
|
159 $ SN, THRESH, TOL, TOLMUL, UNFL |
|
|
160 * .. |
|
|
161 * .. External Functions .. |
|
|
162 LOGICAL LSAME |
|
|
163 DOUBLE PRECISION DLAMCH |
|
|
164 EXTERNAL LSAME, DLAMCH |
|
|
165 * .. |
|
|
166 * .. External Subroutines .. |
|
|
167 EXTERNAL DLARTG, DLAS2, DLASQ1, DLASV2, XERBLA, ZDROT, |
|
|
168 $ ZDSCAL, ZLASR, ZSWAP |
|
|
169 * .. |
|
|
170 * .. Intrinsic Functions .. |
|
|
171 INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT |
|
|
172 * .. |
|
|
173 * .. Executable Statements .. |
|
|
174 * |
|
|
175 * Test the input parameters. |
|
|
176 * |
|
|
177 INFO = 0 |
|
3333
|
178 LOWER = LSAME( UPLO, 'L' ) |
|
|
179 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN |
|
2329
|
180 INFO = -1 |
|
|
181 ELSE IF( N.LT.0 ) THEN |
|
|
182 INFO = -2 |
|
|
183 ELSE IF( NCVT.LT.0 ) THEN |
|
|
184 INFO = -3 |
|
|
185 ELSE IF( NRU.LT.0 ) THEN |
|
|
186 INFO = -4 |
|
|
187 ELSE IF( NCC.LT.0 ) THEN |
|
|
188 INFO = -5 |
|
|
189 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. |
|
|
190 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN |
|
|
191 INFO = -9 |
|
|
192 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN |
|
|
193 INFO = -11 |
|
|
194 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. |
|
|
195 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN |
|
|
196 INFO = -13 |
|
|
197 END IF |
|
|
198 IF( INFO.NE.0 ) THEN |
|
|
199 CALL XERBLA( 'ZBDSQR', -INFO ) |
|
|
200 RETURN |
|
|
201 END IF |
|
|
202 IF( N.EQ.0 ) |
|
|
203 $ RETURN |
|
|
204 IF( N.EQ.1 ) |
|
3333
|
205 $ GO TO 160 |
|
2329
|
206 * |
|
|
207 * ROTATE is true if any singular vectors desired, false otherwise |
|
|
208 * |
|
|
209 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) |
|
|
210 * |
|
|
211 * If no singular vectors desired, use qd algorithm |
|
|
212 * |
|
|
213 IF( .NOT.ROTATE ) THEN |
|
|
214 CALL DLASQ1( N, D, E, RWORK, INFO ) |
|
|
215 RETURN |
|
|
216 END IF |
|
|
217 * |
|
|
218 NM1 = N - 1 |
|
|
219 NM12 = NM1 + NM1 |
|
|
220 NM13 = NM12 + NM1 |
|
3333
|
221 IDIR = 0 |
|
2329
|
222 * |
|
|
223 * Get machine constants |
|
|
224 * |
|
|
225 EPS = DLAMCH( 'Epsilon' ) |
|
|
226 UNFL = DLAMCH( 'Safe minimum' ) |
|
|
227 * |
|
|
228 * If matrix lower bidiagonal, rotate to be upper bidiagonal |
|
|
229 * by applying Givens rotations on the left |
|
|
230 * |
|
3333
|
231 IF( LOWER ) THEN |
|
2329
|
232 DO 10 I = 1, N - 1 |
|
|
233 CALL DLARTG( D( I ), E( I ), CS, SN, R ) |
|
|
234 D( I ) = R |
|
|
235 E( I ) = SN*D( I+1 ) |
|
|
236 D( I+1 ) = CS*D( I+1 ) |
|
|
237 RWORK( I ) = CS |
|
|
238 RWORK( NM1+I ) = SN |
|
|
239 10 CONTINUE |
|
|
240 * |
|
|
241 * Update singular vectors if desired |
|
|
242 * |
|
|
243 IF( NRU.GT.0 ) |
|
|
244 $ CALL ZLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ), |
|
|
245 $ U, LDU ) |
|
|
246 IF( NCC.GT.0 ) |
|
|
247 $ CALL ZLASR( 'L', 'V', 'F', N, NCC, RWORK( 1 ), RWORK( N ), |
|
|
248 $ C, LDC ) |
|
|
249 END IF |
|
|
250 * |
|
|
251 * Compute singular values to relative accuracy TOL |
|
|
252 * (By setting TOL to be negative, algorithm will compute |
|
|
253 * singular values to absolute accuracy ABS(TOL)*norm(input matrix)) |
|
|
254 * |
|
|
255 TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) ) |
|
|
256 TOL = TOLMUL*EPS |
|
|
257 * |
|
|
258 * Compute approximate maximum, minimum singular values |
|
|
259 * |
|
3333
|
260 SMAX = ZERO |
|
|
261 DO 20 I = 1, N |
|
|
262 SMAX = MAX( SMAX, ABS( D( I ) ) ) |
|
2329
|
263 20 CONTINUE |
|
3333
|
264 DO 30 I = 1, N - 1 |
|
|
265 SMAX = MAX( SMAX, ABS( E( I ) ) ) |
|
|
266 30 CONTINUE |
|
2329
|
267 SMINL = ZERO |
|
|
268 IF( TOL.GE.ZERO ) THEN |
|
|
269 * |
|
|
270 * Relative accuracy desired |
|
|
271 * |
|
|
272 SMINOA = ABS( D( 1 ) ) |
|
|
273 IF( SMINOA.EQ.ZERO ) |
|
3333
|
274 $ GO TO 50 |
|
2329
|
275 MU = SMINOA |
|
3333
|
276 DO 40 I = 2, N |
|
2329
|
277 MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) ) |
|
|
278 SMINOA = MIN( SMINOA, MU ) |
|
|
279 IF( SMINOA.EQ.ZERO ) |
|
3333
|
280 $ GO TO 50 |
|
2329
|
281 40 CONTINUE |
|
3333
|
282 50 CONTINUE |
|
2329
|
283 SMINOA = SMINOA / SQRT( DBLE( N ) ) |
|
|
284 THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL ) |
|
|
285 ELSE |
|
|
286 * |
|
|
287 * Absolute accuracy desired |
|
|
288 * |
|
|
289 THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL ) |
|
|
290 END IF |
|
|
291 * |
|
|
292 * Prepare for main iteration loop for the singular values |
|
|
293 * (MAXIT is the maximum number of passes through the inner |
|
|
294 * loop permitted before nonconvergence signalled.) |
|
|
295 * |
|
|
296 MAXIT = MAXITR*N*N |
|
|
297 ITER = 0 |
|
|
298 OLDLL = -1 |
|
|
299 OLDM = -1 |
|
|
300 * |
|
|
301 * M points to last element of unconverged part of matrix |
|
|
302 * |
|
|
303 M = N |
|
|
304 * |
|
|
305 * Begin main iteration loop |
|
|
306 * |
|
3333
|
307 60 CONTINUE |
|
2329
|
308 * |
|
|
309 * Check for convergence or exceeding iteration count |
|
|
310 * |
|
|
311 IF( M.LE.1 ) |
|
3333
|
312 $ GO TO 160 |
|
2329
|
313 IF( ITER.GT.MAXIT ) |
|
3333
|
314 $ GO TO 200 |
|
2329
|
315 * |
|
|
316 * Find diagonal block of matrix to work on |
|
|
317 * |
|
|
318 IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH ) |
|
|
319 $ D( M ) = ZERO |
|
|
320 SMAX = ABS( D( M ) ) |
|
|
321 SMIN = SMAX |
|
3333
|
322 DO 70 LLL = 1, M - 1 |
|
2329
|
323 LL = M - LLL |
|
|
324 ABSS = ABS( D( LL ) ) |
|
|
325 ABSE = ABS( E( LL ) ) |
|
|
326 IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) |
|
|
327 $ D( LL ) = ZERO |
|
|
328 IF( ABSE.LE.THRESH ) |
|
3333
|
329 $ GO TO 80 |
|
2329
|
330 SMIN = MIN( SMIN, ABSS ) |
|
|
331 SMAX = MAX( SMAX, ABSS, ABSE ) |
|
|
332 70 CONTINUE |
|
3333
|
333 LL = 0 |
|
|
334 GO TO 90 |
|
|
335 80 CONTINUE |
|
2329
|
336 E( LL ) = ZERO |
|
|
337 * |
|
|
338 * Matrix splits since E(LL) = 0 |
|
|
339 * |
|
|
340 IF( LL.EQ.M-1 ) THEN |
|
|
341 * |
|
|
342 * Convergence of bottom singular value, return to top of loop |
|
|
343 * |
|
|
344 M = M - 1 |
|
3333
|
345 GO TO 60 |
|
2329
|
346 END IF |
|
3333
|
347 90 CONTINUE |
|
2329
|
348 LL = LL + 1 |
|
|
349 * |
|
|
350 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero |
|
|
351 * |
|
|
352 IF( LL.EQ.M-1 ) THEN |
|
|
353 * |
|
|
354 * 2 by 2 block, handle separately |
|
|
355 * |
|
|
356 CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR, |
|
|
357 $ COSR, SINL, COSL ) |
|
|
358 D( M-1 ) = SIGMX |
|
|
359 E( M-1 ) = ZERO |
|
|
360 D( M ) = SIGMN |
|
|
361 * |
|
|
362 * Compute singular vectors, if desired |
|
|
363 * |
|
|
364 IF( NCVT.GT.0 ) |
|
|
365 $ CALL ZDROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, |
|
|
366 $ COSR, SINR ) |
|
|
367 IF( NRU.GT.0 ) |
|
|
368 $ CALL ZDROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL ) |
|
|
369 IF( NCC.GT.0 ) |
|
|
370 $ CALL ZDROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL, |
|
|
371 $ SINL ) |
|
|
372 M = M - 2 |
|
3333
|
373 GO TO 60 |
|
2329
|
374 END IF |
|
|
375 * |
|
|
376 * If working on new submatrix, choose shift direction |
|
|
377 * (from larger end diagonal element towards smaller) |
|
|
378 * |
|
|
379 IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN |
|
|
380 IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN |
|
|
381 * |
|
|
382 * Chase bulge from top (big end) to bottom (small end) |
|
|
383 * |
|
|
384 IDIR = 1 |
|
|
385 ELSE |
|
|
386 * |
|
|
387 * Chase bulge from bottom (big end) to top (small end) |
|
|
388 * |
|
|
389 IDIR = 2 |
|
|
390 END IF |
|
|
391 END IF |
|
|
392 * |
|
|
393 * Apply convergence tests |
|
|
394 * |
|
|
395 IF( IDIR.EQ.1 ) THEN |
|
|
396 * |
|
|
397 * Run convergence test in forward direction |
|
|
398 * First apply standard test to bottom of matrix |
|
|
399 * |
|
|
400 IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR. |
|
|
401 $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN |
|
|
402 E( M-1 ) = ZERO |
|
3333
|
403 GO TO 60 |
|
2329
|
404 END IF |
|
|
405 * |
|
|
406 IF( TOL.GE.ZERO ) THEN |
|
|
407 * |
|
|
408 * If relative accuracy desired, |
|
|
409 * apply convergence criterion forward |
|
|
410 * |
|
|
411 MU = ABS( D( LL ) ) |
|
|
412 SMINL = MU |
|
3333
|
413 DO 100 LLL = LL, M - 1 |
|
2329
|
414 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN |
|
|
415 E( LLL ) = ZERO |
|
3333
|
416 GO TO 60 |
|
2329
|
417 END IF |
|
|
418 SMINLO = SMINL |
|
|
419 MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) |
|
|
420 SMINL = MIN( SMINL, MU ) |
|
3333
|
421 100 CONTINUE |
|
2329
|
422 END IF |
|
|
423 * |
|
|
424 ELSE |
|
|
425 * |
|
|
426 * Run convergence test in backward direction |
|
|
427 * First apply standard test to top of matrix |
|
|
428 * |
|
|
429 IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR. |
|
|
430 $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN |
|
|
431 E( LL ) = ZERO |
|
3333
|
432 GO TO 60 |
|
2329
|
433 END IF |
|
|
434 * |
|
|
435 IF( TOL.GE.ZERO ) THEN |
|
|
436 * |
|
|
437 * If relative accuracy desired, |
|
|
438 * apply convergence criterion backward |
|
|
439 * |
|
|
440 MU = ABS( D( M ) ) |
|
|
441 SMINL = MU |
|
3333
|
442 DO 110 LLL = M - 1, LL, -1 |
|
2329
|
443 IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN |
|
|
444 E( LLL ) = ZERO |
|
3333
|
445 GO TO 60 |
|
2329
|
446 END IF |
|
|
447 SMINLO = SMINL |
|
|
448 MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) ) |
|
|
449 SMINL = MIN( SMINL, MU ) |
|
3333
|
450 110 CONTINUE |
|
2329
|
451 END IF |
|
|
452 END IF |
|
|
453 OLDLL = LL |
|
|
454 OLDM = M |
|
|
455 * |
|
|
456 * Compute shift. First, test if shifting would ruin relative |
|
|
457 * accuracy, and if so set the shift to zero. |
|
|
458 * |
|
|
459 IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE. |
|
|
460 $ MAX( EPS, HNDRTH*TOL ) ) THEN |
|
|
461 * |
|
|
462 * Use a zero shift to avoid loss of relative accuracy |
|
|
463 * |
|
|
464 SHIFT = ZERO |
|
|
465 ELSE |
|
|
466 * |
|
|
467 * Compute the shift from 2-by-2 block at end of matrix |
|
|
468 * |
|
|
469 IF( IDIR.EQ.1 ) THEN |
|
|
470 SLL = ABS( D( LL ) ) |
|
|
471 CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R ) |
|
|
472 ELSE |
|
|
473 SLL = ABS( D( M ) ) |
|
|
474 CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R ) |
|
|
475 END IF |
|
|
476 * |
|
|
477 * Test if shift negligible, and if so set to zero |
|
|
478 * |
|
|
479 IF( SLL.GT.ZERO ) THEN |
|
|
480 IF( ( SHIFT / SLL )**2.LT.EPS ) |
|
|
481 $ SHIFT = ZERO |
|
|
482 END IF |
|
|
483 END IF |
|
|
484 * |
|
|
485 * Increment iteration count |
|
|
486 * |
|
|
487 ITER = ITER + M - LL |
|
|
488 * |
|
|
489 * If SHIFT = 0, do simplified QR iteration |
|
|
490 * |
|
|
491 IF( SHIFT.EQ.ZERO ) THEN |
|
|
492 IF( IDIR.EQ.1 ) THEN |
|
|
493 * |
|
|
494 * Chase bulge from top to bottom |
|
|
495 * Save cosines and sines for later singular vector updates |
|
|
496 * |
|
|
497 CS = ONE |
|
|
498 OLDCS = ONE |
|
3333
|
499 DO 120 I = LL, M - 1 |
|
2329
|
500 CALL DLARTG( D( I )*CS, E( I ), CS, SN, R ) |
|
3333
|
501 IF( I.GT.LL ) |
|
|
502 $ E( I-1 ) = OLDSN*R |
|
2329
|
503 CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) ) |
|
3333
|
504 RWORK( I-LL+1 ) = CS |
|
|
505 RWORK( I-LL+1+NM1 ) = SN |
|
|
506 RWORK( I-LL+1+NM12 ) = OLDCS |
|
|
507 RWORK( I-LL+1+NM13 ) = OLDSN |
|
|
508 120 CONTINUE |
|
2329
|
509 H = D( M )*CS |
|
|
510 D( M ) = H*OLDCS |
|
|
511 E( M-1 ) = H*OLDSN |
|
|
512 * |
|
|
513 * Update singular vectors |
|
|
514 * |
|
|
515 IF( NCVT.GT.0 ) |
|
|
516 $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), |
|
|
517 $ RWORK( N ), VT( LL, 1 ), LDVT ) |
|
|
518 IF( NRU.GT.0 ) |
|
|
519 $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), |
|
|
520 $ RWORK( NM13+1 ), U( 1, LL ), LDU ) |
|
|
521 IF( NCC.GT.0 ) |
|
|
522 $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), |
|
|
523 $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) |
|
|
524 * |
|
|
525 * Test convergence |
|
|
526 * |
|
|
527 IF( ABS( E( M-1 ) ).LE.THRESH ) |
|
|
528 $ E( M-1 ) = ZERO |
|
|
529 * |
|
|
530 ELSE |
|
|
531 * |
|
|
532 * Chase bulge from bottom to top |
|
|
533 * Save cosines and sines for later singular vector updates |
|
|
534 * |
|
|
535 CS = ONE |
|
|
536 OLDCS = ONE |
|
3333
|
537 DO 130 I = M, LL + 1, -1 |
|
2329
|
538 CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R ) |
|
3333
|
539 IF( I.LT.M ) |
|
|
540 $ E( I ) = OLDSN*R |
|
2329
|
541 CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) ) |
|
3333
|
542 RWORK( I-LL ) = CS |
|
|
543 RWORK( I-LL+NM1 ) = -SN |
|
|
544 RWORK( I-LL+NM12 ) = OLDCS |
|
|
545 RWORK( I-LL+NM13 ) = -OLDSN |
|
|
546 130 CONTINUE |
|
2329
|
547 H = D( LL )*CS |
|
|
548 D( LL ) = H*OLDCS |
|
|
549 E( LL ) = H*OLDSN |
|
|
550 * |
|
|
551 * Update singular vectors |
|
|
552 * |
|
|
553 IF( NCVT.GT.0 ) |
|
|
554 $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), |
|
|
555 $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) |
|
|
556 IF( NRU.GT.0 ) |
|
|
557 $ CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), |
|
|
558 $ RWORK( N ), U( 1, LL ), LDU ) |
|
|
559 IF( NCC.GT.0 ) |
|
|
560 $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ), |
|
|
561 $ RWORK( N ), C( LL, 1 ), LDC ) |
|
|
562 * |
|
|
563 * Test convergence |
|
|
564 * |
|
|
565 IF( ABS( E( LL ) ).LE.THRESH ) |
|
|
566 $ E( LL ) = ZERO |
|
|
567 END IF |
|
|
568 ELSE |
|
|
569 * |
|
|
570 * Use nonzero shift |
|
|
571 * |
|
|
572 IF( IDIR.EQ.1 ) THEN |
|
|
573 * |
|
|
574 * Chase bulge from top to bottom |
|
|
575 * Save cosines and sines for later singular vector updates |
|
|
576 * |
|
|
577 F = ( ABS( D( LL ) )-SHIFT )* |
|
|
578 $ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) ) |
|
|
579 G = E( LL ) |
|
3333
|
580 DO 140 I = LL, M - 1 |
|
2329
|
581 CALL DLARTG( F, G, COSR, SINR, R ) |
|
3333
|
582 IF( I.GT.LL ) |
|
|
583 $ E( I-1 ) = R |
|
2329
|
584 F = COSR*D( I ) + SINR*E( I ) |
|
|
585 E( I ) = COSR*E( I ) - SINR*D( I ) |
|
|
586 G = SINR*D( I+1 ) |
|
|
587 D( I+1 ) = COSR*D( I+1 ) |
|
|
588 CALL DLARTG( F, G, COSL, SINL, R ) |
|
|
589 D( I ) = R |
|
|
590 F = COSL*E( I ) + SINL*D( I+1 ) |
|
|
591 D( I+1 ) = COSL*D( I+1 ) - SINL*E( I ) |
|
3333
|
592 IF( I.LT.M-1 ) THEN |
|
|
593 G = SINL*E( I+1 ) |
|
|
594 E( I+1 ) = COSL*E( I+1 ) |
|
|
595 END IF |
|
|
596 RWORK( I-LL+1 ) = COSR |
|
|
597 RWORK( I-LL+1+NM1 ) = SINR |
|
|
598 RWORK( I-LL+1+NM12 ) = COSL |
|
|
599 RWORK( I-LL+1+NM13 ) = SINL |
|
|
600 140 CONTINUE |
|
2329
|
601 E( M-1 ) = F |
|
|
602 * |
|
|
603 * Update singular vectors |
|
|
604 * |
|
|
605 IF( NCVT.GT.0 ) |
|
|
606 $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ), |
|
|
607 $ RWORK( N ), VT( LL, 1 ), LDVT ) |
|
|
608 IF( NRU.GT.0 ) |
|
|
609 $ CALL ZLASR( 'R', 'V', 'F', NRU, M-LL+1, RWORK( NM12+1 ), |
|
|
610 $ RWORK( NM13+1 ), U( 1, LL ), LDU ) |
|
|
611 IF( NCC.GT.0 ) |
|
|
612 $ CALL ZLASR( 'L', 'V', 'F', M-LL+1, NCC, RWORK( NM12+1 ), |
|
|
613 $ RWORK( NM13+1 ), C( LL, 1 ), LDC ) |
|
|
614 * |
|
|
615 * Test convergence |
|
|
616 * |
|
|
617 IF( ABS( E( M-1 ) ).LE.THRESH ) |
|
|
618 $ E( M-1 ) = ZERO |
|
|
619 * |
|
|
620 ELSE |
|
|
621 * |
|
|
622 * Chase bulge from bottom to top |
|
|
623 * Save cosines and sines for later singular vector updates |
|
|
624 * |
|
|
625 F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT / |
|
|
626 $ D( M ) ) |
|
|
627 G = E( M-1 ) |
|
3333
|
628 DO 150 I = M, LL + 1, -1 |
|
2329
|
629 CALL DLARTG( F, G, COSR, SINR, R ) |
|
3333
|
630 IF( I.LT.M ) |
|
|
631 $ E( I ) = R |
|
2329
|
632 F = COSR*D( I ) + SINR*E( I-1 ) |
|
|
633 E( I-1 ) = COSR*E( I-1 ) - SINR*D( I ) |
|
|
634 G = SINR*D( I-1 ) |
|
|
635 D( I-1 ) = COSR*D( I-1 ) |
|
|
636 CALL DLARTG( F, G, COSL, SINL, R ) |
|
|
637 D( I ) = R |
|
|
638 F = COSL*E( I-1 ) + SINL*D( I-1 ) |
|
|
639 D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 ) |
|
3333
|
640 IF( I.GT.LL+1 ) THEN |
|
|
641 G = SINL*E( I-2 ) |
|
|
642 E( I-2 ) = COSL*E( I-2 ) |
|
|
643 END IF |
|
|
644 RWORK( I-LL ) = COSR |
|
|
645 RWORK( I-LL+NM1 ) = -SINR |
|
|
646 RWORK( I-LL+NM12 ) = COSL |
|
|
647 RWORK( I-LL+NM13 ) = -SINL |
|
|
648 150 CONTINUE |
|
2329
|
649 E( LL ) = F |
|
|
650 * |
|
|
651 * Test convergence |
|
|
652 * |
|
|
653 IF( ABS( E( LL ) ).LE.THRESH ) |
|
|
654 $ E( LL ) = ZERO |
|
|
655 * |
|
|
656 * Update singular vectors if desired |
|
|
657 * |
|
|
658 IF( NCVT.GT.0 ) |
|
|
659 $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCVT, RWORK( NM12+1 ), |
|
|
660 $ RWORK( NM13+1 ), VT( LL, 1 ), LDVT ) |
|
|
661 IF( NRU.GT.0 ) |
|
|
662 $ CALL ZLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ), |
|
|
663 $ RWORK( N ), U( 1, LL ), LDU ) |
|
|
664 IF( NCC.GT.0 ) |
|
|
665 $ CALL ZLASR( 'L', 'V', 'B', M-LL+1, NCC, RWORK( 1 ), |
|
|
666 $ RWORK( N ), C( LL, 1 ), LDC ) |
|
|
667 END IF |
|
|
668 END IF |
|
|
669 * |
|
|
670 * QR iteration finished, go back and check convergence |
|
|
671 * |
|
3333
|
672 GO TO 60 |
|
2329
|
673 * |
|
|
674 * All singular values converged, so make them positive |
|
|
675 * |
|
3333
|
676 160 CONTINUE |
|
|
677 DO 170 I = 1, N |
|
2329
|
678 IF( D( I ).LT.ZERO ) THEN |
|
|
679 D( I ) = -D( I ) |
|
|
680 * |
|
|
681 * Change sign of singular vectors, if desired |
|
|
682 * |
|
|
683 IF( NCVT.GT.0 ) |
|
|
684 $ CALL ZDSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT ) |
|
|
685 END IF |
|
3333
|
686 170 CONTINUE |
|
2329
|
687 * |
|
|
688 * Sort the singular values into decreasing order (insertion sort on |
|
|
689 * singular values, but only one transposition per singular vector) |
|
|
690 * |
|
3333
|
691 DO 190 I = 1, N - 1 |
|
2329
|
692 * |
|
|
693 * Scan for smallest D(I) |
|
|
694 * |
|
|
695 ISUB = 1 |
|
|
696 SMIN = D( 1 ) |
|
3333
|
697 DO 180 J = 2, N + 1 - I |
|
2329
|
698 IF( D( J ).LE.SMIN ) THEN |
|
|
699 ISUB = J |
|
|
700 SMIN = D( J ) |
|
|
701 END IF |
|
3333
|
702 180 CONTINUE |
|
2329
|
703 IF( ISUB.NE.N+1-I ) THEN |
|
|
704 * |
|
|
705 * Swap singular values and vectors |
|
|
706 * |
|
|
707 D( ISUB ) = D( N+1-I ) |
|
|
708 D( N+1-I ) = SMIN |
|
|
709 IF( NCVT.GT.0 ) |
|
|
710 $ CALL ZSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ), |
|
|
711 $ LDVT ) |
|
|
712 IF( NRU.GT.0 ) |
|
|
713 $ CALL ZSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 ) |
|
|
714 IF( NCC.GT.0 ) |
|
|
715 $ CALL ZSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC ) |
|
|
716 END IF |
|
3333
|
717 190 CONTINUE |
|
|
718 GO TO 220 |
|
2329
|
719 * |
|
|
720 * Maximum number of iterations exceeded, failure to converge |
|
|
721 * |
|
3333
|
722 200 CONTINUE |
|
2329
|
723 INFO = 0 |
|
3333
|
724 DO 210 I = 1, N - 1 |
|
2329
|
725 IF( E( I ).NE.ZERO ) |
|
|
726 $ INFO = INFO + 1 |
|
|
727 210 CONTINUE |
|
3333
|
728 220 CONTINUE |
|
2329
|
729 RETURN |
|
|
730 * |
|
|
731 * End of ZBDSQR |
|
|
732 * |
|
|
733 END |
