001:       SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL,
002:      $                   CSR, SNR )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            LDA, LDB
010:       DOUBLE PRECISION   CSL, CSR, SNL, SNR
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
014:      $                   B( LDB, * ), BETA( 2 )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLAGV2 computes the Generalized Schur factorization of a real 2-by-2
021: *  matrix pencil (A,B) where B is upper triangular. This routine
022: *  computes orthogonal (rotation) matrices given by CSL, SNL and CSR,
023: *  SNR such that
024: *
025: *  1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0
026: *     types), then
027: *
028: *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
029: *     [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
030: *
031: *     [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
032: *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],
033: *
034: *  2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,
035: *     then
036: *
037: *     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]
038: *     [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]
039: *
040: *     [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]
041: *     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]
042: *
043: *     where b11 >= b22 > 0.
044: *
045: *
046: *  Arguments
047: *  =========
048: *
049: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, 2)
050: *          On entry, the 2 x 2 matrix A.
051: *          On exit, A is overwritten by the ``A-part'' of the
052: *          generalized Schur form.
053: *
054: *  LDA     (input) INTEGER
055: *          THe leading dimension of the array A.  LDA >= 2.
056: *
057: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, 2)
058: *          On entry, the upper triangular 2 x 2 matrix B.
059: *          On exit, B is overwritten by the ``B-part'' of the
060: *          generalized Schur form.
061: *
062: *  LDB     (input) INTEGER
063: *          THe leading dimension of the array B.  LDB >= 2.
064: *
065: *  ALPHAR  (output) DOUBLE PRECISION array, dimension (2)
066: *  ALPHAI  (output) DOUBLE PRECISION array, dimension (2)
067: *  BETA    (output) DOUBLE PRECISION array, dimension (2)
068: *          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the
069: *          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may
070: *          be zero.
071: *
072: *  CSL     (output) DOUBLE PRECISION
073: *          The cosine of the left rotation matrix.
074: *
075: *  SNL     (output) DOUBLE PRECISION
076: *          The sine of the left rotation matrix.
077: *
078: *  CSR     (output) DOUBLE PRECISION
079: *          The cosine of the right rotation matrix.
080: *
081: *  SNR     (output) DOUBLE PRECISION
082: *          The sine of the right rotation matrix.
083: *
084: *  Further Details
085: *  ===============
086: *
087: *  Based on contributions by
088: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
089: *
090: *  =====================================================================
091: *
092: *     .. Parameters ..
093:       DOUBLE PRECISION   ZERO, ONE
094:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
095: *     ..
096: *     .. Local Scalars ..
097:       DOUBLE PRECISION   ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
098:      $                   R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
099:      $                   WR2
100: *     ..
101: *     .. External Subroutines ..
102:       EXTERNAL           DLAG2, DLARTG, DLASV2, DROT
103: *     ..
104: *     .. External Functions ..
105:       DOUBLE PRECISION   DLAMCH, DLAPY2
106:       EXTERNAL           DLAMCH, DLAPY2
107: *     ..
108: *     .. Intrinsic Functions ..
109:       INTRINSIC          ABS, MAX
110: *     ..
111: *     .. Executable Statements ..
112: *
113:       SAFMIN = DLAMCH( 'S' )
114:       ULP = DLAMCH( 'P' )
115: *
116: *     Scale A
117: *
118:       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
119:      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
120:       ASCALE = ONE / ANORM
121:       A( 1, 1 ) = ASCALE*A( 1, 1 )
122:       A( 1, 2 ) = ASCALE*A( 1, 2 )
123:       A( 2, 1 ) = ASCALE*A( 2, 1 )
124:       A( 2, 2 ) = ASCALE*A( 2, 2 )
125: *
126: *     Scale B
127: *
128:       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
129:      $        SAFMIN )
130:       BSCALE = ONE / BNORM
131:       B( 1, 1 ) = BSCALE*B( 1, 1 )
132:       B( 1, 2 ) = BSCALE*B( 1, 2 )
133:       B( 2, 2 ) = BSCALE*B( 2, 2 )
134: *
135: *     Check if A can be deflated
136: *
137:       IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN
138:          CSL = ONE
139:          SNL = ZERO
140:          CSR = ONE
141:          SNR = ZERO
142:          A( 2, 1 ) = ZERO
143:          B( 2, 1 ) = ZERO
144: *
145: *     Check if B is singular
146: *
147:       ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN
148:          CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
149:          CSR = ONE
150:          SNR = ZERO
151:          CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
152:          CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
153:          A( 2, 1 ) = ZERO
154:          B( 1, 1 ) = ZERO
155:          B( 2, 1 ) = ZERO
156: *
157:       ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN
158:          CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T )
159:          SNR = -SNR
160:          CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
161:          CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
162:          CSL = ONE
163:          SNL = ZERO
164:          A( 2, 1 ) = ZERO
165:          B( 2, 1 ) = ZERO
166:          B( 2, 2 ) = ZERO
167: *
168:       ELSE
169: *
170: *        B is nonsingular, first compute the eigenvalues of (A,B)
171: *
172:          CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2,
173:      $               WI )
174: *
175:          IF( WI.EQ.ZERO ) THEN
176: *
177: *           two real eigenvalues, compute s*A-w*B
178: *
179:             H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 )
180:             H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 )
181:             H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 )
182: *
183:             RR = DLAPY2( H1, H2 )
184:             QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 )
185: *
186:             IF( RR.GT.QQ ) THEN
187: *
188: *              find right rotation matrix to zero 1,1 element of
189: *              (sA - wB)
190: *
191:                CALL DLARTG( H2, H1, CSR, SNR, T )
192: *
193:             ELSE
194: *
195: *              find right rotation matrix to zero 2,1 element of
196: *              (sA - wB)
197: *
198:                CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T )
199: *
200:             END IF
201: *
202:             SNR = -SNR
203:             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
204:             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
205: *
206: *           compute inf norms of A and B
207: *
208:             H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ),
209:      $           ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) )
210:             H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
211:      $           ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
212: *
213:             IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN
214: *
215: *              find left rotation matrix Q to zero out B(2,1)
216: *
217:                CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R )
218: *
219:             ELSE
220: *
221: *              find left rotation matrix Q to zero out A(2,1)
222: *
223:                CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R )
224: *
225:             END IF
226: *
227:             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
228:             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
229: *
230:             A( 2, 1 ) = ZERO
231:             B( 2, 1 ) = ZERO
232: *
233:          ELSE
234: *
235: *           a pair of complex conjugate eigenvalues
236: *           first compute the SVD of the matrix B
237: *
238:             CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR,
239:      $                   CSR, SNL, CSL )
240: *
241: *           Form (A,B) := Q(A,B)Z' where Q is left rotation matrix and
242: *           Z is right rotation matrix computed from DLASV2
243: *
244:             CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL )
245:             CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL )
246:             CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR )
247:             CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR )
248: *
249:             B( 2, 1 ) = ZERO
250:             B( 1, 2 ) = ZERO
251: *
252:          END IF
253: *
254:       END IF
255: *
256: *     Unscaling
257: *
258:       A( 1, 1 ) = ANORM*A( 1, 1 )
259:       A( 2, 1 ) = ANORM*A( 2, 1 )
260:       A( 1, 2 ) = ANORM*A( 1, 2 )
261:       A( 2, 2 ) = ANORM*A( 2, 2 )
262:       B( 1, 1 ) = BNORM*B( 1, 1 )
263:       B( 2, 1 ) = BNORM*B( 2, 1 )
264:       B( 1, 2 ) = BNORM*B( 1, 2 )
265:       B( 2, 2 ) = BNORM*B( 2, 2 )
266: *
267:       IF( WI.EQ.ZERO ) THEN
268:          ALPHAR( 1 ) = A( 1, 1 )
269:          ALPHAR( 2 ) = A( 2, 2 )
270:          ALPHAI( 1 ) = ZERO
271:          ALPHAI( 2 ) = ZERO
272:          BETA( 1 ) = B( 1, 1 )
273:          BETA( 2 ) = B( 2, 2 )
274:       ELSE
275:          ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM
276:          ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM
277:          ALPHAR( 2 ) = ALPHAR( 1 )
278:          ALPHAI( 2 ) = -ALPHAI( 1 )
279:          BETA( 1 ) = ONE
280:          BETA( 2 ) = ONE
281:       END IF
282: *
283:       RETURN
284: *
285: *     End of DLAGV2
286: *
287:       END
288: