LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dbdsqr()

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.
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 243 of file dbdsqr.f.

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