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