223 SUBROUTINE zbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
224 $ ldu, c, ldc, rwork, info )
233 INTEGER info, ldc, ldu, ldvt, n, ncc, ncvt, nru
236 DOUBLE PRECISION d( * ), e( * ), rwork( * )
237 COMPLEX*16 c( ldc, * ), u( ldu, * ), vt( ldvt, * )
243 DOUBLE PRECISION zero
244 parameter( zero = 0.0d0 )
246 parameter( one = 1.0d0 )
247 DOUBLE PRECISION negone
248 parameter( negone = -1.0d0 )
249 DOUBLE PRECISION hndrth
250 parameter( hndrth = 0.01d0 )
252 parameter( ten = 10.0d0 )
253 DOUBLE PRECISION hndrd
254 parameter( hndrd = 100.0d0 )
255 DOUBLE PRECISION meigth
256 parameter( meigth = -0.125d0 )
258 parameter( maxitr = 6 )
261 LOGICAL lower, rotate
262 INTEGER i, idir, isub, iter,
j, ll, lll, m, maxit, nm1,
263 $ nm12, nm13, oldll, oldm
264 DOUBLE PRECISION abse, abss, cosl, cosr, cs, eps, f, g, h, mu,
265 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
266 $ sinr, sll, smax, smin, sminl, sminoa,
267 $ sn, thresh, tol, tolmul, unfl
279 INTRINSIC abs, dble, max, min, sign, sqrt
286 lower =
lsame( uplo,
'L' )
287 IF( .NOT.
lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
289 ELSE IF( n.LT.0 )
THEN
291 ELSE IF( ncvt.LT.0 )
THEN
293 ELSE IF( nru.LT.0 )
THEN
295 ELSE IF( ncc.LT.0 )
THEN
297 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
298 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) )
THEN
300 ELSE IF( ldu.LT.max( 1, nru ) )
THEN
302 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
303 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) )
THEN
307 CALL
xerbla(
'ZBDSQR', -info )
317 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
321 IF( .NOT.rotate )
THEN
322 CALL
dlasq1( n, d, e, rwork, info )
326 IF( info .NE. 2 )
RETURN
338 unfl =
dlamch(
'Safe minimum' )
345 CALL
dlartg( d( i ), e( i ), cs, sn, r )
348 d( i+1 ) = cs*d( i+1 )
356 $ CALL
zlasr(
'R',
'V',
'F', nru, n, rwork( 1 ), rwork( n ),
359 $ CALL
zlasr(
'L',
'V',
'F', n, ncc, rwork( 1 ), rwork( n ),
367 tolmul = max( ten, min( hndrd, eps**meigth ) )
374 smax = max( smax, abs( d( i ) ) )
377 smax = max( smax, abs( e( i ) ) )
380 IF( tol.GE.zero )
THEN
384 sminoa = abs( d( 1 ) )
389 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
390 sminoa = min( sminoa, mu )
395 sminoa = sminoa / sqrt( dble( n ) )
396 thresh = max( tol*sminoa, maxitr*n*n*unfl )
401 thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
430 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
436 abss = abs( d( ll ) )
437 abse = abs( e( ll ) )
438 IF( tol.LT.zero .AND. abss.LE.thresh )
442 smin = min( smin, abss )
443 smax = max( smax, abss, abse )
468 CALL
dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
477 $ CALL
zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
480 $ CALL
zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
482 $ CALL
zdrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
491 IF( ll.GT.oldm .OR. m.LT.oldll )
THEN
492 IF( abs( d( ll ) ).GE.abs( d( m ) ) )
THEN
512 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
513 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) )
THEN
518 IF( tol.GE.zero )
THEN
525 DO 100 lll = ll, m - 1
526 IF( abs( e( lll ) ).LE.tol*mu )
THEN
530 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
531 sminl = min( sminl, mu )
540 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
541 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) )
THEN
546 IF( tol.GE.zero )
THEN
553 DO 110 lll = m - 1, ll, -1
554 IF( abs( e( lll ) ).LE.tol*mu )
THEN
558 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
559 sminl = min( sminl, mu )
569 IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
570 $ max( eps, hndrth*tol ) )
THEN
581 CALL
dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
584 CALL
dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
589 IF( sll.GT.zero )
THEN
590 IF( ( shift / sll )**2.LT.eps )
601 IF( shift.EQ.zero )
THEN
610 CALL
dlartg( d( i )*cs, e( i ), cs, sn, r )
613 CALL
dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
615 rwork( i-ll+1+nm1 ) = sn
616 rwork( i-ll+1+nm12 ) = oldcs
617 rwork( i-ll+1+nm13 ) = oldsn
626 $ CALL
zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
627 $ rwork( n ), vt( ll, 1 ), ldvt )
629 $ CALL
zlasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
630 $ rwork( nm13+1 ), u( 1, ll ), ldu )
632 $ CALL
zlasr(
'L',
'V',
'F', m-ll+1, ncc, rwork( nm12+1 ),
633 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
637 IF( abs( e( m-1 ) ).LE.thresh )
647 DO 130 i = m, ll + 1, -1
648 CALL
dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
651 CALL
dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
653 rwork( i-ll+nm1 ) = -sn
654 rwork( i-ll+nm12 ) = oldcs
655 rwork( i-ll+nm13 ) = -oldsn
664 $ CALL
zlasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
665 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
667 $ CALL
zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
668 $ rwork( n ), u( 1, ll ), ldu )
670 $ CALL
zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
671 $ rwork( n ), c( ll, 1 ), ldc )
675 IF( abs( e( ll ) ).LE.thresh )
687 f = ( abs( d( ll ) )-shift )*
688 $ ( sign( one, d( ll ) )+shift / d( ll ) )
691 CALL
dlartg( f, g, cosr, sinr, r )
694 f = cosr*d( i ) + sinr*e( i )
695 e( i ) = cosr*e( i ) - sinr*d( i )
697 d( i+1 ) = cosr*d( i+1 )
698 CALL
dlartg( f, g, cosl, sinl, r )
700 f = cosl*e( i ) + sinl*d( i+1 )
701 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
704 e( i+1 ) = cosl*e( i+1 )
706 rwork( i-ll+1 ) = cosr
707 rwork( i-ll+1+nm1 ) = sinr
708 rwork( i-ll+1+nm12 ) = cosl
709 rwork( i-ll+1+nm13 ) = sinl
716 $ CALL
zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
717 $ rwork( n ), vt( ll, 1 ), ldvt )
719 $ CALL
zlasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
720 $ rwork( nm13+1 ), u( 1, ll ), ldu )
722 $ CALL
zlasr(
'L',
'V',
'F', m-ll+1, ncc, rwork( nm12+1 ),
723 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
727 IF( abs( e( m-1 ) ).LE.thresh )
735 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
738 DO 150 i = m, ll + 1, -1
739 CALL
dlartg( f, g, cosr, sinr, r )
742 f = cosr*d( i ) + sinr*e( i-1 )
743 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
745 d( i-1 ) = cosr*d( i-1 )
746 CALL
dlartg( f, g, cosl, sinl, r )
748 f = cosl*e( i-1 ) + sinl*d( i-1 )
749 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
752 e( i-2 ) = cosl*e( i-2 )
755 rwork( i-ll+nm1 ) = -sinr
756 rwork( i-ll+nm12 ) = cosl
757 rwork( i-ll+nm13 ) = -sinl
763 IF( abs( e( ll ) ).LE.thresh )
769 $ CALL
zlasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
770 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
772 $ CALL
zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
773 $ rwork( n ), u( 1, ll ), ldu )
775 $ CALL
zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
776 $ rwork( n ), c( ll, 1 ), ldc )
788 IF( d( i ).LT.zero )
THEN
794 $ CALL
zdscal( ncvt, negone, vt( i, 1 ), ldvt )
807 DO 180
j = 2, n + 1 - i
808 IF( d(
j ).LE.smin )
THEN
813 IF( isub.NE.n+1-i )
THEN
817 d( isub ) = d( n+1-i )
820 $ CALL
zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
823 $ CALL
zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
825 $ CALL
zswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
LOGICAL function lsame(CA, CB)
LSAME
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 zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.