223 SUBROUTINE cbdsqr( 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 REAL d( * ), e( * ), rwork( * )
237 COMPLEX c( ldc, * ), u( ldu, * ), vt( ldvt, * )
244 parameter( zero = 0.0e0 )
246 parameter( one = 1.0e0 )
248 parameter( negone = -1.0e0 )
250 parameter( hndrth = 0.01e0 )
252 parameter( ten = 10.0e0 )
254 parameter( hndrd = 100.0e0 )
256 parameter( meigth = -0.125e0 )
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 REAL 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, max, min,
REAL, 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(
'CBDSQR', -info )
317 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
321 IF( .NOT.rotate )
THEN
322 CALL
slasq1( n, d, e, rwork, info )
326 IF( info .NE. 2 ) return
338 unfl =
slamch(
'Safe minimum' )
345 CALL
slartg( d( i ), e( i ), cs, sn, r )
348 d( i+1 ) = cs*d( i+1 )
356 $ CALL
clasr(
'R',
'V',
'F', nru, n, rwork( 1 ), rwork( n ),
359 $ CALL
clasr(
'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(
REAL( 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
slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
477 $ CALL
csrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
480 $ CALL
csrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
482 $ CALL
csrot( 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
slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
584 CALL
slas2( 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
slartg( d( i )*cs, e( i ), cs, sn, r )
613 CALL
slartg( 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
clasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
627 $ rwork( n ), vt( ll, 1 ), ldvt )
629 $ CALL
clasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
630 $ rwork( nm13+1 ), u( 1, ll ), ldu )
632 $ CALL
clasr(
'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
slartg( d( i )*cs, e( i-1 ), cs, sn, r )
651 CALL
slartg( 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
clasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
665 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
667 $ CALL
clasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
668 $ rwork( n ), u( 1, ll ), ldu )
670 $ CALL
clasr(
'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
slartg( 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
slartg( 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
clasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
717 $ rwork( n ), vt( ll, 1 ), ldvt )
719 $ CALL
clasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
720 $ rwork( nm13+1 ), u( 1, ll ), ldu )
722 $ CALL
clasr(
'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
slartg( 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
slartg( 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
clasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
770 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
772 $ CALL
clasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
773 $ rwork( n ), u( 1, ll ), ldu )
775 $ CALL
clasr(
'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
csscal( 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
cswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
823 $ CALL
cswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
825 $ CALL
cswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )