SUBROUTINE SLASQ2( M, Q, E, QQ, EE, EPS, TOL2, SMALL2, SUP, KEND, $ INFO ) * * -- LAPACK routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. INTEGER INFO, KEND, M REAL EPS, SMALL2, SUP, TOL2 * .. * .. Array Arguments .. REAL E( * ), EE( * ), Q( * ), QQ( * ) * .. * * Purpose * ======= * * SLASQ2 computes the singular values of a real N-by-N unreduced * bidiagonal matrix with squared diagonal elements in Q and * squared off-diagonal elements in E. The singular values are * computed to relative accuracy TOL, barring over/underflow or * denormalization. * * Arguments * ========= * * M (input) INTEGER * The number of rows and columns in the matrix. M >= 0. * * Q (output) REAL array, dimension (M) * On normal exit, contains the squared singular values. * * E (workspace) REAL array, dimension (M) * * QQ (input/output) REAL array, dimension (M) * On entry, QQ contains the squared diagonal elements of the * bidiagonal matrix whose SVD is desired. * On exit, QQ is overwritten. * * EE (input/output) REAL array, dimension (M) * On entry, EE(1:N-1) contains the squared off-diagonal * elements of the bidiagonal matrix whose SVD is desired. * On exit, EE is overwritten. * * EPS (input) REAL * Machine epsilon. * * TOL2 (input) REAL * Desired relative accuracy of computed eigenvalues * as defined in SLASQ1. * * SMALL2 (input) REAL * A threshold value as defined in SLASQ1. * * SUP (input/output) REAL * Upper bound for the smallest eigenvalue. * * KEND (input/output) INTEGER * Index where minimum d occurs. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, the algorithm did not converge; i * specifies how many superdiagonals did not converge. * * ===================================================================== * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) REAL FOUR, HALF PARAMETER ( FOUR = 4.0E+0, HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER ICONV, IPHASE, ISP, N, OFF, OFF1 REAL QEMAX, SIGMA, XINF, XX, YY * .. * .. External Subroutines .. EXTERNAL SLASQ3 * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN, NINT, SQRT * .. * .. Executable Statements .. N = M * * Set the default maximum number of iterations * OFF = 0 OFF1 = OFF + 1 SIGMA = ZERO XINF = ZERO ICONV = 0 IPHASE = 2 * * Try deflation at the bottom * * 1x1 deflation * 10 CONTINUE IF( N.LE.2 ) $ GO TO 20 IF( EE( N-1 ).LE.MAX( QQ( N ), XINF, SMALL2 )*TOL2 ) THEN Q( N ) = QQ( N ) N = N - 1 IF( KEND.GT.N ) $ KEND = N SUP = MIN( QQ( N ), QQ( N-1 ) ) GO TO 10 END IF * * 2x2 deflation * IF( EE( N-2 ).LE.MAX( XINF, SMALL2, $ ( QQ( N ) / ( QQ( N )+EE( N-1 )+QQ( N-1 ) ) )*QQ( N-1 ) )* $ TOL2 ) THEN QEMAX = MAX( QQ( N ), QQ( N-1 ), EE( N-1 ) ) IF( QEMAX.NE.ZERO ) THEN IF( QEMAX.EQ.QQ( N-1 ) ) THEN XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX* $ SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) / $ QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) ) ELSE IF( QEMAX.EQ.QQ( N ) ) THEN XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX* $ SQRT( ( ( QQ( N-1 )-QQ( N )+EE( N-1 ) ) / $ QEMAX )**2+FOUR*EE( N-1 ) / QEMAX ) ) ELSE XX = HALF*( QQ( N )+QQ( N-1 )+EE( N-1 )+QEMAX* $ SQRT( ( ( QQ( N )-QQ( N-1 )+EE( N-1 ) ) / $ QEMAX )**2+FOUR*QQ( N-1 ) / QEMAX ) ) END IF YY = ( MAX( QQ( N ), QQ( N-1 ) ) / XX )* $ MIN( QQ( N ), QQ( N-1 ) ) ELSE XX = ZERO YY = ZERO END IF Q( N-1 ) = XX Q( N ) = YY N = N - 2 IF( KEND.GT.N ) $ KEND = N SUP = QQ( N ) GO TO 10 END IF * 20 CONTINUE IF( N.EQ.0 ) THEN * * The lower branch is finished * IF( OFF.EQ.0 ) THEN * * No upper branch; return to SLASQ1 * RETURN ELSE * * Going back to upper branch * XINF = ZERO IF( EE( OFF ).GT.ZERO ) THEN ISP = NINT( EE( OFF ) ) IPHASE = 1 ELSE ISP = -NINT( EE( OFF ) ) IPHASE = 2 END IF SIGMA = E( OFF ) N = OFF - ISP + 1 OFF1 = ISP OFF = OFF1 - 1 IF( N.LE.2 ) $ GO TO 20 IF( IPHASE.EQ.1 ) THEN SUP = MIN( Q( N+OFF ), Q( N-1+OFF ), Q( N-2+OFF ) ) ELSE SUP = MIN( QQ( N+OFF ), QQ( N-1+OFF ), QQ( N-2+OFF ) ) END IF KEND = 0 ICONV = -3 END IF ELSE IF( N.EQ.1 ) THEN * * 1x1 Solver * IF( IPHASE.EQ.1 ) THEN Q( OFF1 ) = Q( OFF1 ) + SIGMA ELSE Q( OFF1 ) = QQ( OFF1 ) + SIGMA END IF N = 0 GO TO 20 * * 2x2 Solver * ELSE IF( N.EQ.2 ) THEN IF( IPHASE.EQ.2 ) THEN QEMAX = MAX( QQ( N+OFF ), QQ( N-1+OFF ), EE( N-1+OFF ) ) IF( QEMAX.NE.ZERO ) THEN IF( QEMAX.EQ.QQ( N-1+OFF ) ) THEN XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+ $ QEMAX*SQRT( ( ( QQ( N+OFF )-QQ( N-1+OFF )+EE( N- $ 1+OFF ) ) / QEMAX )**2+FOUR*EE( OFF+N-1 ) / $ QEMAX ) ) ELSE IF( QEMAX.EQ.QQ( N+OFF ) ) THEN XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+ $ QEMAX*SQRT( ( ( QQ( N-1+OFF )-QQ( N+OFF )+EE( N- $ 1+OFF ) ) / QEMAX )**2+FOUR*EE( N-1+OFF ) / $ QEMAX ) ) ELSE XX = HALF*( QQ( N+OFF )+QQ( N-1+OFF )+EE( N-1+OFF )+ $ QEMAX*SQRT( ( ( QQ( N+OFF )-QQ( N-1+OFF )+EE( N- $ 1+OFF ) ) / QEMAX )**2+FOUR*QQ( N-1+OFF ) / $ QEMAX ) ) END IF YY = ( MAX( QQ( N+OFF ), QQ( N-1+OFF ) ) / XX )* $ MIN( QQ( N+OFF ), QQ( N-1+OFF ) ) ELSE XX = ZERO YY = ZERO END IF ELSE QEMAX = MAX( Q( N+OFF ), Q( N-1+OFF ), E( N-1+OFF ) ) IF( QEMAX.NE.ZERO ) THEN IF( QEMAX.EQ.Q( N-1+OFF ) ) THEN XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+ $ QEMAX*SQRT( ( ( Q( N+OFF )-Q( N-1+OFF )+E( N-1+ $ OFF ) ) / QEMAX )**2+FOUR*E( N-1+OFF ) / $ QEMAX ) ) ELSE IF( QEMAX.EQ.Q( N+OFF ) ) THEN XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+ $ QEMAX*SQRT( ( ( Q( N-1+OFF )-Q( N+OFF )+E( N-1+ $ OFF ) ) / QEMAX )**2+FOUR*E( N-1+OFF ) / $ QEMAX ) ) ELSE XX = HALF*( Q( N+OFF )+Q( N-1+OFF )+E( N-1+OFF )+ $ QEMAX*SQRT( ( ( Q( N+OFF )-Q( N-1+OFF )+E( N-1+ $ OFF ) ) / QEMAX )**2+FOUR*Q( N-1+OFF ) / $ QEMAX ) ) END IF YY = ( MAX( Q( N+OFF ), Q( N-1+OFF ) ) / XX )* $ MIN( Q( N+OFF ), Q( N-1+OFF ) ) ELSE XX = ZERO YY = ZERO END IF END IF Q( N-1+OFF ) = SIGMA + XX Q( N+OFF ) = YY + SIGMA N = 0 GO TO 20 END IF CALL SLASQ3( N, Q( OFF1 ), E( OFF1 ), QQ( OFF1 ), EE( OFF1 ), SUP, $ SIGMA, KEND, OFF, IPHASE, ICONV, EPS, TOL2, SMALL2 ) IF( SUP.LT.ZERO ) THEN INFO = N + OFF RETURN END IF OFF1 = OFF + 1 GO TO 20 * * End of SLASQ2 * END