001: SUBROUTINE SLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL ) 002: * 003: * -- LAPACK auxiliary routine (version 3.2) -- 004: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 005: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 006: * November 2006 007: * 008: * .. Scalar Arguments .. 009: REAL CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN 010: * .. 011: * 012: * Purpose 013: * ======= 014: * 015: * SLASV2 computes the singular value decomposition of a 2-by-2 016: * triangular matrix 017: * [ F G ] 018: * [ 0 H ]. 019: * On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the 020: * smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and 021: * right singular vectors for abs(SSMAX), giving the decomposition 022: * 023: * [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ] 024: * [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ]. 025: * 026: * Arguments 027: * ========= 028: * 029: * F (input) REAL 030: * The (1,1) element of the 2-by-2 matrix. 031: * 032: * G (input) REAL 033: * The (1,2) element of the 2-by-2 matrix. 034: * 035: * H (input) REAL 036: * The (2,2) element of the 2-by-2 matrix. 037: * 038: * SSMIN (output) REAL 039: * abs(SSMIN) is the smaller singular value. 040: * 041: * SSMAX (output) REAL 042: * abs(SSMAX) is the larger singular value. 043: * 044: * SNL (output) REAL 045: * CSL (output) REAL 046: * The vector (CSL, SNL) is a unit left singular vector for the 047: * singular value abs(SSMAX). 048: * 049: * SNR (output) REAL 050: * CSR (output) REAL 051: * The vector (CSR, SNR) is a unit right singular vector for the 052: * singular value abs(SSMAX). 053: * 054: * Further Details 055: * =============== 056: * 057: * Any input parameter may be aliased with any output parameter. 058: * 059: * Barring over/underflow and assuming a guard digit in subtraction, all 060: * output quantities are correct to within a few units in the last 061: * place (ulps). 062: * 063: * In IEEE arithmetic, the code works correctly if one matrix element is 064: * infinite. 065: * 066: * Overflow will not occur unless the largest singular value itself 067: * overflows or is within a few ulps of overflow. (On machines with 068: * partial overflow, like the Cray, overflow may occur if the largest 069: * singular value is within a factor of 2 of overflow.) 070: * 071: * Underflow is harmless if underflow is gradual. Otherwise, results 072: * may correspond to a matrix modified by perturbations of size near 073: * the underflow threshold. 074: * 075: * ===================================================================== 076: * 077: * .. Parameters .. 078: REAL ZERO 079: PARAMETER ( ZERO = 0.0E0 ) 080: REAL HALF 081: PARAMETER ( HALF = 0.5E0 ) 082: REAL ONE 083: PARAMETER ( ONE = 1.0E0 ) 084: REAL TWO 085: PARAMETER ( TWO = 2.0E0 ) 086: REAL FOUR 087: PARAMETER ( FOUR = 4.0E0 ) 088: * .. 089: * .. Local Scalars .. 090: LOGICAL GASMAL, SWAP 091: INTEGER PMAX 092: REAL A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M, 093: $ MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT 094: * .. 095: * .. Intrinsic Functions .. 096: INTRINSIC ABS, SIGN, SQRT 097: * .. 098: * .. External Functions .. 099: REAL SLAMCH 100: EXTERNAL SLAMCH 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.SLAMCH( '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 SLASV2 249: * 250: END 251: