231 SUBROUTINE zbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
232 $ LDU, C, LDC, RWORK, INFO )
240 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
243 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
244 COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
250 DOUBLE PRECISION ZERO
251 parameter( zero = 0.0d0 )
253 parameter( one = 1.0d0 )
254 DOUBLE PRECISION NEGONE
255 parameter( negone = -1.0d0 )
256 DOUBLE PRECISION HNDRTH
257 parameter( hndrth = 0.01d0 )
259 parameter( ten = 10.0d0 )
260 DOUBLE PRECISION HNDRD
261 parameter( hndrd = 100.0d0 )
262 DOUBLE PRECISION MEIGTH
263 parameter( meigth = -0.125d0 )
265 parameter( maxitr = 6 )
268 LOGICAL LOWER, ROTATE
269 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
270 $ maxitdivn, nm1, nm12, nm13, oldll, oldm
271 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
272 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
273 $ sinr, sll, smax, smin, sminoa,
274 $ sn, thresh, tol, tolmul, unfl
278 DOUBLE PRECISION DLAMCH
279 EXTERNAL lsame, dlamch
287 INTRINSIC abs, dble, max, min, sign, sqrt
294 lower = lsame( uplo,
'L' )
295 IF( .NOT.lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
297 ELSE IF( n.LT.0 )
THEN
299 ELSE IF( ncvt.LT.0 )
THEN
301 ELSE IF( nru.LT.0 )
THEN
303 ELSE IF( ncc.LT.0 )
THEN
305 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
306 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) )
THEN
308 ELSE IF( ldu.LT.max( 1, nru ) )
THEN
310 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
311 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) )
THEN
315 CALL xerbla(
'ZBDSQR', -info )
325 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
329 IF( .NOT.rotate )
THEN
330 CALL dlasq1( n, d, e, rwork, info )
334 IF( info .NE. 2 )
RETURN
345 eps = dlamch(
'Epsilon' )
346 unfl = dlamch(
'Safe minimum' )
353 CALL dlartg( d( i ), e( i ), cs, sn, r )
356 d( i+1 ) = cs*d( i+1 )
364 $
CALL zlasr(
'R',
'V',
'F', nru, n, rwork( 1 ),
368 $
CALL zlasr(
'L',
'V',
'F', n, ncc, rwork( 1 ),
377 tolmul = max( ten, min( hndrd, eps**meigth ) )
384 smax = max( smax, abs( d( i ) ) )
387 smax = max( smax, abs( e( i ) ) )
390 IF( tol.GE.zero )
THEN
394 sminoa = abs( d( 1 ) )
399 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
400 sminoa = min( sminoa, mu )
405 sminoa = sminoa / sqrt( dble( n ) )
406 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
411 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
438 iterdivn = iterdivn + 1
439 IF( iterdivn.GE.maxitdivn )
445 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
450 abss = abs( d( ll ) )
451 abse = abs( e( ll ) )
452 IF( tol.LT.zero .AND. abss.LE.thresh )
456 smax = max( smax, abss, abse )
481 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
490 $
CALL zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
493 $
CALL zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
496 $
CALL zdrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
505 IF( ll.GT.oldm .OR. m.LT.oldll )
THEN
506 IF( abs( d( ll ) ).GE.abs( d( m ) ) )
THEN
526 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
527 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) )
THEN
532 IF( tol.GE.zero )
THEN
539 DO 100 lll = ll, m - 1
540 IF( abs( e( lll ) ).LE.tol*mu )
THEN
544 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
545 smin = min( smin, mu )
554 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
555 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) )
THEN
560 IF( tol.GE.zero )
THEN
567 DO 110 lll = m - 1, ll, -1
568 IF( abs( e( lll ) ).LE.tol*mu )
THEN
572 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
573 smin = min( smin, mu )
583 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
584 $ max( eps, hndrth*tol ) )
THEN
595 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
598 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
603 IF( sll.GT.zero )
THEN
604 IF( ( shift / sll )**2.LT.eps )
615 IF( shift.EQ.zero )
THEN
624 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
627 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn,
630 rwork( i-ll+1+nm1 ) = sn
631 rwork( i-ll+1+nm12 ) = oldcs
632 rwork( i-ll+1+nm13 ) = oldsn
641 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
642 $ rwork( n ), vt( ll, 1 ), ldvt )
644 $
CALL zlasr(
'R',
'V',
'F', nru, m-ll+1,
646 $ rwork( nm13+1 ), u( 1, ll ), ldu )
648 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncc,
650 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
654 IF( abs( e( m-1 ) ).LE.thresh )
664 DO 130 i = m, ll + 1, -1
665 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
668 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn,
671 rwork( i-ll+nm1 ) = -sn
672 rwork( i-ll+nm12 ) = oldcs
673 rwork( i-ll+nm13 ) = -oldsn
682 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncvt,
684 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
686 $
CALL zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
687 $ rwork( n ), u( 1, ll ), ldu )
689 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
690 $ rwork( n ), c( ll, 1 ), ldc )
694 IF( abs( e( ll ) ).LE.thresh )
706 f = ( abs( d( ll ) )-shift )*
707 $ ( sign( one, d( ll ) )+shift / d( ll ) )
710 CALL dlartg( f, g, cosr, sinr, r )
713 f = cosr*d( i ) + sinr*e( i )
714 e( i ) = cosr*e( i ) - sinr*d( i )
716 d( i+1 ) = cosr*d( i+1 )
717 CALL dlartg( f, g, cosl, sinl, r )
719 f = cosl*e( i ) + sinl*d( i+1 )
720 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
723 e( i+1 ) = cosl*e( i+1 )
725 rwork( i-ll+1 ) = cosr
726 rwork( i-ll+1+nm1 ) = sinr
727 rwork( i-ll+1+nm12 ) = cosl
728 rwork( i-ll+1+nm13 ) = sinl
735 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
736 $ rwork( n ), vt( ll, 1 ), ldvt )
738 $
CALL zlasr(
'R',
'V',
'F', nru, m-ll+1,
740 $ rwork( nm13+1 ), u( 1, ll ), ldu )
742 $
CALL zlasr(
'L',
'V',
'F', m-ll+1, ncc,
744 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
748 IF( abs( e( m-1 ) ).LE.thresh )
756 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
759 DO 150 i = m, ll + 1, -1
760 CALL dlartg( f, g, cosr, sinr, r )
763 f = cosr*d( i ) + sinr*e( i-1 )
764 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
766 d( i-1 ) = cosr*d( i-1 )
767 CALL dlartg( f, g, cosl, sinl, r )
769 f = cosl*e( i-1 ) + sinl*d( i-1 )
770 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
773 e( i-2 ) = cosl*e( i-2 )
776 rwork( i-ll+nm1 ) = -sinr
777 rwork( i-ll+nm12 ) = cosl
778 rwork( i-ll+nm13 ) = -sinl
784 IF( abs( e( ll ) ).LE.thresh )
790 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncvt,
792 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
794 $
CALL zlasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
795 $ rwork( n ), u( 1, ll ), ldu )
797 $
CALL zlasr(
'L',
'V',
'B', m-ll+1, ncc, rwork( 1 ),
798 $ rwork( n ), c( ll, 1 ), ldc )
810 IF( d( i ).EQ.zero )
THEN
816 IF( d( i ).LT.zero )
THEN
822 $
CALL zdscal( ncvt, negone, vt( i, 1 ), ldvt )
835 DO 180 j = 2, n + 1 - i
836 IF( d( j ).LE.smin )
THEN
841 IF( isub.NE.n+1-i )
THEN
845 d( isub ) = d( n+1-i )
848 $
CALL zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
851 $
CALL zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
853 $
CALL zswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ),