 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 DOUBLE PRECISION CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * DLASV2 computes the singular value decomposition of a 2-by-2 00016 * triangular matrix 00017 * [ F G ] 00018 * [ 0 H ]. 00019 * On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the 00020 * smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and 00021 * right singular vectors for abs(SSMAX), giving the decomposition 00022 * 00023 * [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ] 00024 * [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ]. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * F (input) DOUBLE PRECISION 00030 * The (1,1) element of the 2-by-2 matrix. 00031 * 00032 * G (input) DOUBLE PRECISION 00033 * The (1,2) element of the 2-by-2 matrix. 00034 * 00035 * H (input) DOUBLE PRECISION 00036 * The (2,2) element of the 2-by-2 matrix. 00037 * 00038 * SSMIN (output) DOUBLE PRECISION 00039 * abs(SSMIN) is the smaller singular value. 00040 * 00041 * SSMAX (output) DOUBLE PRECISION 00042 * abs(SSMAX) is the larger singular value. 00043 * 00044 * SNL (output) DOUBLE PRECISION 00045 * CSL (output) DOUBLE PRECISION 00046 * The vector (CSL, SNL) is a unit left singular vector for the 00047 * singular value abs(SSMAX). 00048 * 00049 * SNR (output) DOUBLE PRECISION 00050 * CSR (output) DOUBLE PRECISION 00051 * The vector (CSR, SNR) is a unit right singular vector for the 00052 * singular value abs(SSMAX). 00053 * 00054 * Further Details 00055 * =============== 00056 * 00057 * Any input parameter may be aliased with any output parameter. 00058 * 00059 * Barring over/underflow and assuming a guard digit in subtraction, all 00060 * output quantities are correct to within a few units in the last 00061 * place (ulps). 00062 * 00063 * In IEEE arithmetic, the code works correctly if one matrix element is 00064 * infinite. 00065 * 00066 * Overflow will not occur unless the largest singular value itself 00067 * overflows or is within a few ulps of overflow. (On machines with 00068 * partial overflow, like the Cray, overflow may occur if the largest 00069 * singular value is within a factor of 2 of overflow.) 00070 * 00071 * Underflow is harmless if underflow is gradual. Otherwise, results 00072 * may correspond to a matrix modified by perturbations of size near 00073 * the underflow threshold. 00074 * 00075 * ===================================================================== 00076 * 00077 * .. Parameters .. 00078 DOUBLE PRECISION ZERO 00079 PARAMETER ( ZERO = 0.0D0 ) 00080 DOUBLE PRECISION HALF 00081 PARAMETER ( HALF = 0.5D0 ) 00082 DOUBLE PRECISION ONE 00083 PARAMETER ( ONE = 1.0D0 ) 00084 DOUBLE PRECISION TWO 00085 PARAMETER ( TWO = 2.0D0 ) 00086 DOUBLE PRECISION FOUR 00087 PARAMETER ( FOUR = 4.0D0 ) 00088 * .. 00089 * .. Local Scalars .. 00090 LOGICAL GASMAL, SWAP 00091 INTEGER PMAX 00092 DOUBLE PRECISION A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M, 00093 $ MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT 00094 * .. 00095 * .. Intrinsic Functions .. 00096 INTRINSIC ABS, SIGN, SQRT 00097 * .. 00098 * .. External Functions .. 00099 DOUBLE PRECISION DLAMCH 00100 EXTERNAL DLAMCH 00101 * .. 00102 * .. Executable Statements .. 00103 * 00104 FT = F 00105 FA = ABS( FT ) 00106 HT = H 00107 HA = ABS( H ) 00108 * 00109 * PMAX points to the maximum absolute element of matrix 00110 * PMAX = 1 if F largest in absolute values 00111 * PMAX = 2 if G largest in absolute values 00112 * PMAX = 3 if H largest in absolute values 00113 * 00114 PMAX = 1 00115 SWAP = ( HA.GT.FA ) 00116 IF( SWAP ) THEN 00117 PMAX = 3 00118 TEMP = FT 00119 FT = HT 00120 HT = TEMP 00121 TEMP = FA 00122 FA = HA 00123 HA = TEMP 00124 * 00125 * Now FA .ge. HA 00126 * 00127 END IF 00128 GT = G 00129 GA = ABS( GT ) 00130 IF( GA.EQ.ZERO ) THEN 00131 * 00132 * Diagonal matrix 00133 * 00134 SSMIN = HA 00135 SSMAX = FA 00136 CLT = ONE 00137 CRT = ONE 00138 SLT = ZERO 00139 SRT = ZERO 00140 ELSE 00141 GASMAL = .TRUE. 00142 IF( GA.GT.FA ) THEN 00143 PMAX = 2 00144 IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN 00145 * 00146 * Case of very large GA 00147 * 00148 GASMAL = .FALSE. 00149 SSMAX = GA 00150 IF( HA.GT.ONE ) THEN 00151 SSMIN = FA / ( GA / HA ) 00152 ELSE 00153 SSMIN = ( FA / GA )*HA 00154 END IF 00155 CLT = ONE 00156 SLT = HT / GT 00157 SRT = ONE 00158 CRT = FT / GT 00159 END IF 00160 END IF 00161 IF( GASMAL ) THEN 00162 * 00163 * Normal case 00164 * 00165 D = FA - HA 00166 IF( D.EQ.FA ) THEN 00167 * 00168 * Copes with infinite F or H 00169 * 00170 L = ONE 00171 ELSE 00172 L = D / FA 00173 END IF 00174 * 00175 * Note that 0 .le. L .le. 1 00176 * 00177 M = GT / FT 00178 * 00179 * Note that abs(M) .le. 1/macheps 00180 * 00181 T = TWO - L 00182 * 00183 * Note that T .ge. 1 00184 * 00185 MM = M*M 00186 TT = T*T 00187 S = SQRT( TT+MM ) 00188 * 00189 * Note that 1 .le. S .le. 1 + 1/macheps 00190 * 00191 IF( L.EQ.ZERO ) THEN 00192 R = ABS( M ) 00193 ELSE 00194 R = SQRT( L*L+MM ) 00195 END IF 00196 * 00197 * Note that 0 .le. R .le. 1 + 1/macheps 00198 * 00199 A = HALF*( S+R ) 00200 * 00201 * Note that 1 .le. A .le. 1 + abs(M) 00202 * 00203 SSMIN = HA / A 00204 SSMAX = FA*A 00205 IF( MM.EQ.ZERO ) THEN 00206 * 00207 * Note that M is very tiny 00208 * 00209 IF( L.EQ.ZERO ) THEN 00210 T = SIGN( TWO, FT )*SIGN( ONE, GT ) 00211 ELSE 00212 T = GT / SIGN( D, FT ) + M / T 00213 END IF 00214 ELSE 00215 T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A ) 00216 END IF 00217 L = SQRT( T*T+FOUR ) 00218 CRT = TWO / L 00219 SRT = T / L 00220 CLT = ( CRT+SRT*M ) / A 00221 SLT = ( HT / FT )*SRT / A 00222 END IF 00223 END IF 00224 IF( SWAP ) THEN 00225 CSL = SRT 00226 SNL = CRT 00227 CSR = SLT 00228 SNR = CLT 00229 ELSE 00230 CSL = CLT 00231 SNL = SLT 00232 CSR = CRT 00233 SNR = SRT 00234 END IF 00235 * 00236 * Correct signs of SSMAX and SSMIN 00237 * 00238 IF( PMAX.EQ.1 ) 00239 $ TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F ) 00240 IF( PMAX.EQ.2 ) 00241 $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G ) 00242 IF( PMAX.EQ.3 ) 00243 $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H ) 00244 SSMAX = SIGN( SSMAX, TSIGN ) 00245 SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) ) 00246 RETURN 00247 * 00248 * End of DLASV2 00249 * 00250 END
1.7.3 
