238 SUBROUTINE sbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
239 $ LDU, C, LDC, WORK, INFO )
247 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
250 REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
251 $ vt( ldvt, * ), work( * )
258 parameter( zero = 0.0e0 )
260 parameter( one = 1.0e0 )
262 parameter( negone = -1.0e0 )
264 parameter( hndrth = 0.01e0 )
266 parameter( ten = 10.0e0 )
268 parameter( hndrd = 100.0e0 )
270 parameter( meigth = -0.125e0 )
272 parameter( maxitr = 6 )
275 LOGICAL LOWER, ROTATE
276 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
277 $ maxitdivn, nm1, nm12, nm13, oldll, oldm
278 REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
279 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
280 $ sinr, sll, smax, smin, sminoa,
281 $ sn, thresh, tol, tolmul, unfl
286 EXTERNAL lsame, slamch
294 INTRINSIC abs, max, min, real, sign, sqrt
301 lower = lsame( uplo,
'L' )
302 IF( .NOT.lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
304 ELSE IF( n.LT.0 )
THEN
306 ELSE IF( ncvt.LT.0 )
THEN
308 ELSE IF( nru.LT.0 )
THEN
310 ELSE IF( ncc.LT.0 )
THEN
312 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
313 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) )
THEN
315 ELSE IF( ldu.LT.max( 1, nru ) )
THEN
317 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
318 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) )
THEN
322 CALL xerbla(
'SBDSQR', -info )
332 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
336 IF( .NOT.rotate )
THEN
337 CALL slasq1( n, d, e, work, info )
341 IF( info .NE. 2 )
RETURN
352 eps = slamch(
'Epsilon' )
353 unfl = slamch(
'Safe minimum' )
360 CALL slartg( d( i ), e( i ), cs, sn, r )
363 d( i+1 ) = cs*d( i+1 )
371 $
CALL slasr(
'R',
'V',
'F', nru, n, work( 1 ), work( n ),
375 $
CALL slasr(
'L',
'V',
'F', n, ncc, work( 1 ), work( n ),
384 tolmul = max( ten, min( hndrd, eps**meigth ) )
391 smax = max( smax, abs( d( i ) ) )
394 smax = max( smax, abs( e( i ) ) )
397 IF( tol.GE.zero )
THEN
401 sminoa = abs( d( 1 ) )
406 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
407 sminoa = min( sminoa, mu )
412 sminoa = sminoa / sqrt( real( n ) )
413 thresh = max( tol*sminoa, maxitr*(real( n )*(real( n )*unfl)) )
418 thresh = max( abs( tol )*smax, maxitr*(real( n )*
419 $ (real( n )*unfl)) )
447 iterdivn = iterdivn + 1
448 IF( iterdivn.GE.maxitdivn )
454 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
459 abss = abs( d( ll ) )
460 abse = abs( e( ll ) )
461 IF( tol.LT.zero .AND. abss.LE.thresh )
465 smax = max( smax, abss, abse )
490 CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
499 $
CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
503 $
CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
506 $
CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
515 IF( ll.GT.oldm .OR. m.LT.oldll )
THEN
516 IF( abs( d( ll ) ).GE.abs( d( m ) ) )
THEN
536 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
537 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) )
THEN
542 IF( tol.GE.zero )
THEN
549 DO 100 lll = ll, m - 1
550 IF( abs( e( lll ) ).LE.tol*mu )
THEN
554 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
555 smin = min( smin, mu )
564 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
565 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) )
THEN
570 IF( tol.GE.zero )
THEN
577 DO 110 lll = m - 1, ll, -1
578 IF( abs( e( lll ) ).LE.tol*mu )
THEN
582 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
583 smin = min( smin, mu )
593 IF( tol.GE.zero .AND. real( n )*tol*( smin / smax ).LE.
594 $ max( eps, hndrth*tol ) )
THEN
605 CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
608 CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
613 IF( sll.GT.zero )
THEN
614 IF( ( shift / sll )**2.LT.eps )
625 IF( shift.EQ.zero )
THEN
634 CALL slartg( d( i )*cs, e( i ), cs, sn, r )
637 CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn,
640 work( i-ll+1+nm1 ) = sn
641 work( i-ll+1+nm12 ) = oldcs
642 work( i-ll+1+nm13 ) = oldsn
651 $
CALL slasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
652 $ work( n ), vt( ll, 1 ), ldvt )
654 $
CALL slasr(
'R',
'V',
'F', nru, m-ll+1,
656 $ work( nm13+1 ), u( 1, ll ), ldu )
658 $
CALL slasr(
'L',
'V',
'F', m-ll+1, ncc,
660 $ work( nm13+1 ), c( ll, 1 ), ldc )
664 IF( abs( e( m-1 ) ).LE.thresh )
674 DO 130 i = m, ll + 1, -1
675 CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
678 CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn,
681 work( i-ll+nm1 ) = -sn
682 work( i-ll+nm12 ) = oldcs
683 work( i-ll+nm13 ) = -oldsn
692 $
CALL slasr(
'L',
'V',
'B', m-ll+1, ncvt,
694 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
696 $
CALL slasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
697 $ work( n ), u( 1, ll ), ldu )
699 $
CALL slasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
700 $ work( n ), c( ll, 1 ), ldc )
704 IF( abs( e( ll ) ).LE.thresh )
716 f = ( abs( d( ll ) )-shift )*
717 $ ( sign( one, d( ll ) )+shift / d( ll ) )
720 CALL slartg( f, g, cosr, sinr, r )
723 f = cosr*d( i ) + sinr*e( i )
724 e( i ) = cosr*e( i ) - sinr*d( i )
726 d( i+1 ) = cosr*d( i+1 )
727 CALL slartg( f, g, cosl, sinl, r )
729 f = cosl*e( i ) + sinl*d( i+1 )
730 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
733 e( i+1 ) = cosl*e( i+1 )
735 work( i-ll+1 ) = cosr
736 work( i-ll+1+nm1 ) = sinr
737 work( i-ll+1+nm12 ) = cosl
738 work( i-ll+1+nm13 ) = sinl
745 $
CALL slasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
746 $ work( n ), vt( ll, 1 ), ldvt )
748 $
CALL slasr(
'R',
'V',
'F', nru, m-ll+1,
750 $ work( nm13+1 ), u( 1, ll ), ldu )
752 $
CALL slasr(
'L',
'V',
'F', m-ll+1, ncc,
754 $ work( nm13+1 ), c( ll, 1 ), ldc )
758 IF( abs( e( m-1 ) ).LE.thresh )
766 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
769 DO 150 i = m, ll + 1, -1
770 CALL slartg( f, g, cosr, sinr, r )
773 f = cosr*d( i ) + sinr*e( i-1 )
774 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
776 d( i-1 ) = cosr*d( i-1 )
777 CALL slartg( f, g, cosl, sinl, r )
779 f = cosl*e( i-1 ) + sinl*d( i-1 )
780 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
783 e( i-2 ) = cosl*e( i-2 )
786 work( i-ll+nm1 ) = -sinr
787 work( i-ll+nm12 ) = cosl
788 work( i-ll+nm13 ) = -sinl
794 IF( abs( e( ll ) ).LE.thresh )
800 $
CALL slasr(
'L',
'V',
'B', m-ll+1, ncvt,
802 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
804 $
CALL slasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
805 $ work( n ), u( 1, ll ), ldu )
807 $
CALL slasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
808 $ work( n ), c( ll, 1 ), ldc )
820 IF( d( i ).EQ.zero )
THEN
826 IF( d( i ).LT.zero )
THEN
833 $
CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
846 DO 180 j = 2, n + 1 - i
847 IF( d( j ).LE.smin )
THEN
852 IF( isub.NE.n+1-i )
THEN
856 d( isub ) = d( n+1-i )
859 $
CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
862 $
CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
864 $
CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ),