001:       SUBROUTINE CLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
002:      $                   SNV, CSQ, SNQ )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            UPPER
011:       REAL               A1, A3, B1, B3, CSQ, CSU, CSV
012:       COMPLEX            A2, B2, SNQ, SNU, SNV
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  CLAGS2 computes 2-by-2 unitary matrices U, V and Q, such
019: *  that if ( UPPER ) then
020: *
021: *            U'*A*Q = U'*( A1 A2 )*Q = ( x  0  )
022: *                        ( 0  A3 )     ( x  x  )
023: *  and
024: *            V'*B*Q = V'*( B1 B2 )*Q = ( x  0  )
025: *                        ( 0  B3 )     ( x  x  )
026: *
027: *  or if ( .NOT.UPPER ) then
028: *
029: *            U'*A*Q = U'*( A1 0  )*Q = ( x  x  )
030: *                        ( A2 A3 )     ( 0  x  )
031: *  and
032: *            V'*B*Q = V'*( B1 0  )*Q = ( x  x  )
033: *                        ( B2 B3 )     ( 0  x  )
034: *  where
035: *
036: *    U = (     CSU      SNU ), V = (     CSV     SNV ),
037: *        ( -CONJG(SNU)  CSU )      ( -CONJG(SNV) CSV )
038: *
039: *    Q = (     CSQ      SNQ )
040: *        ( -CONJG(SNQ)  CSQ )
041: *
042: *  Z' denotes the conjugate transpose of Z.
043: *
044: *  The rows of the transformed A and B are parallel. Moreover, if the
045: *  input 2-by-2 matrix A is not zero, then the transformed (1,1) entry
046: *  of A is not zero. If the input matrices A and B are both not zero,
047: *  then the transformed (2,2) element of B is not zero, except when the
048: *  first rows of input A and B are parallel and the second rows are
049: *  zero.
050: *
051: *  Arguments
052: *  =========
053: *
054: *  UPPER   (input) LOGICAL
055: *          = .TRUE.: the input matrices A and B are upper triangular.
056: *          = .FALSE.: the input matrices A and B are lower triangular.
057: *
058: *  A1      (input) REAL
059: *  A2      (input) COMPLEX
060: *  A3      (input) REAL
061: *          On entry, A1, A2 and A3 are elements of the input 2-by-2
062: *          upper (lower) triangular matrix A.
063: *
064: *  B1      (input) REAL
065: *  B2      (input) COMPLEX
066: *  B3      (input) REAL
067: *          On entry, B1, B2 and B3 are elements of the input 2-by-2
068: *          upper (lower) triangular matrix B.
069: *
070: *  CSU     (output) REAL
071: *  SNU     (output) COMPLEX
072: *          The desired unitary matrix U.
073: *
074: *  CSV     (output) REAL
075: *  SNV     (output) COMPLEX
076: *          The desired unitary matrix V.
077: *
078: *  CSQ     (output) REAL
079: *  SNQ     (output) COMPLEX
080: *          The desired unitary matrix Q.
081: *
082: *  =====================================================================
083: *
084: *     .. Parameters ..
085:       REAL               ZERO, ONE
086:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
087: *     ..
088: *     .. Local Scalars ..
089:       REAL               A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
090:      $                   AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2, SNL,
091:      $                   SNR, UA11R, UA22R, VB11R, VB22R
092:       COMPLEX            B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
093:      $                   VB12, VB21, VB22
094: *     ..
095: *     .. External Subroutines ..
096:       EXTERNAL           CLARTG, SLASV2
097: *     ..
098: *     .. Intrinsic Functions ..
099:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, REAL
100: *     ..
101: *     .. Statement Functions ..
102:       REAL               ABS1
103: *     ..
104: *     .. Statement Function definitions ..
105:       ABS1( T ) = ABS( REAL( T ) ) + ABS( AIMAG( 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 SLASV2( 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 CLARTG( -CMPLX( VB11R ), CONJG( VB12 ), CSQ, SNQ,
154:      $                      R )
155:             ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN
156:                CALL CLARTG( -CMPLX( UA11R ), CONJG( UA12 ), CSQ, SNQ,
157:      $                      R )
158:             ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 /
159:      $               ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN
160:                CALL CLARTG( -CMPLX( UA11R ), CONJG( UA12 ), CSQ, SNQ,
161:      $                      R )
162:             ELSE
163:                CALL CLARTG( -CMPLX( VB11R ), CONJG( 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 = -CONJG( D1 )*SNL*A1
178:             UA22 = -CONJG( D1 )*SNL*A2 + CSL*A3
179: *
180:             VB21 = -CONJG( D1 )*SNR*B1
181:             VB22 = -CONJG( 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 CLARTG( -CONJG( VB21 ), CONJG( VB22 ), CSQ, SNQ, R )
190:             ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN
191:                CALL CLARTG( -CONJG( UA21 ), CONJG( UA22 ), CSQ, SNQ, R )
192:             ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 /
193:      $               ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN
194:                CALL CLARTG( -CONJG( UA21 ), CONJG( UA22 ), CSQ, SNQ, R )
195:             ELSE
196:                CALL CLARTG( -CONJG( VB21 ), CONJG( VB22 ), CSQ, SNQ, R )
197:             END IF
198: *
199:             CSU = SNL
200:             SNU = D1*CSL
201:             CSV = SNR
202:             SNV = D1*CSR
203: *
204:          END IF
205: *
206:       ELSE
207: *
208: *        Input matrices A and B are lower triangular matrices
209: *
210: *        Form matrix C = A*adj(B) = ( a 0 )
211: *                                   ( c d )
212: *
213:          A = A1*B3
214:          D = A3*B1
215:          C = A2*B3 - A3*B2
216:          FC = ABS( C )
217: *
218: *        Transform complex 2-by-2 matrix C to real matrix by unitary
219: *        diagonal matrix diag(d1,1).
220: *
221:          D1 = ONE
222:          IF( FC.NE.ZERO )
223:      $      D1 = C / FC
224: *
225: *        The SVD of real 2 by 2 triangular C
226: *
227: *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
228: *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
229: *
230:          CALL SLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL )
231: *
232:          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
233:      $        THEN
234: *
235: *           Compute the (2,1) and (2,2) elements of U'*A and V'*B,
236: *           and (2,1) element of |U|'*|A| and |V|'*|B|.
237: *
238:             UA21 = -D1*SNR*A1 + CSR*A2
239:             UA22R = CSR*A3
240: *
241:             VB21 = -D1*SNL*B1 + CSL*B2
242:             VB22R = CSL*B3
243: *
244:             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 )
245:             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 )
246: *
247: *           zero (2,1) elements of U'*A and V'*B.
248: *
249:             IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN
250:                CALL CLARTG( CMPLX( VB22R ), VB21, CSQ, SNQ, R )
251:             ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN
252:                CALL CLARTG( CMPLX( UA22R ), UA21, CSQ, SNQ, R )
253:             ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 /
254:      $               ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN
255:                CALL CLARTG( CMPLX( UA22R ), UA21, CSQ, SNQ, R )
256:             ELSE
257:                CALL CLARTG( CMPLX( VB22R ), VB21, CSQ, SNQ, R )
258:             END IF
259: *
260:             CSU = CSR
261:             SNU = -CONJG( D1 )*SNR
262:             CSV = CSL
263:             SNV = -CONJG( D1 )*SNL
264: *
265:          ELSE
266: *
267: *           Compute the (1,1) and (1,2) elements of U'*A and V'*B,
268: *           and (1,1) element of |U|'*|A| and |V|'*|B|.
269: *
270:             UA11 = CSR*A1 + CONJG( D1 )*SNR*A2
271:             UA12 = CONJG( D1 )*SNR*A3
272: *
273:             VB11 = CSL*B1 + CONJG( D1 )*SNL*B2
274:             VB12 = CONJG( D1 )*SNL*B3
275: *
276:             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 )
277:             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 )
278: *
279: *           zero (1,1) elements of U'*A and V'*B, and then swap.
280: *
281:             IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN
282:                CALL CLARTG( VB12, VB11, CSQ, SNQ, R )
283:             ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN
284:                CALL CLARTG( UA12, UA11, CSQ, SNQ, R )
285:             ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 /
286:      $               ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN
287:                CALL CLARTG( UA12, UA11, CSQ, SNQ, R )
288:             ELSE
289:                CALL CLARTG( VB12, VB11, CSQ, SNQ, R )
290:             END IF
291: *
292:             CSU = SNR
293:             SNU = CONJG( D1 )*CSR
294:             CSV = SNL
295:             SNV = CONJG( D1 )*CSL
296: *
297:          END IF
298: *
299:       END IF
300: *
301:       RETURN
302: *
303: *     End of CLAGS2
304: *
305:       END
306: