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          ABSDBLEDCMPLXDCONJGDIMAG
 98 *     ..
 99 *     .. Statement Functions ..
100       DOUBLE PRECISION   ABS1
101 *     ..
102 *     .. Statement Function definitions ..
103       ABS1( T ) = ABSDBLE( T ) ) + ABSDIMAG( 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          IFABS( 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          IFABS( 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