LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sbdsqr()

subroutine sbdsqr ( character  UPLO,
integer  N,
integer  NCVT,
integer  NRU,
integer  NCC,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldvt, * )  VT,
integer  LDVT,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( * )  WORK,
integer  INFO 
)

SBDSQR

Download SBDSQR + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SBDSQR computes the singular values and, optionally, the right and/or
 left singular vectors from the singular value decomposition (SVD) of
 a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
 zero-shift QR algorithm.  The SVD of B has the form

    B = Q * S * P**T

 where S is the diagonal matrix of singular values, Q is an orthogonal
 matrix of left singular vectors, and P is an orthogonal matrix of
 right singular vectors.  If left singular vectors are requested, this
 subroutine actually returns U*Q instead of Q, and, if right singular
 vectors are requested, this subroutine returns P**T*VT instead of
 P**T, for given real input matrices U and VT.  When U and VT are the
 orthogonal matrices that reduce a general matrix A to bidiagonal
 form:  A = U*B*VT, as computed by SGEBRD, then

    A = (U*Q) * S * (P**T*VT)

 is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
 for a given real input matrix C.

 See "Computing  Small Singular Values of Bidiagonal Matrices With
 Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
 LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
 no. 5, pp. 873-912, Sept 1990) and
 "Accurate singular values and differential qd algorithms," by
 B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
 Department, University of California at Berkeley, July 1992
 for a detailed description of the algorithm.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  B is upper bidiagonal;
          = 'L':  B is lower bidiagonal.
[in]N
          N is INTEGER
          The order of the matrix B.  N >= 0.
[in]NCVT
          NCVT is INTEGER
          The number of columns of the matrix VT. NCVT >= 0.
[in]NRU
          NRU is INTEGER
          The number of rows of the matrix U. NRU >= 0.
[in]NCC
          NCC is INTEGER
          The number of columns of the matrix C. NCC >= 0.
[in,out]D
          D is REAL array, dimension (N)
          On entry, the n diagonal elements of the bidiagonal matrix B.
          On exit, if INFO=0, the singular values of B in decreasing
          order.
[in,out]E
          E is REAL array, dimension (N-1)
          On entry, the N-1 offdiagonal elements of the bidiagonal
          matrix B.
          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
          will contain the diagonal and superdiagonal elements of a
          bidiagonal matrix orthogonally equivalent to the one given
          as input.
[in,out]VT
          VT is REAL array, dimension (LDVT, NCVT)
          On entry, an N-by-NCVT matrix VT.
          On exit, VT is overwritten by P**T * VT.
          Not referenced if NCVT = 0.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.
          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
[in,out]U
          U is REAL array, dimension (LDU, N)
          On entry, an NRU-by-N matrix U.
          On exit, U is overwritten by U * Q.
          Not referenced if NRU = 0.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= max(1,NRU).
[in,out]C
          C is REAL array, dimension (LDC, NCC)
          On entry, an N-by-NCC matrix C.
          On exit, C is overwritten by Q**T * C.
          Not referenced if NCC = 0.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C.
          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
[out]WORK
          WORK is REAL array, dimension (4*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  If INFO = -i, the i-th argument had an illegal value
          > 0:
             if NCVT = NRU = NCC = 0,
                = 1, a split was marked by a positive value in E
                = 2, current block of Z not diagonalized after 30*N
                     iterations (in inner while loop)
                = 3, termination criterion of outer while loop not met
                     (program created more than N unreduced blocks)
             else NCVT = NRU = NCC = 0,
                   the algorithm did not converge; D and E contain the
                   elements of a bidiagonal matrix which is orthogonally
                   similar to the input matrix B;  if INFO = i, i
                   elements of E have not converged to zero.
Internal Parameters:
  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
          TOLMUL controls the convergence criterion of the QR loop.
          If it is positive, TOLMUL*EPS is the desired relative
             precision in the computed singular values.
          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
             desired absolute accuracy in the computed singular
             values (corresponds to relative accuracy
             abs(TOLMUL*EPS) in the largest singular value.
          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
             between 10 (for fast convergence) and .1/EPS
             (for there to be some accuracy in the results).
          Default is to lose at either one eighth or 2 of the
             available decimal digits in each computed singular value
             (whichever is smaller).

  MAXITR  INTEGER, default = 6
          MAXITR controls the maximum number of passes of the
          algorithm through its inner loop. The algorithms stops
          (and so fails to converge) if the number of passes
          through the inner loop exceeds MAXITR*N**2.
Note:
  Bug report from Cezary Dendek.
  On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
  removed since it can overflow pretty easily (for N larger or equal
  than 18,919). We instead use MAXITDIVN = MAXITR*N.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017

Definition at line 242 of file sbdsqr.f.

242 *
243 * -- LAPACK computational routine (version 3.7.1) --
244 * -- LAPACK is a software package provided by Univ. of Tennessee, --
245 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246 * June 2017
247 *
248 * .. Scalar Arguments ..
249  CHARACTER uplo
250  INTEGER info, ldc, ldu, ldvt, n, ncc, ncvt, nru
251 * ..
252 * .. Array Arguments ..
253  REAL c( ldc, * ), d( * ), e( * ), u( ldu, * ),
254  $ vt( ldvt, * ), work( * )
255 * ..
256 *
257 * =====================================================================
258 *
259 * .. Parameters ..
260  REAL zero
261  parameter( zero = 0.0e0 )
262  REAL one
263  parameter( one = 1.0e0 )
264  REAL negone
265  parameter( negone = -1.0e0 )
266  REAL hndrth
267  parameter( hndrth = 0.01e0 )
268  REAL ten
269  parameter( ten = 10.0e0 )
270  REAL hndrd
271  parameter( hndrd = 100.0e0 )
272  REAL meigth
273  parameter( meigth = -0.125e0 )
274  INTEGER maxitr
275  parameter( maxitr = 6 )
276 * ..
277 * .. Local Scalars ..
278  LOGICAL lower, rotate
279  INTEGER i, idir, isub, iter, iterdivn, j, ll, lll, m,
280  $ maxitdivn, nm1, nm12, nm13, oldll, oldm
281  REAL abse, abss, cosl, cosr, cs, eps, f, g, h, mu,
282  $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
283  $ sinr, sll, smax, smin, sminl, sminoa,
284  $ sn, thresh, tol, tolmul, unfl
285 * ..
286 * .. External Functions ..
287  LOGICAL lsame
288  REAL slamch
289  EXTERNAL lsame, slamch
290 * ..
291 * .. External Subroutines ..
292  EXTERNAL slartg, slas2, slasq1, slasr, slasv2, srot,
293  $ sscal, sswap, xerbla
294 * ..
295 * .. Intrinsic Functions ..
296  INTRINSIC abs, max, min, REAL, sign, sqrt
297 * ..
298 * .. Executable Statements ..
299 *
300 * Test the input parameters.
301 *
302  info = 0
303  lower = lsame( uplo, 'L' )
304  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
305  info = -1
306  ELSE IF( n.LT.0 ) THEN
307  info = -2
308  ELSE IF( ncvt.LT.0 ) THEN
309  info = -3
310  ELSE IF( nru.LT.0 ) THEN
311  info = -4
312  ELSE IF( ncc.LT.0 ) THEN
313  info = -5
314  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
315  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
316  info = -9
317  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
318  info = -11
319  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
320  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
321  info = -13
322  END IF
323  IF( info.NE.0 ) THEN
324  CALL xerbla( 'SBDSQR', -info )
325  RETURN
326  END IF
327  IF( n.EQ.0 )
328  $ RETURN
329  IF( n.EQ.1 )
330  $ GO TO 160
331 *
332 * ROTATE is true if any singular vectors desired, false otherwise
333 *
334  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
335 *
336 * If no singular vectors desired, use qd algorithm
337 *
338  IF( .NOT.rotate ) THEN
339  CALL slasq1( n, d, e, work, info )
340 *
341 * If INFO equals 2, dqds didn't finish, try to finish
342 *
343  IF( info .NE. 2 ) RETURN
344  info = 0
345  END IF
346 *
347  nm1 = n - 1
348  nm12 = nm1 + nm1
349  nm13 = nm12 + nm1
350  idir = 0
351 *
352 * Get machine constants
353 *
354  eps = slamch( 'Epsilon' )
355  unfl = slamch( 'Safe minimum' )
356 *
357 * If matrix lower bidiagonal, rotate to be upper bidiagonal
358 * by applying Givens rotations on the left
359 *
360  IF( lower ) THEN
361  DO 10 i = 1, n - 1
362  CALL slartg( d( i ), e( i ), cs, sn, r )
363  d( i ) = r
364  e( i ) = sn*d( i+1 )
365  d( i+1 ) = cs*d( i+1 )
366  work( i ) = cs
367  work( nm1+i ) = sn
368  10 CONTINUE
369 *
370 * Update singular vectors if desired
371 *
372  IF( nru.GT.0 )
373  $ CALL slasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
374  $ ldu )
375  IF( ncc.GT.0 )
376  $ CALL slasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), c,
377  $ ldc )
378  END IF
379 *
380 * Compute singular values to relative accuracy TOL
381 * (By setting TOL to be negative, algorithm will compute
382 * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
383 *
384  tolmul = max( ten, min( hndrd, eps**meigth ) )
385  tol = tolmul*eps
386 *
387 * Compute approximate maximum, minimum singular values
388 *
389  smax = zero
390  DO 20 i = 1, n
391  smax = max( smax, abs( d( i ) ) )
392  20 CONTINUE
393  DO 30 i = 1, n - 1
394  smax = max( smax, abs( e( i ) ) )
395  30 CONTINUE
396  sminl = zero
397  IF( tol.GE.zero ) THEN
398 *
399 * Relative accuracy desired
400 *
401  sminoa = abs( d( 1 ) )
402  IF( sminoa.EQ.zero )
403  $ GO TO 50
404  mu = sminoa
405  DO 40 i = 2, n
406  mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
407  sminoa = min( sminoa, mu )
408  IF( sminoa.EQ.zero )
409  $ GO TO 50
410  40 CONTINUE
411  50 CONTINUE
412  sminoa = sminoa / sqrt( REAL( N ) )
413  thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
414  ELSE
415 *
416 * Absolute accuracy desired
417 *
418  thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
419  END IF
420 *
421 * Prepare for main iteration loop for the singular values
422 * (MAXIT is the maximum number of passes through the inner
423 * loop permitted before nonconvergence signalled.)
424 *
425  maxitdivn = maxitr*n
426  iterdivn = 0
427  iter = -1
428  oldll = -1
429  oldm = -1
430 *
431 * M points to last element of unconverged part of matrix
432 *
433  m = n
434 *
435 * Begin main iteration loop
436 *
437  60 CONTINUE
438 *
439 * Check for convergence or exceeding iteration count
440 *
441  IF( m.LE.1 )
442  $ GO TO 160
443 *
444  IF( iter.GE.n ) THEN
445  iter = iter - n
446  iterdivn = iterdivn + 1
447  IF( iterdivn.GE.maxitdivn )
448  $ GO TO 200
449  END IF
450 *
451 * Find diagonal block of matrix to work on
452 *
453  IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
454  $ d( m ) = zero
455  smax = abs( d( m ) )
456  smin = smax
457  DO 70 lll = 1, m - 1
458  ll = m - lll
459  abss = abs( d( ll ) )
460  abse = abs( e( ll ) )
461  IF( tol.LT.zero .AND. abss.LE.thresh )
462  $ d( ll ) = zero
463  IF( abse.LE.thresh )
464  $ GO TO 80
465  smin = min( smin, abss )
466  smax = max( smax, abss, abse )
467  70 CONTINUE
468  ll = 0
469  GO TO 90
470  80 CONTINUE
471  e( ll ) = zero
472 *
473 * Matrix splits since E(LL) = 0
474 *
475  IF( ll.EQ.m-1 ) THEN
476 *
477 * Convergence of bottom singular value, return to top of loop
478 *
479  m = m - 1
480  GO TO 60
481  END IF
482  90 CONTINUE
483  ll = ll + 1
484 *
485 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
486 *
487  IF( ll.EQ.m-1 ) THEN
488 *
489 * 2 by 2 block, handle separately
490 *
491  CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
492  $ cosr, sinl, cosl )
493  d( m-1 ) = sigmx
494  e( m-1 ) = zero
495  d( m ) = sigmn
496 *
497 * Compute singular vectors, if desired
498 *
499  IF( ncvt.GT.0 )
500  $ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
501  $ sinr )
502  IF( nru.GT.0 )
503  $ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
504  IF( ncc.GT.0 )
505  $ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
506  $ sinl )
507  m = m - 2
508  GO TO 60
509  END IF
510 *
511 * If working on new submatrix, choose shift direction
512 * (from larger end diagonal element towards smaller)
513 *
514  IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
515  IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
516 *
517 * Chase bulge from top (big end) to bottom (small end)
518 *
519  idir = 1
520  ELSE
521 *
522 * Chase bulge from bottom (big end) to top (small end)
523 *
524  idir = 2
525  END IF
526  END IF
527 *
528 * Apply convergence tests
529 *
530  IF( idir.EQ.1 ) THEN
531 *
532 * Run convergence test in forward direction
533 * First apply standard test to bottom of matrix
534 *
535  IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
536  $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
537  e( m-1 ) = zero
538  GO TO 60
539  END IF
540 *
541  IF( tol.GE.zero ) THEN
542 *
543 * If relative accuracy desired,
544 * apply convergence criterion forward
545 *
546  mu = abs( d( ll ) )
547  sminl = mu
548  DO 100 lll = ll, m - 1
549  IF( abs( e( lll ) ).LE.tol*mu ) THEN
550  e( lll ) = zero
551  GO TO 60
552  END IF
553  mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
554  sminl = min( sminl, mu )
555  100 CONTINUE
556  END IF
557 *
558  ELSE
559 *
560 * Run convergence test in backward direction
561 * First apply standard test to top of matrix
562 *
563  IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
564  $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
565  e( ll ) = zero
566  GO TO 60
567  END IF
568 *
569  IF( tol.GE.zero ) THEN
570 *
571 * If relative accuracy desired,
572 * apply convergence criterion backward
573 *
574  mu = abs( d( m ) )
575  sminl = mu
576  DO 110 lll = m - 1, ll, -1
577  IF( abs( e( lll ) ).LE.tol*mu ) THEN
578  e( lll ) = zero
579  GO TO 60
580  END IF
581  mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
582  sminl = min( sminl, mu )
583  110 CONTINUE
584  END IF
585  END IF
586  oldll = ll
587  oldm = m
588 *
589 * Compute shift. First, test if shifting would ruin relative
590 * accuracy, and if so set the shift to zero.
591 *
592  IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
593  $ max( eps, hndrth*tol ) ) THEN
594 *
595 * Use a zero shift to avoid loss of relative accuracy
596 *
597  shift = zero
598  ELSE
599 *
600 * Compute the shift from 2-by-2 block at end of matrix
601 *
602  IF( idir.EQ.1 ) THEN
603  sll = abs( d( ll ) )
604  CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
605  ELSE
606  sll = abs( d( m ) )
607  CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
608  END IF
609 *
610 * Test if shift negligible, and if so set to zero
611 *
612  IF( sll.GT.zero ) THEN
613  IF( ( shift / sll )**2.LT.eps )
614  $ shift = zero
615  END IF
616  END IF
617 *
618 * Increment iteration count
619 *
620  iter = iter + m - ll
621 *
622 * If SHIFT = 0, do simplified QR iteration
623 *
624  IF( shift.EQ.zero ) THEN
625  IF( idir.EQ.1 ) THEN
626 *
627 * Chase bulge from top to bottom
628 * Save cosines and sines for later singular vector updates
629 *
630  cs = one
631  oldcs = one
632  DO 120 i = ll, m - 1
633  CALL slartg( d( i )*cs, e( i ), cs, sn, r )
634  IF( i.GT.ll )
635  $ e( i-1 ) = oldsn*r
636  CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
637  work( i-ll+1 ) = cs
638  work( i-ll+1+nm1 ) = sn
639  work( i-ll+1+nm12 ) = oldcs
640  work( i-ll+1+nm13 ) = oldsn
641  120 CONTINUE
642  h = d( m )*cs
643  d( m ) = h*oldcs
644  e( m-1 ) = h*oldsn
645 *
646 * Update singular vectors
647 *
648  IF( ncvt.GT.0 )
649  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
650  $ work( n ), vt( ll, 1 ), ldvt )
651  IF( nru.GT.0 )
652  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
653  $ work( nm13+1 ), u( 1, ll ), ldu )
654  IF( ncc.GT.0 )
655  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
656  $ work( nm13+1 ), c( ll, 1 ), ldc )
657 *
658 * Test convergence
659 *
660  IF( abs( e( m-1 ) ).LE.thresh )
661  $ e( m-1 ) = zero
662 *
663  ELSE
664 *
665 * Chase bulge from bottom to top
666 * Save cosines and sines for later singular vector updates
667 *
668  cs = one
669  oldcs = one
670  DO 130 i = m, ll + 1, -1
671  CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
672  IF( i.LT.m )
673  $ e( i ) = oldsn*r
674  CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
675  work( i-ll ) = cs
676  work( i-ll+nm1 ) = -sn
677  work( i-ll+nm12 ) = oldcs
678  work( i-ll+nm13 ) = -oldsn
679  130 CONTINUE
680  h = d( ll )*cs
681  d( ll ) = h*oldcs
682  e( ll ) = h*oldsn
683 *
684 * Update singular vectors
685 *
686  IF( ncvt.GT.0 )
687  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
688  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
689  IF( nru.GT.0 )
690  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
691  $ work( n ), u( 1, ll ), ldu )
692  IF( ncc.GT.0 )
693  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
694  $ work( n ), c( ll, 1 ), ldc )
695 *
696 * Test convergence
697 *
698  IF( abs( e( ll ) ).LE.thresh )
699  $ e( ll ) = zero
700  END IF
701  ELSE
702 *
703 * Use nonzero shift
704 *
705  IF( idir.EQ.1 ) THEN
706 *
707 * Chase bulge from top to bottom
708 * Save cosines and sines for later singular vector updates
709 *
710  f = ( abs( d( ll ) )-shift )*
711  $ ( sign( one, d( ll ) )+shift / d( ll ) )
712  g = e( ll )
713  DO 140 i = ll, m - 1
714  CALL slartg( f, g, cosr, sinr, r )
715  IF( i.GT.ll )
716  $ e( i-1 ) = r
717  f = cosr*d( i ) + sinr*e( i )
718  e( i ) = cosr*e( i ) - sinr*d( i )
719  g = sinr*d( i+1 )
720  d( i+1 ) = cosr*d( i+1 )
721  CALL slartg( f, g, cosl, sinl, r )
722  d( i ) = r
723  f = cosl*e( i ) + sinl*d( i+1 )
724  d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
725  IF( i.LT.m-1 ) THEN
726  g = sinl*e( i+1 )
727  e( i+1 ) = cosl*e( i+1 )
728  END IF
729  work( i-ll+1 ) = cosr
730  work( i-ll+1+nm1 ) = sinr
731  work( i-ll+1+nm12 ) = cosl
732  work( i-ll+1+nm13 ) = sinl
733  140 CONTINUE
734  e( m-1 ) = f
735 *
736 * Update singular vectors
737 *
738  IF( ncvt.GT.0 )
739  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
740  $ work( n ), vt( ll, 1 ), ldvt )
741  IF( nru.GT.0 )
742  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
743  $ work( nm13+1 ), u( 1, ll ), ldu )
744  IF( ncc.GT.0 )
745  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
746  $ work( nm13+1 ), c( ll, 1 ), ldc )
747 *
748 * Test convergence
749 *
750  IF( abs( e( m-1 ) ).LE.thresh )
751  $ e( m-1 ) = zero
752 *
753  ELSE
754 *
755 * Chase bulge from bottom to top
756 * Save cosines and sines for later singular vector updates
757 *
758  f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
759  $ d( m ) )
760  g = e( m-1 )
761  DO 150 i = m, ll + 1, -1
762  CALL slartg( f, g, cosr, sinr, r )
763  IF( i.LT.m )
764  $ e( i ) = r
765  f = cosr*d( i ) + sinr*e( i-1 )
766  e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
767  g = sinr*d( i-1 )
768  d( i-1 ) = cosr*d( i-1 )
769  CALL slartg( f, g, cosl, sinl, r )
770  d( i ) = r
771  f = cosl*e( i-1 ) + sinl*d( i-1 )
772  d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
773  IF( i.GT.ll+1 ) THEN
774  g = sinl*e( i-2 )
775  e( i-2 ) = cosl*e( i-2 )
776  END IF
777  work( i-ll ) = cosr
778  work( i-ll+nm1 ) = -sinr
779  work( i-ll+nm12 ) = cosl
780  work( i-ll+nm13 ) = -sinl
781  150 CONTINUE
782  e( ll ) = f
783 *
784 * Test convergence
785 *
786  IF( abs( e( ll ) ).LE.thresh )
787  $ e( ll ) = zero
788 *
789 * Update singular vectors if desired
790 *
791  IF( ncvt.GT.0 )
792  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
793  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
794  IF( nru.GT.0 )
795  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
796  $ work( n ), u( 1, ll ), ldu )
797  IF( ncc.GT.0 )
798  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
799  $ work( n ), c( ll, 1 ), ldc )
800  END IF
801  END IF
802 *
803 * QR iteration finished, go back and check convergence
804 *
805  GO TO 60
806 *
807 * All singular values converged, so make them positive
808 *
809  160 CONTINUE
810  DO 170 i = 1, n
811  IF( d( i ).LT.zero ) THEN
812  d( i ) = -d( i )
813 *
814 * Change sign of singular vectors, if desired
815 *
816  IF( ncvt.GT.0 )
817  $ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
818  END IF
819  170 CONTINUE
820 *
821 * Sort the singular values into decreasing order (insertion sort on
822 * singular values, but only one transposition per singular vector)
823 *
824  DO 190 i = 1, n - 1
825 *
826 * Scan for smallest D(I)
827 *
828  isub = 1
829  smin = d( 1 )
830  DO 180 j = 2, n + 1 - i
831  IF( d( j ).LE.smin ) THEN
832  isub = j
833  smin = d( j )
834  END IF
835  180 CONTINUE
836  IF( isub.NE.n+1-i ) THEN
837 *
838 * Swap singular values and vectors
839 *
840  d( isub ) = d( n+1-i )
841  d( n+1-i ) = smin
842  IF( ncvt.GT.0 )
843  $ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
844  $ ldvt )
845  IF( nru.GT.0 )
846  $ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
847  IF( ncc.GT.0 )
848  $ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
849  END IF
850  190 CONTINUE
851  GO TO 220
852 *
853 * Maximum number of iterations exceeded, failure to converge
854 *
855  200 CONTINUE
856  info = 0
857  DO 210 i = 1, n - 1
858  IF( e( i ).NE.zero )
859  $ info = info + 1
860  210 CONTINUE
861  220 CONTINUE
862  RETURN
863 *
864 * End of SBDSQR
865 *
subroutine slasq1(N, D, E, WORK, INFO)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition: slasq1.f:110
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:109
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:94
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: slasr.f:201
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: slasv2.f:140
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
Here is the call graph for this function:
Here is the caller graph for this function: