001:       SUBROUTINE ZLAGS2( 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:       DOUBLE PRECISION   A1, A3, B1, B3, CSQ, CSU, CSV
012:       COMPLEX*16         A2, B2, SNQ, SNU, SNV
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  ZLAGS2 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) DOUBLE PRECISION
059: *  A2      (input) COMPLEX*16
060: *  A3      (input) DOUBLE PRECISION
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) DOUBLE PRECISION
065: *  B2      (input) COMPLEX*16
066: *  B3      (input) DOUBLE PRECISION
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) DOUBLE PRECISION
071: *  SNU     (output) COMPLEX*16
072: *          The desired unitary matrix U.
073: *
074: *  CSV     (output) DOUBLE PRECISION
075: *  SNV     (output) COMPLEX*16
076: *          The desired unitary matrix V.
077: *
078: *  CSQ     (output) DOUBLE PRECISION
079: *  SNQ     (output) COMPLEX*16
080: *          The desired unitary matrix Q.
081: *
082: *  =====================================================================
083: *
084: *     .. Parameters ..
085:       DOUBLE PRECISION   ZERO, ONE
086:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
087: *     ..
088: *     .. Local Scalars ..
089:       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11, 
090:      $                   AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2, 
091:      $                   SNL, SNR, UA11R, UA22R, VB11R, VB22R
092:       COMPLEX*16         B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
093:      $                   VB12, VB21, VB22
094: *     ..
095: *     .. External Subroutines ..
096:       EXTERNAL           DLASV2, ZLARTG
097: *     ..
098: *     .. Intrinsic Functions ..
099:       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
310: