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 )
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR