LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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:111
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: