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