Annotation of rpl/lapack/lapack/dlags2.f, revision 1.13
1.12 bertrand 1: *> \brief \b DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAGS2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlags2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlags2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlags2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
22: * SNV, CSQ, SNQ )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL UPPER
26: * DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
27: * $ SNU, SNV
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
37: *> that if ( UPPER ) then
38: *>
39: *> U**T *A*Q = U**T *( A1 A2 )*Q = ( x 0 )
40: *> ( 0 A3 ) ( x x )
41: *> and
42: *> V**T*B*Q = V**T *( B1 B2 )*Q = ( x 0 )
43: *> ( 0 B3 ) ( x x )
44: *>
45: *> or if ( .NOT.UPPER ) then
46: *>
47: *> U**T *A*Q = U**T *( A1 0 )*Q = ( x x )
48: *> ( A2 A3 ) ( 0 x )
49: *> and
50: *> V**T*B*Q = V**T*( B1 0 )*Q = ( x x )
51: *> ( B2 B3 ) ( 0 x )
52: *>
53: *> The rows of the transformed A and B are parallel, where
54: *>
55: *> U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
56: *> ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
57: *>
58: *> Z**T denotes the transpose of Z.
59: *>
60: *> \endverbatim
61: *
62: * Arguments:
63: * ==========
64: *
65: *> \param[in] UPPER
66: *> \verbatim
67: *> UPPER is LOGICAL
68: *> = .TRUE.: the input matrices A and B are upper triangular.
69: *> = .FALSE.: the input matrices A and B are lower triangular.
70: *> \endverbatim
71: *>
72: *> \param[in] A1
73: *> \verbatim
74: *> A1 is DOUBLE PRECISION
75: *> \endverbatim
76: *>
77: *> \param[in] A2
78: *> \verbatim
79: *> A2 is DOUBLE PRECISION
80: *> \endverbatim
81: *>
82: *> \param[in] A3
83: *> \verbatim
84: *> A3 is DOUBLE PRECISION
85: *> On entry, A1, A2 and A3 are elements of the input 2-by-2
86: *> upper (lower) triangular matrix A.
87: *> \endverbatim
88: *>
89: *> \param[in] B1
90: *> \verbatim
91: *> B1 is DOUBLE PRECISION
92: *> \endverbatim
93: *>
94: *> \param[in] B2
95: *> \verbatim
96: *> B2 is DOUBLE PRECISION
97: *> \endverbatim
98: *>
99: *> \param[in] B3
100: *> \verbatim
101: *> B3 is DOUBLE PRECISION
102: *> On entry, B1, B2 and B3 are elements of the input 2-by-2
103: *> upper (lower) triangular matrix B.
104: *> \endverbatim
105: *>
106: *> \param[out] CSU
107: *> \verbatim
108: *> CSU is DOUBLE PRECISION
109: *> \endverbatim
110: *>
111: *> \param[out] SNU
112: *> \verbatim
113: *> SNU is DOUBLE PRECISION
114: *> The desired orthogonal matrix U.
115: *> \endverbatim
116: *>
117: *> \param[out] CSV
118: *> \verbatim
119: *> CSV is DOUBLE PRECISION
120: *> \endverbatim
121: *>
122: *> \param[out] SNV
123: *> \verbatim
124: *> SNV is DOUBLE PRECISION
125: *> The desired orthogonal matrix V.
126: *> \endverbatim
127: *>
128: *> \param[out] CSQ
129: *> \verbatim
130: *> CSQ is DOUBLE PRECISION
131: *> \endverbatim
132: *>
133: *> \param[out] SNQ
134: *> \verbatim
135: *> SNQ is DOUBLE PRECISION
136: *> The desired orthogonal matrix Q.
137: *> \endverbatim
138: *
139: * Authors:
140: * ========
141: *
142: *> \author Univ. of Tennessee
143: *> \author Univ. of California Berkeley
144: *> \author Univ. of Colorado Denver
145: *> \author NAG Ltd.
146: *
1.12 bertrand 147: *> \date September 2012
1.9 bertrand 148: *
149: *> \ingroup doubleOTHERauxiliary
150: *
151: * =====================================================================
1.1 bertrand 152: SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
153: $ SNV, CSQ, SNQ )
154: *
1.12 bertrand 155: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 156: * -- LAPACK is a software package provided by Univ. of Tennessee, --
157: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 bertrand 158: * September 2012
1.1 bertrand 159: *
160: * .. Scalar Arguments ..
161: LOGICAL UPPER
162: DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
163: $ SNU, SNV
164: * ..
165: *
166: * =====================================================================
167: *
168: * .. Parameters ..
169: DOUBLE PRECISION ZERO
170: PARAMETER ( ZERO = 0.0D+0 )
171: * ..
172: * .. Local Scalars ..
173: DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
174: $ AVB21, AVB22, B, C, CSL, CSR, D, R, S1, S2,
175: $ SNL, SNR, UA11, UA11R, UA12, UA21, UA22, UA22R,
176: $ VB11, VB11R, VB12, VB21, VB22, VB22R
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL DLARTG, DLASV2
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC ABS
183: * ..
184: * .. Executable Statements ..
185: *
186: IF( UPPER ) THEN
187: *
188: * Input matrices A and B are upper triangular matrices
189: *
190: * Form matrix C = A*adj(B) = ( a b )
191: * ( 0 d )
192: *
193: A = A1*B3
194: D = A3*B1
195: B = A2*B1 - A1*B2
196: *
197: * The SVD of real 2-by-2 triangular C
198: *
199: * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
200: * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
201: *
202: CALL DLASV2( A, B, D, S1, S2, SNR, CSR, SNL, CSL )
203: *
204: IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
205: $ THEN
206: *
1.8 bertrand 207: * Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
208: * and (1,2) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 209: *
210: UA11R = CSL*A1
211: UA12 = CSL*A2 + SNL*A3
212: *
213: VB11R = CSR*B1
214: VB12 = CSR*B2 + SNR*B3
215: *
216: AUA12 = ABS( CSL )*ABS( A2 ) + ABS( SNL )*ABS( A3 )
217: AVB12 = ABS( CSR )*ABS( B2 ) + ABS( SNR )*ABS( B3 )
218: *
1.8 bertrand 219: * zero (1,2) elements of U**T *A and V**T *B
1.1 bertrand 220: *
221: IF( ( ABS( UA11R )+ABS( UA12 ) ).NE.ZERO ) THEN
222: IF( AUA12 / ( ABS( UA11R )+ABS( UA12 ) ).LE.AVB12 /
223: $ ( ABS( VB11R )+ABS( VB12 ) ) ) THEN
224: CALL DLARTG( -UA11R, UA12, CSQ, SNQ, R )
225: ELSE
226: CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
227: END IF
228: ELSE
229: CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
230: END IF
231: *
232: CSU = CSL
233: SNU = -SNL
234: CSV = CSR
235: SNV = -SNR
236: *
237: ELSE
238: *
1.8 bertrand 239: * Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
240: * and (2,2) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 241: *
242: UA21 = -SNL*A1
243: UA22 = -SNL*A2 + CSL*A3
244: *
245: VB21 = -SNR*B1
246: VB22 = -SNR*B2 + CSR*B3
247: *
248: AUA22 = ABS( SNL )*ABS( A2 ) + ABS( CSL )*ABS( A3 )
249: AVB22 = ABS( SNR )*ABS( B2 ) + ABS( CSR )*ABS( B3 )
250: *
1.8 bertrand 251: * zero (2,2) elements of U**T*A and V**T*B, and then swap.
1.1 bertrand 252: *
253: IF( ( ABS( UA21 )+ABS( UA22 ) ).NE.ZERO ) THEN
254: IF( AUA22 / ( ABS( UA21 )+ABS( UA22 ) ).LE.AVB22 /
255: $ ( ABS( VB21 )+ABS( VB22 ) ) ) THEN
256: CALL DLARTG( -UA21, UA22, CSQ, SNQ, R )
257: ELSE
258: CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
259: END IF
260: ELSE
261: CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
262: END IF
263: *
264: CSU = SNL
265: SNU = CSL
266: CSV = SNR
267: SNV = CSR
268: *
269: END IF
270: *
271: ELSE
272: *
273: * Input matrices A and B are lower triangular matrices
274: *
275: * Form matrix C = A*adj(B) = ( a 0 )
276: * ( c d )
277: *
278: A = A1*B3
279: D = A3*B1
280: C = A2*B3 - A3*B2
281: *
282: * The SVD of real 2-by-2 triangular C
283: *
284: * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
285: * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
286: *
287: CALL DLASV2( A, C, D, S1, S2, SNR, CSR, SNL, CSL )
288: *
289: IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
290: $ THEN
291: *
1.8 bertrand 292: * Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
293: * and (2,1) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 294: *
295: UA21 = -SNR*A1 + CSR*A2
296: UA22R = CSR*A3
297: *
298: VB21 = -SNL*B1 + CSL*B2
299: VB22R = CSL*B3
300: *
301: AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS( A2 )
302: AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS( B2 )
303: *
1.8 bertrand 304: * zero (2,1) elements of U**T *A and V**T *B.
1.1 bertrand 305: *
306: IF( ( ABS( UA21 )+ABS( UA22R ) ).NE.ZERO ) THEN
307: IF( AUA21 / ( ABS( UA21 )+ABS( UA22R ) ).LE.AVB21 /
308: $ ( ABS( VB21 )+ABS( VB22R ) ) ) THEN
309: CALL DLARTG( UA22R, UA21, CSQ, SNQ, R )
310: ELSE
311: CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
312: END IF
313: ELSE
314: CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
315: END IF
316: *
317: CSU = CSR
318: SNU = -SNR
319: CSV = CSL
320: SNV = -SNL
321: *
322: ELSE
323: *
1.8 bertrand 324: * Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
325: * and (1,1) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 326: *
327: UA11 = CSR*A1 + SNR*A2
328: UA12 = SNR*A3
329: *
330: VB11 = CSL*B1 + SNL*B2
331: VB12 = SNL*B3
332: *
333: AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS( A2 )
334: AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS( B2 )
335: *
1.8 bertrand 336: * zero (1,1) elements of U**T*A and V**T*B, and then swap.
1.1 bertrand 337: *
338: IF( ( ABS( UA11 )+ABS( UA12 ) ).NE.ZERO ) THEN
339: IF( AUA11 / ( ABS( UA11 )+ABS( UA12 ) ).LE.AVB11 /
340: $ ( ABS( VB11 )+ABS( VB12 ) ) ) THEN
341: CALL DLARTG( UA12, UA11, CSQ, SNQ, R )
342: ELSE
343: CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
344: END IF
345: ELSE
346: CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
347: END IF
348: *
349: CSU = SNR
350: SNU = CSR
351: CSV = SNL
352: SNV = CSL
353: *
354: END IF
355: *
356: END IF
357: *
358: RETURN
359: *
360: * End of DLAGS2
361: *
362: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlags2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack