|
2329
|
1 SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL ) |
|
|
2 * |
|
3333
|
3 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2329
|
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 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 |
