![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlanv2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlanv2.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
2: *
3: * -- LAPACK driver routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
10: * ..
11: *
12: * Purpose
13: * =======
14: *
15: * DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
16: * matrix in standard form:
17: *
18: * [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
19: * [ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
20: *
21: * where either
22: * 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
23: * 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
24: * conjugate eigenvalues.
25: *
26: * Arguments
27: * =========
28: *
29: * A (input/output) DOUBLE PRECISION
30: * B (input/output) DOUBLE PRECISION
31: * C (input/output) DOUBLE PRECISION
32: * D (input/output) DOUBLE PRECISION
33: * On entry, the elements of the input matrix.
34: * On exit, they are overwritten by the elements of the
35: * standardised Schur form.
36: *
37: * RT1R (output) DOUBLE PRECISION
38: * RT1I (output) DOUBLE PRECISION
39: * RT2R (output) DOUBLE PRECISION
40: * RT2I (output) DOUBLE PRECISION
41: * The real and imaginary parts of the eigenvalues. If the
42: * eigenvalues are a complex conjugate pair, RT1I > 0.
43: *
44: * CS (output) DOUBLE PRECISION
45: * SN (output) DOUBLE PRECISION
46: * Parameters of the rotation matrix.
47: *
48: * Further Details
49: * ===============
50: *
51: * Modified by V. Sima, Research Institute for Informatics, Bucharest,
52: * Romania, to reduce the risk of cancellation errors,
53: * when computing real eigenvalues, and to ensure, if possible, that
54: * abs(RT1R) >= abs(RT2R).
55: *
56: * =====================================================================
57: *
58: * .. Parameters ..
59: DOUBLE PRECISION ZERO, HALF, ONE
60: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
61: DOUBLE PRECISION MULTPL
62: PARAMETER ( MULTPL = 4.0D+0 )
63: * ..
64: * .. Local Scalars ..
65: DOUBLE PRECISION AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
66: $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
67: * ..
68: * .. External Functions ..
69: DOUBLE PRECISION DLAMCH, DLAPY2
70: EXTERNAL DLAMCH, DLAPY2
71: * ..
72: * .. Intrinsic Functions ..
73: INTRINSIC ABS, MAX, MIN, SIGN, SQRT
74: * ..
75: * .. Executable Statements ..
76: *
77: EPS = DLAMCH( 'P' )
78: IF( C.EQ.ZERO ) THEN
79: CS = ONE
80: SN = ZERO
81: GO TO 10
82: *
83: ELSE IF( B.EQ.ZERO ) THEN
84: *
85: * Swap rows and columns
86: *
87: CS = ZERO
88: SN = ONE
89: TEMP = D
90: D = A
91: A = TEMP
92: B = -C
93: C = ZERO
94: GO TO 10
95: ELSE IF( ( A-D ).EQ.ZERO .AND. SIGN( ONE, B ).NE.SIGN( ONE, C ) )
96: $ THEN
97: CS = ONE
98: SN = ZERO
99: 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 = DLAPY2( 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 = DLAPY2( 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 DLANV2
205: *
206: END