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