001:       SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
002:      $                   LDU, C, LDC, RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               D( * ), E( * ), RWORK( * )
015:       COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CBDSQR computes the singular values and, optionally, the right and/or
022: *  left singular vectors from the singular value decomposition (SVD) of
023: *  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
024: *  zero-shift QR algorithm.  The SVD of B has the form
025: *  
026: *     B = Q * S * P**H
027: *  
028: *  where S is the diagonal matrix of singular values, Q is an orthogonal
029: *  matrix of left singular vectors, and P is an orthogonal matrix of
030: *  right singular vectors.  If left singular vectors are requested, this
031: *  subroutine actually returns U*Q instead of Q, and, if right singular
032: *  vectors are requested, this subroutine returns P**H*VT instead of
033: *  P**H, for given complex input matrices U and VT.  When U and VT are
034: *  the unitary matrices that reduce a general matrix A to bidiagonal
035: *  form: A = U*B*VT, as computed by CGEBRD, then
036: *  
037: *     A = (U*Q) * S * (P**H*VT)
038: *  
039: *  is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
040: *  for a given complex input matrix C.
041: *
042: *  See "Computing  Small Singular Values of Bidiagonal Matrices With
043: *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
044: *  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
045: *  no. 5, pp. 873-912, Sept 1990) and
046: *  "Accurate singular values and differential qd algorithms," by
047: *  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
048: *  Department, University of California at Berkeley, July 1992
049: *  for a detailed description of the algorithm.
050: *
051: *  Arguments
052: *  =========
053: *
054: *  UPLO    (input) CHARACTER*1
055: *          = 'U':  B is upper bidiagonal;
056: *          = 'L':  B is lower bidiagonal.
057: *
058: *  N       (input) INTEGER
059: *          The order of the matrix B.  N >= 0.
060: *
061: *  NCVT    (input) INTEGER
062: *          The number of columns of the matrix VT. NCVT >= 0.
063: *
064: *  NRU     (input) INTEGER
065: *          The number of rows of the matrix U. NRU >= 0.
066: *
067: *  NCC     (input) INTEGER
068: *          The number of columns of the matrix C. NCC >= 0.
069: *
070: *  D       (input/output) REAL array, dimension (N)
071: *          On entry, the n diagonal elements of the bidiagonal matrix B.
072: *          On exit, if INFO=0, the singular values of B in decreasing
073: *          order.
074: *
075: *  E       (input/output) REAL array, dimension (N-1)
076: *          On entry, the N-1 offdiagonal elements of the bidiagonal
077: *          matrix B.
078: *          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
079: *          will contain the diagonal and superdiagonal elements of a
080: *          bidiagonal matrix orthogonally equivalent to the one given
081: *          as input.
082: *
083: *  VT      (input/output) COMPLEX array, dimension (LDVT, NCVT)
084: *          On entry, an N-by-NCVT matrix VT.
085: *          On exit, VT is overwritten by P**H * VT.
086: *          Not referenced if NCVT = 0.
087: *
088: *  LDVT    (input) INTEGER
089: *          The leading dimension of the array VT.
090: *          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
091: *
092: *  U       (input/output) COMPLEX array, dimension (LDU, N)
093: *          On entry, an NRU-by-N matrix U.
094: *          On exit, U is overwritten by U * Q.
095: *          Not referenced if NRU = 0.
096: *
097: *  LDU     (input) INTEGER
098: *          The leading dimension of the array U.  LDU >= max(1,NRU).
099: *
100: *  C       (input/output) COMPLEX 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) REAL 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  REAL, 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:       REAL               ZERO
148:       PARAMETER          ( ZERO = 0.0E0 )
149:       REAL               ONE
150:       PARAMETER          ( ONE = 1.0E0 )
151:       REAL               NEGONE
152:       PARAMETER          ( NEGONE = -1.0E0 )
153:       REAL               HNDRTH
154:       PARAMETER          ( HNDRTH = 0.01E0 )
155:       REAL               TEN
156:       PARAMETER          ( TEN = 10.0E0 )
157:       REAL               HNDRD
158:       PARAMETER          ( HNDRD = 100.0E0 )
159:       REAL               MEIGTH
160:       PARAMETER          ( MEIGTH = -0.125E0 )
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:       REAL               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:       REAL               SLAMCH
176:       EXTERNAL           LSAME, SLAMCH
177: *     ..
178: *     .. External Subroutines ..
179:       EXTERNAL           CLASR, CSROT, CSSCAL, CSWAP, SLARTG, SLAS2,
180:      $                   SLASQ1, SLASV2, XERBLA
181: *     ..
182: *     .. Intrinsic Functions ..
183:       INTRINSIC          ABS, MAX, MIN, REAL, 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( 'CBDSQR', -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 SLASQ1( 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 = SLAMCH( 'Epsilon' )
238:       UNFL = SLAMCH( '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 SLARTG( 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 CLASR( 'R', 'V', 'F', NRU, N, RWORK( 1 ), RWORK( N ),
257:      $                  U, LDU )
258:          IF( NCC.GT.0 )
259:      $      CALL CLASR( '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( REAL( 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 SLASV2( 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 CSROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
378:      $                  COSR, SINR )
379:          IF( NRU.GT.0 )
380:      $      CALL CSROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
381:          IF( NCC.GT.0 )
382:      $      CALL CSROT( 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 SLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
482:          ELSE
483:             SLL = ABS( D( M ) )
484:             CALL SLAS2( 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 SLARTG( D( I )*CS, E( I ), CS, SN, R )
511:                IF( I.GT.LL )
512:      $            E( I-1 ) = OLDSN*R
513:                CALL SLARTG( 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 CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
527:      $                     RWORK( N ), VT( LL, 1 ), LDVT )
528:             IF( NRU.GT.0 )
529:      $         CALL CLASR( '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 CLASR( '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 SLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
549:                IF( I.LT.M )
550:      $            E( I ) = OLDSN*R
551:                CALL SLARTG( 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 CLASR( '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 CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
568:      $                     RWORK( N ), U( 1, LL ), LDU )
569:             IF( NCC.GT.0 )
570:      $         CALL CLASR( '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 SLARTG( 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 SLARTG( 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 CLASR( 'L', 'V', 'F', M-LL+1, NCVT, RWORK( 1 ),
617:      $                     RWORK( N ), VT( LL, 1 ), LDVT )
618:             IF( NRU.GT.0 )
619:      $         CALL CLASR( '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 CLASR( '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 SLARTG( 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 SLARTG( 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 CLASR( '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 CLASR( 'R', 'V', 'B', NRU, M-LL+1, RWORK( 1 ),
673:      $                     RWORK( N ), U( 1, LL ), LDU )
674:             IF( NCC.GT.0 )
675:      $         CALL CLASR( '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 CSSCAL( 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 CSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
721:      $                     LDVT )
722:             IF( NRU.GT.0 )
723:      $         CALL CSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
724:             IF( NCC.GT.0 )
725:      $         CALL CSWAP( 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 CBDSQR
742: *
743:       END
744: