1 /* ./src_f77/dlags2.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
dlags2_(logical * upper,doublereal * a1,doublereal * a2,doublereal * a3,doublereal * b1,doublereal * b2,doublereal * b3,doublereal * csu,doublereal * snu,doublereal * csv,doublereal * snv,doublereal * csq,doublereal * snq)8 /* Subroutine */ int dlags2_(logical *upper, doublereal *a1, doublereal *a2,
9 doublereal *a3, doublereal *b1, doublereal *b2, doublereal *b3,
10 doublereal *csu, doublereal *snu, doublereal *csv, doublereal *snv,
11 doublereal *csq, doublereal *snq)
12 {
13 /* System generated locals */
14 doublereal d__1;
15
16 /* Local variables */
17 static doublereal a, b, c__, d__, r__, s1, s2, ua11, ua12, ua21, ua22,
18 vb11, vb12, vb21, vb22, csl, csr, snl, snr, aua11, aua12, aua21,
19 aua22, avb11, avb12, avb21, avb22, ua11r, ua22r, vb11r, vb22r;
20 extern /* Subroutine */ int dlasv2_(doublereal *, doublereal *,
21 doublereal *, doublereal *, doublereal *, doublereal *,
22 doublereal *, doublereal *, doublereal *), dlartg_(doublereal *,
23 doublereal *, doublereal *, doublereal *, doublereal *);
24
25
26 /* -- LAPACK auxiliary routine (version 3.0) -- */
27 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
28 /* Courant Institute, Argonne National Lab, and Rice University */
29 /* September 30, 1994 */
30
31 /* .. Scalar Arguments .. */
32 /* .. */
33
34 /* Purpose */
35 /* ======= */
36
37 /* DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such */
38 /* that if ( UPPER ) then */
39
40 /* U'*A*Q = U'*( A1 A2 )*Q = ( x 0 ) */
41 /* ( 0 A3 ) ( x x ) */
42 /* and */
43 /* V'*B*Q = V'*( B1 B2 )*Q = ( x 0 ) */
44 /* ( 0 B3 ) ( x x ) */
45
46 /* or if ( .NOT.UPPER ) then */
47
48 /* U'*A*Q = U'*( A1 0 )*Q = ( x x ) */
49 /* ( A2 A3 ) ( 0 x ) */
50 /* and */
51 /* V'*B*Q = V'*( B1 0 )*Q = ( x x ) */
52 /* ( B2 B3 ) ( 0 x ) */
53
54 /* The rows of the transformed A and B are parallel, where */
55
56 /* U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ ) */
57 /* ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ ) */
58
59 /* Z' denotes the transpose of Z. */
60
61
62 /* Arguments */
63 /* ========= */
64
65 /* UPPER (input) LOGICAL */
66 /* = .TRUE.: the input matrices A and B are upper triangular. */
67 /* = .FALSE.: the input matrices A and B are lower triangular. */
68
69 /* A1 (input) DOUBLE PRECISION */
70 /* A2 (input) DOUBLE PRECISION */
71 /* A3 (input) DOUBLE PRECISION */
72 /* On entry, A1, A2 and A3 are elements of the input 2-by-2 */
73 /* upper (lower) triangular matrix A. */
74
75 /* B1 (input) DOUBLE PRECISION */
76 /* B2 (input) DOUBLE PRECISION */
77 /* B3 (input) DOUBLE PRECISION */
78 /* On entry, B1, B2 and B3 are elements of the input 2-by-2 */
79 /* upper (lower) triangular matrix B. */
80
81 /* CSU (output) DOUBLE PRECISION */
82 /* SNU (output) DOUBLE PRECISION */
83 /* The desired orthogonal matrix U. */
84
85 /* CSV (output) DOUBLE PRECISION */
86 /* SNV (output) DOUBLE PRECISION */
87 /* The desired orthogonal matrix V. */
88
89 /* CSQ (output) DOUBLE PRECISION */
90 /* SNQ (output) DOUBLE PRECISION */
91 /* The desired orthogonal matrix Q. */
92
93 /* ===================================================================== */
94
95 /* .. Parameters .. */
96 /* .. */
97 /* .. Local Scalars .. */
98 /* .. */
99 /* .. External Subroutines .. */
100 /* .. */
101 /* .. Intrinsic Functions .. */
102 /* .. */
103 /* .. Executable Statements .. */
104
105 if (*upper) {
106
107 /* Input matrices A and B are upper triangular matrices */
108
109 /* Form matrix C = A*adj(B) = ( a b ) */
110 /* ( 0 d ) */
111
112 a = *a1 * *b3;
113 d__ = *a3 * *b1;
114 b = *a2 * *b1 - *a1 * *b2;
115
116 /* The SVD of real 2-by-2 triangular C */
117
118 /* ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 ) */
119 /* ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T ) */
120
121 dlasv2_(&a, &b, &d__, &s1, &s2, &snr, &csr, &snl, &csl);
122
123 if (abs(csl) >= abs(snl) || abs(csr) >= abs(snr)) {
124
125 /* Compute the (1,1) and (1,2) elements of U'*A and V'*B, */
126 /* and (1,2) element of |U|'*|A| and |V|'*|B|. */
127
128 ua11r = csl * *a1;
129 ua12 = csl * *a2 + snl * *a3;
130
131 vb11r = csr * *b1;
132 vb12 = csr * *b2 + snr * *b3;
133
134 aua12 = abs(csl) * abs(*a2) + abs(snl) * abs(*a3);
135 avb12 = abs(csr) * abs(*b2) + abs(snr) * abs(*b3);
136
137 /* zero (1,2) elements of U'*A and V'*B */
138
139 if (abs(ua11r) + abs(ua12) != 0.) {
140 if (aua12 / (abs(ua11r) + abs(ua12)) <= avb12 / (abs(vb11r) +
141 abs(vb12))) {
142 d__1 = -ua11r;
143 dlartg_(&d__1, &ua12, csq, snq, &r__);
144 } else {
145 d__1 = -vb11r;
146 dlartg_(&d__1, &vb12, csq, snq, &r__);
147 }
148 } else {
149 d__1 = -vb11r;
150 dlartg_(&d__1, &vb12, csq, snq, &r__);
151 }
152
153 *csu = csl;
154 *snu = -snl;
155 *csv = csr;
156 *snv = -snr;
157
158 } else {
159
160 /* Compute the (2,1) and (2,2) elements of U'*A and V'*B, */
161 /* and (2,2) element of |U|'*|A| and |V|'*|B|. */
162
163 ua21 = -snl * *a1;
164 ua22 = -snl * *a2 + csl * *a3;
165
166 vb21 = -snr * *b1;
167 vb22 = -snr * *b2 + csr * *b3;
168
169 aua22 = abs(snl) * abs(*a2) + abs(csl) * abs(*a3);
170 avb22 = abs(snr) * abs(*b2) + abs(csr) * abs(*b3);
171
172 /* zero (2,2) elements of U'*A and V'*B, and then swap. */
173
174 if (abs(ua21) + abs(ua22) != 0.) {
175 if (aua22 / (abs(ua21) + abs(ua22)) <= avb22 / (abs(vb21) +
176 abs(vb22))) {
177 d__1 = -ua21;
178 dlartg_(&d__1, &ua22, csq, snq, &r__);
179 } else {
180 d__1 = -vb21;
181 dlartg_(&d__1, &vb22, csq, snq, &r__);
182 }
183 } else {
184 d__1 = -vb21;
185 dlartg_(&d__1, &vb22, csq, snq, &r__);
186 }
187
188 *csu = snl;
189 *snu = csl;
190 *csv = snr;
191 *snv = csr;
192
193 }
194
195 } else {
196
197 /* Input matrices A and B are lower triangular matrices */
198
199 /* Form matrix C = A*adj(B) = ( a 0 ) */
200 /* ( c d ) */
201
202 a = *a1 * *b3;
203 d__ = *a3 * *b1;
204 c__ = *a2 * *b3 - *a3 * *b2;
205
206 /* The SVD of real 2-by-2 triangular C */
207
208 /* ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 ) */
209 /* ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T ) */
210
211 dlasv2_(&a, &c__, &d__, &s1, &s2, &snr, &csr, &snl, &csl);
212
213 if (abs(csr) >= abs(snr) || abs(csl) >= abs(snl)) {
214
215 /* Compute the (2,1) and (2,2) elements of U'*A and V'*B, */
216 /* and (2,1) element of |U|'*|A| and |V|'*|B|. */
217
218 ua21 = -snr * *a1 + csr * *a2;
219 ua22r = csr * *a3;
220
221 vb21 = -snl * *b1 + csl * *b2;
222 vb22r = csl * *b3;
223
224 aua21 = abs(snr) * abs(*a1) + abs(csr) * abs(*a2);
225 avb21 = abs(snl) * abs(*b1) + abs(csl) * abs(*b2);
226
227 /* zero (2,1) elements of U'*A and V'*B. */
228
229 if (abs(ua21) + abs(ua22r) != 0.) {
230 if (aua21 / (abs(ua21) + abs(ua22r)) <= avb21 / (abs(vb21) +
231 abs(vb22r))) {
232 dlartg_(&ua22r, &ua21, csq, snq, &r__);
233 } else {
234 dlartg_(&vb22r, &vb21, csq, snq, &r__);
235 }
236 } else {
237 dlartg_(&vb22r, &vb21, csq, snq, &r__);
238 }
239
240 *csu = csr;
241 *snu = -snr;
242 *csv = csl;
243 *snv = -snl;
244
245 } else {
246
247 /* Compute the (1,1) and (1,2) elements of U'*A and V'*B, */
248 /* and (1,1) element of |U|'*|A| and |V|'*|B|. */
249
250 ua11 = csr * *a1 + snr * *a2;
251 ua12 = snr * *a3;
252
253 vb11 = csl * *b1 + snl * *b2;
254 vb12 = snl * *b3;
255
256 aua11 = abs(csr) * abs(*a1) + abs(snr) * abs(*a2);
257 avb11 = abs(csl) * abs(*b1) + abs(snl) * abs(*b2);
258
259 /* zero (1,1) elements of U'*A and V'*B, and then swap. */
260
261 if (abs(ua11) + abs(ua12) != 0.) {
262 if (aua11 / (abs(ua11) + abs(ua12)) <= avb11 / (abs(vb11) +
263 abs(vb12))) {
264 dlartg_(&ua12, &ua11, csq, snq, &r__);
265 } else {
266 dlartg_(&vb12, &vb11, csq, snq, &r__);
267 }
268 } else {
269 dlartg_(&vb12, &vb11, csq, snq, &r__);
270 }
271
272 *csu = snr;
273 *snu = csr;
274 *csv = snl;
275 *snv = csl;
276
277 }
278
279 }
280
281 return 0;
282
283 /* End of DLAGS2 */
284
285 } /* dlags2_ */
286
287