1 SUBROUTINE SLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL ) 2* 3* -- LAPACK auxiliary routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* October 31, 1992 7* 8* .. Scalar Arguments .. 9 REAL CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN 10* .. 11* 12* Purpose 13* ======= 14* 15* SLASV2 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) REAL 30* The (1,1) element of the 2-by-2 matrix. 31* 32* G (input) REAL 33* The (1,2) element of the 2-by-2 matrix. 34* 35* H (input) REAL 36* The (2,2) element of the 2-by-2 matrix. 37* 38* SSMIN (output) REAL 39* abs(SSMIN) is the smaller singular value. 40* 41* SSMAX (output) REAL 42* abs(SSMAX) is the larger singular value. 43* 44* SNL (output) REAL 45* CSL (output) REAL 46* The vector (CSL, SNL) is a unit left singular vector for the 47* singular value abs(SSMAX). 48* 49* SNR (output) REAL 50* CSR (output) REAL 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 REAL ZERO 79 PARAMETER ( ZERO = 0.0E0 ) 80 REAL HALF 81 PARAMETER ( HALF = 0.5E0 ) 82 REAL ONE 83 PARAMETER ( ONE = 1.0E0 ) 84 REAL TWO 85 PARAMETER ( TWO = 2.0E0 ) 86 REAL FOUR 87 PARAMETER ( FOUR = 4.0E0 ) 88* .. 89* .. Local Scalars .. 90 LOGICAL GASMAL, SWAP 91 INTEGER PMAX 92 REAL 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 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