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