LAPACK 3.12.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, 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 smin = 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 DO 70 lll = 1, m - 1
455 ll = m - lll
456 abss = abs( d( ll ) )
457 abse = abs( e( ll ) )
458 IF( tol.LT.zero .AND. abss.LE.thresh )
459 $ d( ll ) = zero
460 IF( abse.LE.thresh )
461 $ GO TO 80
462 smax = max( smax, abss, abse )
463 70 CONTINUE
464 ll = 0
465 GO TO 90
466 80 CONTINUE
467 e( ll ) = zero
468*
469* Matrix splits since E(LL) = 0
470*
471 IF( ll.EQ.m-1 ) THEN
472*
473* Convergence of bottom singular value, return to top of loop
474*
475 m = m - 1
476 GO TO 60
477 END IF
478 90 CONTINUE
479 ll = ll + 1
480*
481* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
482*
483 IF( ll.EQ.m-1 ) THEN
484*
485* 2 by 2 block, handle separately
486*
487 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
488 $ cosr, sinl, cosl )
489 d( m-1 ) = sigmx
490 e( m-1 ) = zero
491 d( m ) = sigmn
492*
493* Compute singular vectors, if desired
494*
495 IF( ncvt.GT.0 )
496 $ CALL drot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
497 $ sinr )
498 IF( nru.GT.0 )
499 $ CALL drot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
500 IF( ncc.GT.0 )
501 $ CALL drot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
502 $ sinl )
503 m = m - 2
504 GO TO 60
505 END IF
506*
507* If working on new submatrix, choose shift direction
508* (from larger end diagonal element towards smaller)
509*
510 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
511 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
512*
513* Chase bulge from top (big end) to bottom (small end)
514*
515 idir = 1
516 ELSE
517*
518* Chase bulge from bottom (big end) to top (small end)
519*
520 idir = 2
521 END IF
522 END IF
523*
524* Apply convergence tests
525*
526 IF( idir.EQ.1 ) THEN
527*
528* Run convergence test in forward direction
529* First apply standard test to bottom of matrix
530*
531 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
532 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
533 e( m-1 ) = zero
534 GO TO 60
535 END IF
536*
537 IF( tol.GE.zero ) THEN
538*
539* If relative accuracy desired,
540* apply convergence criterion forward
541*
542 mu = abs( d( ll ) )
543 smin = mu
544 DO 100 lll = ll, m - 1
545 IF( abs( e( lll ) ).LE.tol*mu ) THEN
546 e( lll ) = zero
547 GO TO 60
548 END IF
549 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
550 smin = min( smin, mu )
551 100 CONTINUE
552 END IF
553*
554 ELSE
555*
556* Run convergence test in backward direction
557* First apply standard test to top of matrix
558*
559 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
560 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
561 e( ll ) = zero
562 GO TO 60
563 END IF
564*
565 IF( tol.GE.zero ) THEN
566*
567* If relative accuracy desired,
568* apply convergence criterion backward
569*
570 mu = abs( d( m ) )
571 smin = mu
572 DO 110 lll = m - 1, ll, -1
573 IF( abs( e( lll ) ).LE.tol*mu ) THEN
574 e( lll ) = zero
575 GO TO 60
576 END IF
577 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
578 smin = min( smin, mu )
579 110 CONTINUE
580 END IF
581 END IF
582 oldll = ll
583 oldm = m
584*
585* Compute shift. First, test if shifting would ruin relative
586* accuracy, and if so set the shift to zero.
587*
588 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
589 $ max( eps, hndrth*tol ) ) THEN
590*
591* Use a zero shift to avoid loss of relative accuracy
592*
593 shift = zero
594 ELSE
595*
596* Compute the shift from 2-by-2 block at end of matrix
597*
598 IF( idir.EQ.1 ) THEN
599 sll = abs( d( ll ) )
600 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
601 ELSE
602 sll = abs( d( m ) )
603 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
604 END IF
605*
606* Test if shift negligible, and if so set to zero
607*
608 IF( sll.GT.zero ) THEN
609 IF( ( shift / sll )**2.LT.eps )
610 $ shift = zero
611 END IF
612 END IF
613*
614* Increment iteration count
615*
616 iter = iter + m - ll
617*
618* If SHIFT = 0, do simplified QR iteration
619*
620 IF( shift.EQ.zero ) THEN
621 IF( idir.EQ.1 ) THEN
622*
623* Chase bulge from top to bottom
624* Save cosines and sines for later singular vector updates
625*
626 cs = one
627 oldcs = one
628 DO 120 i = ll, m - 1
629 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
630 IF( i.GT.ll )
631 $ e( i-1 ) = oldsn*r
632 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
633 work( i-ll+1 ) = cs
634 work( i-ll+1+nm1 ) = sn
635 work( i-ll+1+nm12 ) = oldcs
636 work( i-ll+1+nm13 ) = oldsn
637 120 CONTINUE
638 h = d( m )*cs
639 d( m ) = h*oldcs
640 e( m-1 ) = h*oldsn
641*
642* Update singular vectors
643*
644 IF( ncvt.GT.0 )
645 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
646 $ work( n ), vt( ll, 1 ), ldvt )
647 IF( nru.GT.0 )
648 $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
649 $ work( nm13+1 ), u( 1, ll ), ldu )
650 IF( ncc.GT.0 )
651 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
652 $ work( nm13+1 ), c( ll, 1 ), ldc )
653*
654* Test convergence
655*
656 IF( abs( e( m-1 ) ).LE.thresh )
657 $ e( m-1 ) = zero
658*
659 ELSE
660*
661* Chase bulge from bottom to top
662* Save cosines and sines for later singular vector updates
663*
664 cs = one
665 oldcs = one
666 DO 130 i = m, ll + 1, -1
667 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
668 IF( i.LT.m )
669 $ e( i ) = oldsn*r
670 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
671 work( i-ll ) = cs
672 work( i-ll+nm1 ) = -sn
673 work( i-ll+nm12 ) = oldcs
674 work( i-ll+nm13 ) = -oldsn
675 130 CONTINUE
676 h = d( ll )*cs
677 d( ll ) = h*oldcs
678 e( ll ) = h*oldsn
679*
680* Update singular vectors
681*
682 IF( ncvt.GT.0 )
683 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
684 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
685 IF( nru.GT.0 )
686 $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
687 $ work( n ), u( 1, ll ), ldu )
688 IF( ncc.GT.0 )
689 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
690 $ work( n ), c( ll, 1 ), ldc )
691*
692* Test convergence
693*
694 IF( abs( e( ll ) ).LE.thresh )
695 $ e( ll ) = zero
696 END IF
697 ELSE
698*
699* Use nonzero shift
700*
701 IF( idir.EQ.1 ) THEN
702*
703* Chase bulge from top to bottom
704* Save cosines and sines for later singular vector updates
705*
706 f = ( abs( d( ll ) )-shift )*
707 $ ( sign( one, d( ll ) )+shift / d( ll ) )
708 g = e( ll )
709 DO 140 i = ll, m - 1
710 CALL dlartg( f, g, cosr, sinr, r )
711 IF( i.GT.ll )
712 $ e( i-1 ) = r
713 f = cosr*d( i ) + sinr*e( i )
714 e( i ) = cosr*e( i ) - sinr*d( i )
715 g = sinr*d( i+1 )
716 d( i+1 ) = cosr*d( i+1 )
717 CALL dlartg( f, g, cosl, sinl, r )
718 d( i ) = r
719 f = cosl*e( i ) + sinl*d( i+1 )
720 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
721 IF( i.LT.m-1 ) THEN
722 g = sinl*e( i+1 )
723 e( i+1 ) = cosl*e( i+1 )
724 END IF
725 work( i-ll+1 ) = cosr
726 work( i-ll+1+nm1 ) = sinr
727 work( i-ll+1+nm12 ) = cosl
728 work( i-ll+1+nm13 ) = sinl
729 140 CONTINUE
730 e( m-1 ) = f
731*
732* Update singular vectors
733*
734 IF( ncvt.GT.0 )
735 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
736 $ work( n ), vt( ll, 1 ), ldvt )
737 IF( nru.GT.0 )
738 $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
739 $ work( nm13+1 ), u( 1, ll ), ldu )
740 IF( ncc.GT.0 )
741 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
742 $ work( nm13+1 ), c( ll, 1 ), ldc )
743*
744* Test convergence
745*
746 IF( abs( e( m-1 ) ).LE.thresh )
747 $ e( m-1 ) = zero
748*
749 ELSE
750*
751* Chase bulge from bottom to top
752* Save cosines and sines for later singular vector updates
753*
754 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
755 $ d( m ) )
756 g = e( m-1 )
757 DO 150 i = m, ll + 1, -1
758 CALL dlartg( f, g, cosr, sinr, r )
759 IF( i.LT.m )
760 $ e( i ) = r
761 f = cosr*d( i ) + sinr*e( i-1 )
762 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
763 g = sinr*d( i-1 )
764 d( i-1 ) = cosr*d( i-1 )
765 CALL dlartg( f, g, cosl, sinl, r )
766 d( i ) = r
767 f = cosl*e( i-1 ) + sinl*d( i-1 )
768 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
769 IF( i.GT.ll+1 ) THEN
770 g = sinl*e( i-2 )
771 e( i-2 ) = cosl*e( i-2 )
772 END IF
773 work( i-ll ) = cosr
774 work( i-ll+nm1 ) = -sinr
775 work( i-ll+nm12 ) = cosl
776 work( i-ll+nm13 ) = -sinl
777 150 CONTINUE
778 e( ll ) = f
779*
780* Test convergence
781*
782 IF( abs( e( ll ) ).LE.thresh )
783 $ e( ll ) = zero
784*
785* Update singular vectors if desired
786*
787 IF( ncvt.GT.0 )
788 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
789 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
790 IF( nru.GT.0 )
791 $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
792 $ work( n ), u( 1, ll ), ldu )
793 IF( ncc.GT.0 )
794 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
795 $ work( n ), c( ll, 1 ), ldc )
796 END IF
797 END IF
798*
799* QR iteration finished, go back and check convergence
800*
801 GO TO 60
802*
803* All singular values converged, so make them positive
804*
805 160 CONTINUE
806 DO 170 i = 1, n
807 IF( d( i ).LT.zero ) THEN
808 d( i ) = -d( i )
809*
810* Change sign of singular vectors, if desired
811*
812 IF( ncvt.GT.0 )
813 $ CALL dscal( ncvt, negone, vt( i, 1 ), ldvt )
814 END IF
815 170 CONTINUE
816*
817* Sort the singular values into decreasing order (insertion sort on
818* singular values, but only one transposition per singular vector)
819*
820 DO 190 i = 1, n - 1
821*
822* Scan for smallest D(I)
823*
824 isub = 1
825 smin = d( 1 )
826 DO 180 j = 2, n + 1 - i
827 IF( d( j ).LE.smin ) THEN
828 isub = j
829 smin = d( j )
830 END IF
831 180 CONTINUE
832 IF( isub.NE.n+1-i ) THEN
833*
834* Swap singular values and vectors
835*
836 d( isub ) = d( n+1-i )
837 d( n+1-i ) = smin
838 IF( ncvt.GT.0 )
839 $ CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
840 $ ldvt )
841 IF( nru.GT.0 )
842 $ CALL dswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
843 IF( ncc.GT.0 )
844 $ CALL dswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
845 END IF
846 190 CONTINUE
847 GO TO 220
848*
849* Maximum number of iterations exceeded, failure to converge
850*
851 200 CONTINUE
852 info = 0
853 DO 210 i = 1, n - 1
854 IF( e( i ).NE.zero )
855 $ info = info + 1
856 210 CONTINUE
857 220 CONTINUE
858 RETURN
859*
860* End of DBDSQR
861*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:105
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 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:136
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: