220 SUBROUTINE cbdsqr( 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 REAL D( * ), E( * ), RWORK( * )
233 COMPLEX C( LDC, * ), U( LDU, * ), VT( LDVT, * )
240 parameter( zero = 0.0e0 )
242 parameter( one = 1.0e0 )
244 parameter( negone = -1.0e0 )
246 parameter( hndrth = 0.01e0 )
248 parameter( ten = 10.0e0 )
250 parameter( hndrd = 100.0e0 )
252 parameter( meigth = -0.125e0 )
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 REAL 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
268 EXTERNAL lsame, slamch
275 INTRINSIC abs, max, min, real, 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(
'CBDSQR', -info )
313 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
317 IF( .NOT.rotate )
THEN
318 CALL slasq1( n, d, e, rwork, info )
322 IF( info .NE. 2 )
RETURN
333 eps = slamch(
'Epsilon' )
334 unfl = slamch(
'Safe minimum' )
341 CALL slartg( d( i ), e( i ), cs, sn, r )
344 d( i+1 ) = cs*d( i+1 )
352 $
CALL clasr(
'R',
'V',
'F', nru, n, rwork( 1 ), rwork( n ),
355 $
CALL clasr(
'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( real( 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 slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
473 $
CALL csrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
476 $
CALL csrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
478 $
CALL csrot( 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 slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
580 CALL slas2( 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 slartg( d( i )*cs, e( i ), cs, sn, r )
609 CALL slartg( 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 clasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
623 $ rwork( n ), vt( ll, 1 ), ldvt )
625 $
CALL clasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
626 $ rwork( nm13+1 ), u( 1, ll ), ldu )
628 $
CALL clasr(
'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 slartg( d( i )*cs, e( i-1 ), cs, sn, r )
647 CALL slartg( 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 clasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
661 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
663 $
CALL clasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
664 $ rwork( n ), u( 1, ll ), ldu )
666 $
CALL clasr(
'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 slartg( 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 slartg( 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 clasr(
'L',
'V',
'F', m-ll+1, ncvt, rwork( 1 ),
713 $ rwork( n ), vt( ll, 1 ), ldvt )
715 $
CALL clasr(
'R',
'V',
'F', nru, m-ll+1, rwork( nm12+1 ),
716 $ rwork( nm13+1 ), u( 1, ll ), ldu )
718 $
CALL clasr(
'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 slartg( 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 slartg( 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 clasr(
'L',
'V',
'B', m-ll+1, ncvt, rwork( nm12+1 ),
766 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
768 $
CALL clasr(
'R',
'V',
'B', nru, m-ll+1, rwork( 1 ),
769 $ rwork( n ), u( 1, ll ), ldu )
771 $
CALL clasr(
'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 csscal( 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 cswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
819 $
CALL cswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
821 $
CALL cswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slasq1(N, D, E, WORK, INFO)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
subroutine csscal(N, SA, CX, INCX)
CSSCAL
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR