230 SUBROUTINE sbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
231 $ ldu, c, ldc, work, info )
240 INTEGER info, ldc, ldu, ldvt, n, ncc, ncvt, nru
243 REAL c( ldc, * ), d( * ), e( * ), u( ldu, * ),
244 $ vt( ldvt, * ), work( * )
251 parameter( zero = 0.0e0 )
253 parameter( one = 1.0e0 )
255 parameter( negone = -1.0e0 )
257 parameter( hndrth = 0.01e0 )
259 parameter( ten = 10.0e0 )
261 parameter( hndrd = 100.0e0 )
263 parameter( meigth = -0.125e0 )
265 parameter( maxitr = 6 )
268 LOGICAL lower, rotate
269 INTEGER i, idir, isub, iter,
j, ll, lll, m, maxit, nm1,
270 $ nm12, nm13, oldll, oldm
271 REAL abse, abss, cosl, cosr, cs, eps, f, g, h, mu,
272 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
273 $ sinr, sll, smax, smin, sminl, sminoa,
274 $ sn, thresh, tol, tolmul, unfl
286 INTRINSIC abs, max, min,
REAL, sign, sqrt
293 lower =
lsame( uplo,
'L' )
294 IF( .NOT.
lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
296 ELSE IF( n.LT.0 )
THEN
298 ELSE IF( ncvt.LT.0 )
THEN
300 ELSE IF( nru.LT.0 )
THEN
302 ELSE IF( ncc.LT.0 )
THEN
304 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
305 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) )
THEN
307 ELSE IF( ldu.LT.max( 1, nru ) )
THEN
309 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
310 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) )
THEN
314 CALL
xerbla(
'SBDSQR', -info )
324 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
328 IF( .NOT.rotate )
THEN
329 CALL
slasq1( n, d, e, work, info )
333 IF( info .NE. 2 )
RETURN
345 unfl =
slamch(
'Safe minimum' )
352 CALL
slartg( d( i ), e( i ), cs, sn, r )
355 d( i+1 ) = cs*d( i+1 )
363 $ CALL
slasr(
'R',
'V',
'F', nru, n, work( 1 ), work( n ), u,
366 $ CALL
slasr(
'L',
'V',
'F', n, ncc, work( 1 ), work( n ), c,
374 tolmul = max( ten, min( hndrd, eps**meigth ) )
381 smax = max( smax, abs( d( i ) ) )
384 smax = max( smax, abs( e( i ) ) )
387 IF( tol.GE.zero )
THEN
391 sminoa = abs( d( 1 ) )
396 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
397 sminoa = min( sminoa, mu )
402 sminoa = sminoa / sqrt(
REAL( N ) )
403 thresh = max( tol*sminoa, maxitr*n*n*unfl )
408 thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
437 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
443 abss = abs( d( ll ) )
444 abse = abs( e( ll ) )
445 IF( tol.LT.zero .AND. abss.LE.thresh )
449 smin = min( smin, abss )
450 smax = max( smax, abss, abse )
475 CALL
slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
484 $ CALL
srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
487 $ CALL
srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
489 $ CALL
srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
498 IF( ll.GT.oldm .OR. m.LT.oldll )
THEN
499 IF( abs( d( ll ) ).GE.abs( d( m ) ) )
THEN
519 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
520 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) )
THEN
525 IF( tol.GE.zero )
THEN
532 DO 100 lll = ll, m - 1
533 IF( abs( e( lll ) ).LE.tol*mu )
THEN
537 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
538 sminl = min( sminl, mu )
547 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
548 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) )
THEN
553 IF( tol.GE.zero )
THEN
560 DO 110 lll = m - 1, ll, -1
561 IF( abs( e( lll ) ).LE.tol*mu )
THEN
565 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
566 sminl = min( sminl, mu )
576 IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
577 $ max( eps, hndrth*tol ) )
THEN
588 CALL
slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
591 CALL
slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
596 IF( sll.GT.zero )
THEN
597 IF( ( shift / sll )**2.LT.eps )
608 IF( shift.EQ.zero )
THEN
617 CALL
slartg( d( i )*cs, e( i ), cs, sn, r )
620 CALL
slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
622 work( i-ll+1+nm1 ) = sn
623 work( i-ll+1+nm12 ) = oldcs
624 work( i-ll+1+nm13 ) = oldsn
633 $ CALL
slasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
634 $ work( n ), vt( ll, 1 ), ldvt )
636 $ CALL
slasr(
'R',
'V',
'F', nru, m-ll+1, work( nm12+1 ),
637 $ work( nm13+1 ), u( 1, ll ), ldu )
639 $ CALL
slasr(
'L',
'V',
'F', m-ll+1, ncc, work( nm12+1 ),
640 $ work( nm13+1 ), c( ll, 1 ), ldc )
644 IF( abs( e( m-1 ) ).LE.thresh )
654 DO 130 i = m, ll + 1, -1
655 CALL
slartg( d( i )*cs, e( i-1 ), cs, sn, r )
658 CALL
slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
660 work( i-ll+nm1 ) = -sn
661 work( i-ll+nm12 ) = oldcs
662 work( i-ll+nm13 ) = -oldsn
671 $ CALL
slasr(
'L',
'V',
'B', m-ll+1, ncvt, work( nm12+1 ),
672 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
674 $ CALL
slasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
675 $ work( n ), u( 1, ll ), ldu )
677 $ CALL
slasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
678 $ work( n ), c( ll, 1 ), ldc )
682 IF( abs( e( ll ) ).LE.thresh )
694 f = ( abs( d( ll ) )-shift )*
695 $ ( sign( one, d( ll ) )+shift / d( ll ) )
698 CALL
slartg( f, g, cosr, sinr, r )
701 f = cosr*d( i ) + sinr*e( i )
702 e( i ) = cosr*e( i ) - sinr*d( i )
704 d( i+1 ) = cosr*d( i+1 )
705 CALL
slartg( f, g, cosl, sinl, r )
707 f = cosl*e( i ) + sinl*d( i+1 )
708 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
711 e( i+1 ) = cosl*e( i+1 )
713 work( i-ll+1 ) = cosr
714 work( i-ll+1+nm1 ) = sinr
715 work( i-ll+1+nm12 ) = cosl
716 work( i-ll+1+nm13 ) = sinl
723 $ CALL
slasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
724 $ work( n ), vt( ll, 1 ), ldvt )
726 $ CALL
slasr(
'R',
'V',
'F', nru, m-ll+1, work( nm12+1 ),
727 $ work( nm13+1 ), u( 1, ll ), ldu )
729 $ CALL
slasr(
'L',
'V',
'F', m-ll+1, ncc, work( nm12+1 ),
730 $ work( nm13+1 ), c( ll, 1 ), ldc )
734 IF( abs( e( m-1 ) ).LE.thresh )
742 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
745 DO 150 i = m, ll + 1, -1
746 CALL
slartg( f, g, cosr, sinr, r )
749 f = cosr*d( i ) + sinr*e( i-1 )
750 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
752 d( i-1 ) = cosr*d( i-1 )
753 CALL
slartg( f, g, cosl, sinl, r )
755 f = cosl*e( i-1 ) + sinl*d( i-1 )
756 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
759 e( i-2 ) = cosl*e( i-2 )
762 work( i-ll+nm1 ) = -sinr
763 work( i-ll+nm12 ) = cosl
764 work( i-ll+nm13 ) = -sinl
770 IF( abs( e( ll ) ).LE.thresh )
776 $ CALL
slasr(
'L',
'V',
'B', m-ll+1, ncvt, work( nm12+1 ),
777 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
779 $ CALL
slasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
780 $ work( n ), u( 1, ll ), ldu )
782 $ CALL
slasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
783 $ work( n ), c( ll, 1 ), ldc )
795 IF( d( i ).LT.zero )
THEN
801 $ CALL
sscal( ncvt, negone, vt( i, 1 ), ldvt )
814 DO 180
j = 2, n + 1 - i
815 IF( d(
j ).LE.smin )
THEN
820 IF( isub.NE.n+1-i )
THEN
824 d( isub ) = d( n+1-i )
827 $ CALL
sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
830 $ CALL
sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
832 $ CALL
sswap( 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
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
REAL function slamch(CMACH)
SLAMCH
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
subroutine xerbla(SRNAME, INFO)
XERBLA
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 slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT