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