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