Annotation of rpl/lapack/lapack/dlasv2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
! 2: *
! 3: * -- LAPACK auxiliary routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: DOUBLE PRECISION CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
! 10: * ..
! 11: *
! 12: * Purpose
! 13: * =======
! 14: *
! 15: * DLASV2 computes the singular value decomposition of a 2-by-2
! 16: * triangular matrix
! 17: * [ F G ]
! 18: * [ 0 H ].
! 19: * On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
! 20: * smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
! 21: * right singular vectors for abs(SSMAX), giving the decomposition
! 22: *
! 23: * [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
! 24: * [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
! 25: *
! 26: * Arguments
! 27: * =========
! 28: *
! 29: * F (input) DOUBLE PRECISION
! 30: * The (1,1) element of the 2-by-2 matrix.
! 31: *
! 32: * G (input) DOUBLE PRECISION
! 33: * The (1,2) element of the 2-by-2 matrix.
! 34: *
! 35: * H (input) DOUBLE PRECISION
! 36: * The (2,2) element of the 2-by-2 matrix.
! 37: *
! 38: * SSMIN (output) DOUBLE PRECISION
! 39: * abs(SSMIN) is the smaller singular value.
! 40: *
! 41: * SSMAX (output) DOUBLE PRECISION
! 42: * abs(SSMAX) is the larger singular value.
! 43: *
! 44: * SNL (output) DOUBLE PRECISION
! 45: * CSL (output) DOUBLE PRECISION
! 46: * The vector (CSL, SNL) is a unit left singular vector for the
! 47: * singular value abs(SSMAX).
! 48: *
! 49: * SNR (output) DOUBLE PRECISION
! 50: * CSR (output) DOUBLE PRECISION
! 51: * The vector (CSR, SNR) is a unit right singular vector for the
! 52: * singular value abs(SSMAX).
! 53: *
! 54: * Further Details
! 55: * ===============
! 56: *
! 57: * Any input parameter may be aliased with any output parameter.
! 58: *
! 59: * Barring over/underflow and assuming a guard digit in subtraction, all
! 60: * output quantities are correct to within a few units in the last
! 61: * place (ulps).
! 62: *
! 63: * In IEEE arithmetic, the code works correctly if one matrix element is
! 64: * infinite.
! 65: *
! 66: * Overflow will not occur unless the largest singular value itself
! 67: * overflows or is within a few ulps of overflow. (On machines with
! 68: * partial overflow, like the Cray, overflow may occur if the largest
! 69: * singular value is within a factor of 2 of overflow.)
! 70: *
! 71: * Underflow is harmless if underflow is gradual. Otherwise, results
! 72: * may correspond to a matrix modified by perturbations of size near
! 73: * the underflow threshold.
! 74: *
! 75: * =====================================================================
! 76: *
! 77: * .. Parameters ..
! 78: DOUBLE PRECISION ZERO
! 79: PARAMETER ( ZERO = 0.0D0 )
! 80: DOUBLE PRECISION HALF
! 81: PARAMETER ( HALF = 0.5D0 )
! 82: DOUBLE PRECISION ONE
! 83: PARAMETER ( ONE = 1.0D0 )
! 84: DOUBLE PRECISION TWO
! 85: PARAMETER ( TWO = 2.0D0 )
! 86: DOUBLE PRECISION FOUR
! 87: PARAMETER ( FOUR = 4.0D0 )
! 88: * ..
! 89: * .. Local Scalars ..
! 90: LOGICAL GASMAL, SWAP
! 91: INTEGER PMAX
! 92: DOUBLE PRECISION A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
! 93: $ MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
! 94: * ..
! 95: * .. Intrinsic Functions ..
! 96: INTRINSIC ABS, SIGN, SQRT
! 97: * ..
! 98: * .. External Functions ..
! 99: DOUBLE PRECISION DLAMCH
! 100: EXTERNAL DLAMCH
! 101: * ..
! 102: * .. Executable Statements ..
! 103: *
! 104: FT = F
! 105: FA = ABS( FT )
! 106: HT = H
! 107: HA = ABS( H )
! 108: *
! 109: * PMAX points to the maximum absolute element of matrix
! 110: * PMAX = 1 if F largest in absolute values
! 111: * PMAX = 2 if G largest in absolute values
! 112: * PMAX = 3 if H largest in absolute values
! 113: *
! 114: PMAX = 1
! 115: SWAP = ( HA.GT.FA )
! 116: IF( SWAP ) THEN
! 117: PMAX = 3
! 118: TEMP = FT
! 119: FT = HT
! 120: HT = TEMP
! 121: TEMP = FA
! 122: FA = HA
! 123: HA = TEMP
! 124: *
! 125: * Now FA .ge. HA
! 126: *
! 127: END IF
! 128: GT = G
! 129: GA = ABS( GT )
! 130: IF( GA.EQ.ZERO ) THEN
! 131: *
! 132: * Diagonal matrix
! 133: *
! 134: SSMIN = HA
! 135: SSMAX = FA
! 136: CLT = ONE
! 137: CRT = ONE
! 138: SLT = ZERO
! 139: SRT = ZERO
! 140: ELSE
! 141: GASMAL = .TRUE.
! 142: IF( GA.GT.FA ) THEN
! 143: PMAX = 2
! 144: IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
! 145: *
! 146: * Case of very large GA
! 147: *
! 148: GASMAL = .FALSE.
! 149: SSMAX = GA
! 150: IF( HA.GT.ONE ) THEN
! 151: SSMIN = FA / ( GA / HA )
! 152: ELSE
! 153: SSMIN = ( FA / GA )*HA
! 154: END IF
! 155: CLT = ONE
! 156: SLT = HT / GT
! 157: SRT = ONE
! 158: CRT = FT / GT
! 159: END IF
! 160: END IF
! 161: IF( GASMAL ) THEN
! 162: *
! 163: * Normal case
! 164: *
! 165: D = FA - HA
! 166: IF( D.EQ.FA ) THEN
! 167: *
! 168: * Copes with infinite F or H
! 169: *
! 170: L = ONE
! 171: ELSE
! 172: L = D / FA
! 173: END IF
! 174: *
! 175: * Note that 0 .le. L .le. 1
! 176: *
! 177: M = GT / FT
! 178: *
! 179: * Note that abs(M) .le. 1/macheps
! 180: *
! 181: T = TWO - L
! 182: *
! 183: * Note that T .ge. 1
! 184: *
! 185: MM = M*M
! 186: TT = T*T
! 187: S = SQRT( TT+MM )
! 188: *
! 189: * Note that 1 .le. S .le. 1 + 1/macheps
! 190: *
! 191: IF( L.EQ.ZERO ) THEN
! 192: R = ABS( M )
! 193: ELSE
! 194: R = SQRT( L*L+MM )
! 195: END IF
! 196: *
! 197: * Note that 0 .le. R .le. 1 + 1/macheps
! 198: *
! 199: A = HALF*( S+R )
! 200: *
! 201: * Note that 1 .le. A .le. 1 + abs(M)
! 202: *
! 203: SSMIN = HA / A
! 204: SSMAX = FA*A
! 205: IF( MM.EQ.ZERO ) THEN
! 206: *
! 207: * Note that M is very tiny
! 208: *
! 209: IF( L.EQ.ZERO ) THEN
! 210: T = SIGN( TWO, FT )*SIGN( ONE, GT )
! 211: ELSE
! 212: T = GT / SIGN( D, FT ) + M / T
! 213: END IF
! 214: ELSE
! 215: T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
! 216: END IF
! 217: L = SQRT( T*T+FOUR )
! 218: CRT = TWO / L
! 219: SRT = T / L
! 220: CLT = ( CRT+SRT*M ) / A
! 221: SLT = ( HT / FT )*SRT / A
! 222: END IF
! 223: END IF
! 224: IF( SWAP ) THEN
! 225: CSL = SRT
! 226: SNL = CRT
! 227: CSR = SLT
! 228: SNR = CLT
! 229: ELSE
! 230: CSL = CLT
! 231: SNL = SLT
! 232: CSR = CRT
! 233: SNR = SRT
! 234: END IF
! 235: *
! 236: * Correct signs of SSMAX and SSMIN
! 237: *
! 238: IF( PMAX.EQ.1 )
! 239: $ TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
! 240: IF( PMAX.EQ.2 )
! 241: $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
! 242: IF( PMAX.EQ.3 )
! 243: $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
! 244: SSMAX = SIGN( SSMAX, TSIGN )
! 245: SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
! 246: RETURN
! 247: *
! 248: * End of DLASV2
! 249: *
! 250: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasv2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack