LAPACK  3.10.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-1))
[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.

Definition at line 239 of file dbdsqr.f.

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