001:       SUBROUTINE SLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
002: *
003: *  -- LAPACK driver routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       REAL               A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
015: *  matrix in standard form:
016: *
017: *       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
018: *       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
019: *
020: *  where either
021: *  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
022: *  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
023: *  conjugate eigenvalues.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  A       (input/output) REAL            
029: *  B       (input/output) REAL            
030: *  C       (input/output) REAL            
031: *  D       (input/output) REAL            
032: *          On entry, the elements of the input matrix.
033: *          On exit, they are overwritten by the elements of the
034: *          standardised Schur form.
035: *
036: *  RT1R    (output) REAL 
037: *  RT1I    (output) REAL            
038: *  RT2R    (output) REAL            
039: *  RT2I    (output) REAL            
040: *          The real and imaginary parts of the eigenvalues. If the
041: *          eigenvalues are a complex conjugate pair, RT1I > 0.
042: *
043: *  CS      (output) REAL            
044: *  SN      (output) REAL            
045: *          Parameters of the rotation matrix.
046: *
047: *  Further Details
048: *  ===============
049: *
050: *  Modified by V. Sima, Research Institute for Informatics, Bucharest,
051: *  Romania, to reduce the risk of cancellation errors,
052: *  when computing real eigenvalues, and to ensure, if possible, that
053: *  abs(RT1R) >= abs(RT2R).
054: *
055: *  =====================================================================
056: *
057: *     .. Parameters ..
058:       REAL               ZERO, HALF, ONE
059:       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
060:       REAL               MULTPL
061:       PARAMETER          ( MULTPL = 4.0E+0 )
062: *     ..
063: *     .. Local Scalars ..
064:       REAL               AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
065:      $                   SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
066: *     ..
067: *     .. External Functions ..
068:       REAL               SLAMCH, SLAPY2
069:       EXTERNAL           SLAMCH, SLAPY2
070: *     ..
071: *     .. Intrinsic Functions ..
072:       INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
073: *     ..
074: *     .. Executable Statements ..
075: *
076:       EPS = SLAMCH( 'P' )
077:       IF( C.EQ.ZERO ) THEN
078:          CS = ONE
079:          SN = ZERO
080:          GO TO 10
081: *
082:       ELSE IF( B.EQ.ZERO ) THEN
083: *
084: *        Swap rows and columns
085: *
086:          CS = ZERO
087:          SN = ONE
088:          TEMP = D
089:          D = A
090:          A = TEMP
091:          B = -C
092:          C = ZERO
093:          GO TO 10
094:       ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
095:      $   SIGN( ONE, C ) ) THEN
096:          CS = ONE
097:          SN = ZERO
098:          GO TO 10
099:       ELSE
100: *
101:          TEMP = A - D
102:          P = HALF*TEMP
103:          BCMAX = MAX( ABS( B ), ABS( C ) )
104:          BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C )
105:          SCALE = MAX( ABS( P ), BCMAX )
106:          Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
107: *
108: *        If Z is of the order of the machine accuracy, postpone the
109: *        decision on the nature of eigenvalues
110: *
111:          IF( Z.GE.MULTPL*EPS ) THEN
112: *
113: *           Real eigenvalues. Compute A and D.
114: *
115:             Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
116:             A = D + Z
117:             D = D - ( BCMAX / Z )*BCMIS
118: *
119: *           Compute B and the rotation matrix
120: *
121:             TAU = SLAPY2( C, Z )
122:             CS = Z / TAU
123:             SN = C / TAU
124:             B = B - C
125:             C = ZERO
126:          ELSE
127: *
128: *           Complex eigenvalues, or real (almost) equal eigenvalues.
129: *           Make diagonal elements equal.
130: *
131:             SIGMA = B + C
132:             TAU = SLAPY2( SIGMA, TEMP )
133:             CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
134:             SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
135: *
136: *           Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
137: *                   [ CC  DD ]   [ C  D ] [ SN  CS ]
138: *
139:             AA = A*CS + B*SN
140:             BB = -A*SN + B*CS
141:             CC = C*CS + D*SN
142:             DD = -C*SN + D*CS
143: *
144: *           Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
145: *                   [ C  D ]   [-SN  CS ] [ CC  DD ]
146: *
147:             A = AA*CS + CC*SN
148:             B = BB*CS + DD*SN
149:             C = -AA*SN + CC*CS
150:             D = -BB*SN + DD*CS
151: *
152:             TEMP = HALF*( A+D )
153:             A = TEMP
154:             D = TEMP
155: *
156:             IF( C.NE.ZERO ) THEN
157:                IF( B.NE.ZERO ) THEN
158:                   IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
159: *
160: *                    Real eigenvalues: reduce to upper triangular form
161: *
162:                      SAB = SQRT( ABS( B ) )
163:                      SAC = SQRT( ABS( C ) )
164:                      P = SIGN( SAB*SAC, C )
165:                      TAU = ONE / SQRT( ABS( B+C ) )
166:                      A = TEMP + P
167:                      D = TEMP - P
168:                      B = B - C
169:                      C = ZERO
170:                      CS1 = SAB*TAU
171:                      SN1 = SAC*TAU
172:                      TEMP = CS*CS1 - SN*SN1
173:                      SN = CS*SN1 + SN*CS1
174:                      CS = TEMP
175:                   END IF
176:                ELSE
177:                   B = -C
178:                   C = ZERO
179:                   TEMP = CS
180:                   CS = -SN
181:                   SN = TEMP
182:                END IF
183:             END IF
184:          END IF
185: *
186:       END IF
187: *
188:    10 CONTINUE
189: *
190: *     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
191: *
192:       RT1R = A
193:       RT2R = D
194:       IF( C.EQ.ZERO ) THEN
195:          RT1I = ZERO
196:          RT2I = ZERO
197:       ELSE
198:          RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
199:          RT2I = -RT1I
200:       END IF
201:       RETURN
202: *
203: *     End of SLANV2
204: *
205:       END
206: