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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     January 2007
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
015:      $                   VT( LDVT, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SBDSQR 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**T
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**T*VT instead of
033: *  P**T, for given real input matrices U and VT.  When U and VT are the
034: *  orthogonal matrices that reduce a general matrix A to bidiagonal
035: *  form:  A = U*B*VT, as computed by SGEBRD, then
036: * 
037: *     A = (U*Q) * S * (P**T*VT)
038: * 
039: *  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
040: *  for a given real 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) REAL array, dimension (LDVT, NCVT)
084: *          On entry, an N-by-NCVT matrix VT.
085: *          On exit, VT is overwritten by P**T * 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) REAL 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) REAL 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) REAL 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  REAL, 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:       REAL               ZERO
155:       PARAMETER          ( ZERO = 0.0E0 )
156:       REAL               ONE
157:       PARAMETER          ( ONE = 1.0E0 )
158:       REAL               NEGONE
159:       PARAMETER          ( NEGONE = -1.0E0 )
160:       REAL               HNDRTH
161:       PARAMETER          ( HNDRTH = 0.01E0 )
162:       REAL               TEN
163:       PARAMETER          ( TEN = 10.0E0 )
164:       REAL               HNDRD
165:       PARAMETER          ( HNDRD = 100.0E0 )
166:       REAL               MEIGTH
167:       PARAMETER          ( MEIGTH = -0.125E0 )
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:       REAL               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:       REAL               SLAMCH
183:       EXTERNAL           LSAME, SLAMCH
184: *     ..
185: *     .. External Subroutines ..
186:       EXTERNAL           SLARTG, SLAS2, SLASQ1, SLASR, SLASV2, SROT,
187:      $                   SSCAL, SSWAP, XERBLA
188: *     ..
189: *     .. Intrinsic Functions ..
190:       INTRINSIC          ABS, MAX, MIN, REAL, 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( 'SBDSQR', -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 SLASQ1( 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 = SLAMCH( 'Epsilon' )
245:       UNFL = SLAMCH( '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 SLARTG( 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 SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
264:      $                  LDU )
265:          IF( NCC.GT.0 )
266:      $      CALL SLASR( '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( REAL( 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 SLASV2( 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 SROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
385:      $                 SINR )
386:          IF( NRU.GT.0 )
387:      $      CALL SROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
388:          IF( NCC.GT.0 )
389:      $      CALL SROT( 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 SLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
489:          ELSE
490:             SLL = ABS( D( M ) )
491:             CALL SLAS2( 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 SLARTG( D( I )*CS, E( I ), CS, SN, R )
518:                IF( I.GT.LL )
519:      $            E( I-1 ) = OLDSN*R
520:                CALL SLARTG( 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 SLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
534:      $                     WORK( N ), VT( LL, 1 ), LDVT )
535:             IF( NRU.GT.0 )
536:      $         CALL SLASR( '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 SLASR( '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 SLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
556:                IF( I.LT.M )
557:      $            E( I ) = OLDSN*R
558:                CALL SLARTG( 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 SLASR( '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 SLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
575:      $                     WORK( N ), U( 1, LL ), LDU )
576:             IF( NCC.GT.0 )
577:      $         CALL SLASR( '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 SLARTG( 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 SLARTG( 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 SLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
624:      $                     WORK( N ), VT( LL, 1 ), LDVT )
625:             IF( NRU.GT.0 )
626:      $         CALL SLASR( '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 SLASR( '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 SLARTG( 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 SLARTG( 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 SLASR( '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 SLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
680:      $                     WORK( N ), U( 1, LL ), LDU )
681:             IF( NCC.GT.0 )
682:      $         CALL SLASR( '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 SSCAL( 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 SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
728:      $                     LDVT )
729:             IF( NRU.GT.0 )
730:      $         CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
731:             IF( NCC.GT.0 )
732:      $         CALL SSWAP( 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 SBDSQR
749: *
750:       END
751: