LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dbdsqr ( character  UPLO,
integer  N,
integer  NCVT,
integer  NRU,
integer  NCC,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( * )  WORK,
integer  INFO 
)

DBDSQR

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

Purpose:
 DBDSQR 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 DGEBRD, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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  DOUBLE PRECISION, 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 232 of file dbdsqr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: