220 SUBROUTINE zbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
221 $ LDU, C, LDC, RWORK, INFO )
229 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
232 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
233 COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
239 DOUBLE PRECISION ZERO
240 parameter( zero = 0.0d0 )
242 parameter( one = 1.0d0 )
243 DOUBLE PRECISION NEGONE
244 parameter( negone = -1.0d0 )
245 DOUBLE PRECISION HNDRTH
246 parameter( hndrth = 0.01d0 )
248 parameter( ten = 10.0d0 )
249 DOUBLE PRECISION HNDRD
250 parameter( hndrd = 100.0d0 )
251 DOUBLE PRECISION MEIGTH
252 parameter( meigth = -0.125d0 )
254 parameter( maxitr = 6 )
257 LOGICAL LOWER, ROTATE
258 INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
259 $ nm12, nm13, oldll, oldm
260 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
261 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
262 $ sinr, sll, smax, smin, sminl, sminoa,
263 $ sn, thresh, tol, tolmul, unfl
267 DOUBLE PRECISION DLAMCH
268 EXTERNAL lsame, dlamch
275 INTRINSIC abs, dble, max, min, sign, sqrt
282 lower = lsame( uplo,
'L' )
283 IF( .NOT.lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
285 ELSE IF( n.LT.0 )
THEN
287 ELSE IF( ncvt.LT.0 )
THEN
289 ELSE IF( nru.LT.0 )
THEN
291 ELSE IF( ncc.LT.0 )
THEN
293 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
294 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) )
THEN
296 ELSE IF( ldu.LT.max( 1, nru ) )
THEN
298 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
299 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) )
THEN
303 CALL xerbla(
'ZBDSQR', -info )
313 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
317 IF( .NOT.rotate )
THEN
318 CALL dlasq1( n, d, e, rwork, info )
322 IF( info .NE. 2 )
RETURN
333 eps = dlamch(
'Epsilon' )
334 unfl = dlamch(
'Safe minimum' )
341 CALL dlartg( d( i ), e( i ), cs, sn, r )
344 d( i+1 ) = cs*d( i+1 )
352 $
CALL zlasr(
'R',
'V',
'F', nru, n, rwork( 1 ), rwork( n ),
355 $
CALL zlasr(
'L',
'V',
'F', n, ncc, rwork( 1 ), rwork( n ),
363 tolmul = max( ten, min( hndrd, eps**meigth ) )
370 smax = max( smax, abs( d( i ) ) )
373 smax = max( smax, abs( e( i ) ) )
376 IF( tol.GE.zero )
THEN
380 sminoa = abs( d( 1 ) )
385 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
386 sminoa = min( sminoa, mu )
391 sminoa = sminoa / sqrt( dble( n ) )
392 thresh = max( tol*sminoa, maxitr*n*n*unfl )
397 thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
426 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
432 abss = abs( d( ll ) )
433 abse = abs( e( ll ) )
434 IF( tol.LT.zero .AND. abss.LE.thresh )
438 smin = min( smin, abss )
439 smax = max( smax, abss, abse )
464 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
473 $
CALL zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
476 $
CALL zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
478 $
CALL zdrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
487 IF( ll.GT.oldm .OR. m.LT.oldll )
THEN
488 IF( abs( d( ll ) ).GE.abs( d( m ) ) )
THEN
508 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
509 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) )
THEN
514 IF( tol.GE.zero )
THEN
521 DO 100 lll = ll, m - 1
522 IF( abs( e( lll ) ).LE.tol*mu )
THEN
526 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
527 sminl = min( sminl, mu )
536 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
537 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) )
THEN
542 IF( tol.GE.zero )
THEN
549 DO 110 lll = m - 1, ll, -1
550 IF( abs( e( lll ) ).LE.tol*mu )
THEN
554 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
555 sminl = min( sminl, mu )
565 IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
566 $ max( eps, hndrth*tol ) )
THEN
577 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
580 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
585 IF( sll.GT.zero )
THEN
586 IF( ( shift / sll )**2.LT.eps )
597 IF( shift.EQ.zero )
THEN
606 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
609 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
611 rwork( i-ll+1+nm1 ) = sn
612 rwork( i-ll+1+nm12 ) = oldcs
613 rwork( i-ll+1+nm13 ) = oldsn
622 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
623 $ rwork( n ), vt( ll, 1 ), ldvt )
625 $
CALL zlasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
626 $ rwork( nm13+1 ), u( 1, ll ), ldu )
628 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncc, rwork( nm12+1 ),
629 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
633 IF( abs( e( m-1 ) ).LE.thresh )
643 DO 130 i = m, ll + 1, -1
644 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
647 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
649 rwork( i-ll+nm1 ) = -sn
650 rwork( i-ll+nm12 ) = oldcs
651 rwork( i-ll+nm13 ) = -oldsn
660 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
661 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
663 $
CALL zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
664 $ rwork( n ), u( 1, ll ), ldu )
666 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
667 $ rwork( n ), c( ll, 1 ), ldc )
671 IF( abs( e( ll ) ).LE.thresh )
683 f = ( abs( d( ll ) )-shift )*
684 $ ( sign( one, d( ll ) )+shift / d( ll ) )
687 CALL dlartg( f, g, cosr, sinr, r )
690 f = cosr*d( i ) + sinr*e( i )
691 e( i ) = cosr*e( i ) - sinr*d( i )
693 d( i+1 ) = cosr*d( i+1 )
694 CALL dlartg( f, g, cosl, sinl, r )
696 f = cosl*e( i ) + sinl*d( i+1 )
697 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
700 e( i+1 ) = cosl*e( i+1 )
702 rwork( i-ll+1 ) = cosr
703 rwork( i-ll+1+nm1 ) = sinr
704 rwork( i-ll+1+nm12 ) = cosl
705 rwork( i-ll+1+nm13 ) = sinl
712 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
713 $ rwork( n ), vt( ll, 1 ), ldvt )
715 $
CALL zlasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
716 $ rwork( nm13+1 ), u( 1, ll ), ldu )
718 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncc, rwork( nm12+1 ),
719 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
723 IF( abs( e( m-1 ) ).LE.thresh )
731 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
734 DO 150 i = m, ll + 1, -1
735 CALL dlartg( f, g, cosr, sinr, r )
738 f = cosr*d( i ) + sinr*e( i-1 )
739 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
741 d( i-1 ) = cosr*d( i-1 )
742 CALL dlartg( f, g, cosl, sinl, r )
744 f = cosl*e( i-1 ) + sinl*d( i-1 )
745 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
748 e( i-2 ) = cosl*e( i-2 )
751 rwork( i-ll+nm1 ) = -sinr
752 rwork( i-ll+nm12 ) = cosl
753 rwork( i-ll+nm13 ) = -sinl
759 IF( abs( e( ll ) ).LE.thresh )
765 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
766 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
768 $
CALL zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
769 $ rwork( n ), u( 1, ll ), ldu )
771 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
772 $ rwork( n ), c( ll, 1 ), ldc )
784 IF( d( i ).LT.zero )
THEN
790 $
CALL zdscal( ncvt, negone, vt( i, 1 ), ldvt )
803 DO 180 j = 2, n + 1 - i
804 IF( d( j ).LE.smin )
THEN
809 IF( isub.NE.n+1-i )
THEN
813 d( isub ) = d( n+1-i )
816 $
CALL zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
819 $
CALL zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
821 $
CALL zswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlasq1(N, D, E, WORK, INFO)
DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zdrot(N, ZX, INCX, ZY, INCY, C, S)
ZDROT
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR