1: SUBROUTINE ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
2: $ SNV, CSQ, SNQ )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: LOGICAL UPPER
11: DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV
12: COMPLEX*16 A2, B2, SNQ, SNU, SNV
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * ZLAGS2 computes 2-by-2 unitary matrices U, V and Q, such
19: * that if ( UPPER ) then
20: *
21: * U'*A*Q = U'*( A1 A2 )*Q = ( x 0 )
22: * ( 0 A3 ) ( x x )
23: * and
24: * V'*B*Q = V'*( B1 B2 )*Q = ( x 0 )
25: * ( 0 B3 ) ( x x )
26: *
27: * or if ( .NOT.UPPER ) then
28: *
29: * U'*A*Q = U'*( A1 0 )*Q = ( x x )
30: * ( A2 A3 ) ( 0 x )
31: * and
32: * V'*B*Q = V'*( B1 0 )*Q = ( x x )
33: * ( B2 B3 ) ( 0 x )
34: * where
35: *
36: * U = ( CSU SNU ), V = ( CSV SNV ),
37: * ( -CONJG(SNU) CSU ) ( -CONJG(SNV) CSV )
38: *
39: * Q = ( CSQ SNQ )
40: * ( -CONJG(SNQ) CSQ )
41: *
42: * Z' denotes the conjugate transpose of Z.
43: *
44: * The rows of the transformed A and B are parallel. Moreover, if the
45: * input 2-by-2 matrix A is not zero, then the transformed (1,1) entry
46: * of A is not zero. If the input matrices A and B are both not zero,
47: * then the transformed (2,2) element of B is not zero, except when the
48: * first rows of input A and B are parallel and the second rows are
49: * zero.
50: *
51: * Arguments
52: * =========
53: *
54: * UPPER (input) LOGICAL
55: * = .TRUE.: the input matrices A and B are upper triangular.
56: * = .FALSE.: the input matrices A and B are lower triangular.
57: *
58: * A1 (input) DOUBLE PRECISION
59: * A2 (input) COMPLEX*16
60: * A3 (input) DOUBLE PRECISION
61: * On entry, A1, A2 and A3 are elements of the input 2-by-2
62: * upper (lower) triangular matrix A.
63: *
64: * B1 (input) DOUBLE PRECISION
65: * B2 (input) COMPLEX*16
66: * B3 (input) DOUBLE PRECISION
67: * On entry, B1, B2 and B3 are elements of the input 2-by-2
68: * upper (lower) triangular matrix B.
69: *
70: * CSU (output) DOUBLE PRECISION
71: * SNU (output) COMPLEX*16
72: * The desired unitary matrix U.
73: *
74: * CSV (output) DOUBLE PRECISION
75: * SNV (output) COMPLEX*16
76: * The desired unitary matrix V.
77: *
78: * CSQ (output) DOUBLE PRECISION
79: * SNQ (output) COMPLEX*16
80: * The desired unitary matrix Q.
81: *
82: * =====================================================================
83: *
84: * .. Parameters ..
85: DOUBLE PRECISION ZERO, ONE
86: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
87: * ..
88: * .. Local Scalars ..
89: DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11,
90: $ AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2,
91: $ SNL, SNR, UA11R, UA22R, VB11R, VB22R
92: COMPLEX*16 B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
93: $ VB12, VB21, VB22
94: * ..
95: * .. External Subroutines ..
96: EXTERNAL DLASV2, ZLARTG
97: * ..
98: * .. Intrinsic Functions ..
99: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG
100: * ..
101: * .. Statement Functions ..
102: DOUBLE PRECISION ABS1
103: * ..
104: * .. Statement Function definitions ..
105: ABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) )
106: * ..
107: * .. Executable Statements ..
108: *
109: IF( UPPER ) THEN
110: *
111: * Input matrices A and B are upper triangular matrices
112: *
113: * Form matrix C = A*adj(B) = ( a b )
114: * ( 0 d )
115: *
116: A = A1*B3
117: D = A3*B1
118: B = A2*B1 - A1*B2
119: FB = ABS( B )
120: *
121: * Transform complex 2-by-2 matrix C to real matrix by unitary
122: * diagonal matrix diag(1,D1).
123: *
124: D1 = ONE
125: IF( FB.NE.ZERO )
126: $ D1 = B / FB
127: *
128: * The SVD of real 2 by 2 triangular C
129: *
130: * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
131: * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
132: *
133: CALL DLASV2( A, FB, D, S1, S2, SNR, CSR, SNL, CSL )
134: *
135: IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
136: $ THEN
137: *
138: * Compute the (1,1) and (1,2) elements of U'*A and V'*B,
139: * and (1,2) element of |U|'*|A| and |V|'*|B|.
140: *
141: UA11R = CSL*A1
142: UA12 = CSL*A2 + D1*SNL*A3
143: *
144: VB11R = CSR*B1
145: VB12 = CSR*B2 + D1*SNR*B3
146: *
147: AUA12 = ABS( CSL )*ABS1( A2 ) + ABS( SNL )*ABS( A3 )
148: AVB12 = ABS( CSR )*ABS1( B2 ) + ABS( SNR )*ABS( B3 )
149: *
150: * zero (1,2) elements of U'*A and V'*B
151: *
152: IF( ( ABS( UA11R )+ABS1( UA12 ) ).EQ.ZERO ) THEN
153: CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
154: $ R )
155: ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN
156: CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
157: $ R )
158: ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 /
159: $ ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN
160: CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ,
161: $ R )
162: ELSE
163: CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ,
164: $ R )
165: END IF
166: *
167: CSU = CSL
168: SNU = -D1*SNL
169: CSV = CSR
170: SNV = -D1*SNR
171: *
172: ELSE
173: *
174: * Compute the (2,1) and (2,2) elements of U'*A and V'*B,
175: * and (2,2) element of |U|'*|A| and |V|'*|B|.
176: *
177: UA21 = -DCONJG( D1 )*SNL*A1
178: UA22 = -DCONJG( D1 )*SNL*A2 + CSL*A3
179: *
180: VB21 = -DCONJG( D1 )*SNR*B1
181: VB22 = -DCONJG( D1 )*SNR*B2 + CSR*B3
182: *
183: AUA22 = ABS( SNL )*ABS1( A2 ) + ABS( CSL )*ABS( A3 )
184: AVB22 = ABS( SNR )*ABS1( B2 ) + ABS( CSR )*ABS( B3 )
185: *
186: * zero (2,2) elements of U'*A and V'*B, and then swap.
187: *
188: IF( ( ABS1( UA21 )+ABS1( UA22 ) ).EQ.ZERO ) THEN
189: CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
190: $ R )
191: ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN
192: CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
193: $ R )
194: ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 /
195: $ ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN
196: CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ,
197: $ R )
198: ELSE
199: CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ,
200: $ R )
201: END IF
202: *
203: CSU = SNL
204: SNU = D1*CSL
205: CSV = SNR
206: SNV = D1*CSR
207: *
208: END IF
209: *
210: ELSE
211: *
212: * Input matrices A and B are lower triangular matrices
213: *
214: * Form matrix C = A*adj(B) = ( a 0 )
215: * ( c d )
216: *
217: A = A1*B3
218: D = A3*B1
219: C = A2*B3 - A3*B2
220: FC = ABS( C )
221: *
222: * Transform complex 2-by-2 matrix C to real matrix by unitary
223: * diagonal matrix diag(d1,1).
224: *
225: D1 = ONE
226: IF( FC.NE.ZERO )
227: $ D1 = C / FC
228: *
229: * The SVD of real 2 by 2 triangular C
230: *
231: * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
232: * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
233: *
234: CALL DLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL )
235: *
236: IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
237: $ THEN
238: *
239: * Compute the (2,1) and (2,2) elements of U'*A and V'*B,
240: * and (2,1) element of |U|'*|A| and |V|'*|B|.
241: *
242: UA21 = -D1*SNR*A1 + CSR*A2
243: UA22R = CSR*A3
244: *
245: VB21 = -D1*SNL*B1 + CSL*B2
246: VB22R = CSL*B3
247: *
248: AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 )
249: AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 )
250: *
251: * zero (2,1) elements of U'*A and V'*B.
252: *
253: IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN
254: CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
255: ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN
256: CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
257: ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 /
258: $ ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN
259: CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R )
260: ELSE
261: CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R )
262: END IF
263: *
264: CSU = CSR
265: SNU = -DCONJG( D1 )*SNR
266: CSV = CSL
267: SNV = -DCONJG( D1 )*SNL
268: *
269: ELSE
270: *
271: * Compute the (1,1) and (1,2) elements of U'*A and V'*B,
272: * and (1,1) element of |U|'*|A| and |V|'*|B|.
273: *
274: UA11 = CSR*A1 + DCONJG( D1 )*SNR*A2
275: UA12 = DCONJG( D1 )*SNR*A3
276: *
277: VB11 = CSL*B1 + DCONJG( D1 )*SNL*B2
278: VB12 = DCONJG( D1 )*SNL*B3
279: *
280: AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 )
281: AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 )
282: *
283: * zero (1,1) elements of U'*A and V'*B, and then swap.
284: *
285: IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN
286: CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
287: ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN
288: CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
289: ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 /
290: $ ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN
291: CALL ZLARTG( UA12, UA11, CSQ, SNQ, R )
292: ELSE
293: CALL ZLARTG( VB12, VB11, CSQ, SNQ, R )
294: END IF
295: *
296: CSU = SNR
297: SNU = DCONJG( D1 )*CSR
298: CSV = SNL
299: SNV = DCONJG( D1 )*CSL
300: *
301: END IF
302: *
303: END IF
304: *
305: RETURN
306: *
307: * End of ZLAGS2
308: *
309: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlags2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack