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