1       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
  2      $                   CSR, SNR )
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.2.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 *     June 2010
  8 *
  9 *     .. Scalar Arguments ..
 10       INTEGER            LDA, LDB
 11       DOUBLE PRECISION   CSL, CSR, SNL, SNR
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
 15      $                   B( LDB, * ), BETA( 2 )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
 22 *  matrix pencil (A,B) where B is upper triangular. This routine
 23 *  computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
 24 *  SNR such that
 25 *
 26 *  1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
 27 *     types), then
 28 *
 29 *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
 30 *     [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
 31 *
 32 *     [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
 33 *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
 34 *
 35 *  2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
 36 *     then
 37 *
 38 *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
 39 *     [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
 40 *
 41 *     [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
 42 *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
 43 *
 44 *     where b11 >= b22 > 0.
 45 *
 46 *
 47 *  Arguments
 48 *  =========
 49 *
 50 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, 2)
 51 *          On entry, the 2 x 2 matrix A.
 52 *          On exit, A is overwritten by the ``A-part'' of the
 53 *          generalized Schur form.
 54 *
 55 *  LDA     (input) INTEGER
 56 *          THe leading dimension of the array A.  LDA >= 2.
 57 *
 58 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, 2)
 59 *          On entry, the upper triangular 2 x 2 matrix B.
 60 *          On exit, B is overwritten by the ``B-part'' of the
 61 *          generalized Schur form.
 62 *
 63 *  LDB     (input) INTEGER
 64 *          THe leading dimension of the array B.  LDB >= 2.
 65 *
 66 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (2)
 67 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (2)
 68 *  BETA    (output) DOUBLE PRECISION array, dimension (2)
 69 *          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
 70 *          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
 71 *          be zero.
 72 *
 73 *  CSL     (output) DOUBLE PRECISION
 74 *          The cosine of the left rotation matrix.
 75 *
 76 *  SNL     (output) DOUBLE PRECISION
 77 *          The sine of the left rotation matrix.
 78 *
 79 *  CSR     (output) DOUBLE PRECISION
 80 *          The cosine of the right rotation matrix.
 81 *
 82 *  SNR     (output) DOUBLE PRECISION
 83 *          The sine of the right rotation matrix.
 84 *
 85 *  Further Details
 86 *  ===============
 87 *
 88 *  Based on contributions by
 89 *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 90 *
 91 *  =====================================================================
 92 *
 93 *     .. Parameters ..
 94       DOUBLE PRECISION   ZERO, ONE
 95       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 96 *     ..
 97 *     .. Local Scalars ..
 98       DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
 99      $                   R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
100      $                   WR2
101 *     ..
102 *     .. External Subroutines ..
103       EXTERNAL           DLAG2, DLARTG, DLASV2, DROT
104 *     ..
105 *     .. External Functions ..
106       DOUBLE PRECISION   DLAMCH, DLAPY2
107       EXTERNAL           DLAMCH, DLAPY2
108 *     ..
109 *     .. Intrinsic Functions ..
110       INTRINSIC          ABSMAX
111 *     ..
112 *     .. Executable Statements ..
113 *
114       SAFMIN = DLAMCH( 'S' )
115       ULP = DLAMCH( 'P' )
116 *
117 *     Scale A
118 *
119       ANORM = MAXABS( A( 11 ) )+ABS( A( 21 ) ),
120      $        ABS( A( 12 ) )+ABS( A( 22 ) ), SAFMIN )
121       ASCALE = ONE / ANORM
122       A( 11 ) = ASCALE*A( 11 )
123       A( 12 ) = ASCALE*A( 12 )
124       A( 21 ) = ASCALE*A( 21 )
125       A( 22 ) = ASCALE*A( 22 )
126 *
127 *     Scale B
128 *
129       BNORM = MAXABS( B( 11 ) ), ABS( B( 12 ) )+ABS( B( 22 ) ),
130      $        SAFMIN )
131       BSCALE = ONE / BNORM
132       B( 11 ) = BSCALE*B( 11 )
133       B( 12 ) = BSCALE*B( 12 )
134       B( 22 ) = BSCALE*B( 22 )
135 *
136 *     Check if A can be deflated
137 *
138       IFABS( A( 21 ) ).LE.ULP ) THEN
139          CSL = ONE
140          SNL = ZERO
141          CSR = ONE
142          SNR = ZERO
143          A( 21 ) = ZERO
144          B( 21 ) = ZERO
145          WI = ZERO
146 *
147 *     Check if B is singular
148 *
149       ELSE IFABS( B( 11 ) ).LE.ULP ) THEN
150          CALL DLARTG( A( 11 ), A( 21 ), CSL, SNL, R )
151          CSR = ONE
152          SNR = ZERO
153          CALL DROT( 2, A( 11 ), LDA, A( 21 ), LDA, CSL, SNL )
154          CALL DROT( 2, B( 11 ), LDB, B( 21 ), LDB, CSL, SNL )
155          A( 21 ) = ZERO
156          B( 11 ) = ZERO
157          B( 21 ) = ZERO
158          WI = ZERO
159 *
160       ELSE IFABS( B( 22 ) ).LE.ULP ) THEN
161          CALL DLARTG( A( 22 ), A( 21 ), CSR, SNR, T )
162          SNR = -SNR
163          CALL DROT( 2, A( 11 ), 1, A( 12 ), 1, CSR, SNR )
164          CALL DROT( 2, B( 11 ), 1, B( 12 ), 1, CSR, SNR )
165          CSL = ONE
166          SNL = ZERO
167          A( 21 ) = ZERO
168          B( 21 ) = ZERO
169          B( 22 ) = ZERO
170          WI = ZERO
171 *
172       ELSE
173 *
174 *        B is nonsingular, first compute the eigenvalues of (A,B)
175 *
176          CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
177      $               WI )
178 *
179          IF( WI.EQ.ZERO ) THEN
180 *
181 *           two real eigenvalues, compute s*A-w*B
182 *
183             H1 = SCALE1*A( 11 ) - WR1*B( 11 )
184             H2 = SCALE1*A( 12 ) - WR1*B( 12 )
185             H3 = SCALE1*A( 22 ) - WR1*B( 22 )
186 *
187             RR = DLAPY2( H1, H2 )
188             QQ = DLAPY2( SCALE1*A( 21 ), H3 )
189 *
190             IF( RR.GT.QQ ) THEN
191 *
192 *              find right rotation matrix to zero 1,1 element of
193 *              (sA - wB)
194 *
195                CALL DLARTG( H2, H1, CSR, SNR, T )
196 *
197             ELSE
198 *
199 *              find right rotation matrix to zero 2,1 element of
200 *              (sA - wB)
201 *
202                CALL DLARTG( H3, SCALE1*A( 21 ), CSR, SNR, T )
203 *
204             END IF
205 *
206             SNR = -SNR
207             CALL DROT( 2, A( 11 ), 1, A( 12 ), 1, CSR, SNR )
208             CALL DROT( 2, B( 11 ), 1, B( 12 ), 1, CSR, SNR )
209 *
210 *           compute inf norms of A and B
211 *
212             H1 = MAXABS( A( 11 ) )+ABS( A( 12 ) ),
213      $           ABS( A( 21 ) )+ABS( A( 22 ) ) )
214             H2 = MAXABS( B( 11 ) )+ABS( B( 12 ) ),
215      $           ABS( B( 21 ) )+ABS( B( 22 ) ) )
216 *
217             IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
218 *
219 *              find left rotation matrix Q to zero out B(2,1)
220 *
221                CALL DLARTG( B( 11 ), B( 21 ), CSL, SNL, R )
222 *
223             ELSE
224 *
225 *              find left rotation matrix Q to zero out A(2,1)
226 *
227                CALL DLARTG( A( 11 ), A( 21 ), CSL, SNL, R )
228 *
229             END IF
230 *
231             CALL DROT( 2, A( 11 ), LDA, A( 21 ), LDA, CSL, SNL )
232             CALL DROT( 2, B( 11 ), LDB, B( 21 ), LDB, CSL, SNL )
233 *
234             A( 21 ) = ZERO
235             B( 21 ) = ZERO
236 *
237          ELSE
238 *
239 *           a pair of complex conjugate eigenvalues
240 *           first compute the SVD of the matrix B
241 *
242             CALL DLASV2( B( 11 ), B( 12 ), B( 22 ), R, T, SNR,
243      $                   CSR, SNL, CSL )
244 *
245 *           Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and
246 *           Z is right rotation matrix computed from DLASV2
247 *
248             CALL DROT( 2, A( 11 ), LDA, A( 21 ), LDA, CSL, SNL )
249             CALL DROT( 2, B( 11 ), LDB, B( 21 ), LDB, CSL, SNL )
250             CALL DROT( 2, A( 11 ), 1, A( 12 ), 1, CSR, SNR )
251             CALL DROT( 2, B( 11 ), 1, B( 12 ), 1, CSR, SNR )
252 *
253             B( 21 ) = ZERO
254             B( 12 ) = ZERO
255 *
256          END IF
257 *
258       END IF
259 *
260 *     Unscaling
261 *
262       A( 11 ) = ANORM*A( 11 )
263       A( 21 ) = ANORM*A( 21 )
264       A( 12 ) = ANORM*A( 12 )
265       A( 22 ) = ANORM*A( 22 )
266       B( 11 ) = BNORM*B( 11 )
267       B( 21 ) = BNORM*B( 21 )
268       B( 12 ) = BNORM*B( 12 )
269       B( 22 ) = BNORM*B( 22 )
270 *
271       IF( WI.EQ.ZERO ) THEN
272          ALPHAR( 1 ) = A( 11 )
273          ALPHAR( 2 ) = A( 22 )
274          ALPHAI( 1 ) = ZERO
275          ALPHAI( 2 ) = ZERO
276          BETA( 1 ) = B( 11 )
277          BETA( 2 ) = B( 22 )
278       ELSE
279          ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
280          ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
281          ALPHAR( 2 ) = ALPHAR( 1 )
282          ALPHAI( 2 ) = -ALPHAI( 1 )
283          BETA( 1 ) = ONE
284          BETA( 2 ) = ONE
285       END IF
286 *
287       RETURN
288 *
289 *     End of DLAGV2
290 *
291       END