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