1 SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
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
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
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