|
2329
|
1 SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK driver routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN |
|
|
9 * .. |
|
|
10 * |
|
|
11 * Purpose |
|
|
12 * ======= |
|
|
13 * |
|
|
14 * DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric |
|
|
15 * matrix in standard form: |
|
|
16 * |
|
|
17 * [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ] |
|
|
18 * [ C D ] [ SN CS ] [ CC DD ] [-SN CS ] |
|
|
19 * |
|
|
20 * where either |
|
|
21 * 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or |
|
|
22 * 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex |
|
|
23 * conjugate eigenvalues. |
|
|
24 * |
|
|
25 * Arguments |
|
|
26 * ========= |
|
|
27 * |
|
|
28 * A (input/output) DOUBLE PRECISION |
|
|
29 * B (input/output) DOUBLE PRECISION |
|
|
30 * C (input/output) DOUBLE PRECISION |
|
|
31 * D (input/output) DOUBLE PRECISION |
|
|
32 * On entry, the elements of the input matrix. |
|
|
33 * On exit, they are overwritten by the elements of the |
|
|
34 * standardised Schur form. |
|
|
35 * |
|
|
36 * RT1R (output) DOUBLE PRECISION |
|
|
37 * RT1I (output) DOUBLE PRECISION |
|
|
38 * RT2R (output) DOUBLE PRECISION |
|
|
39 * RT2I (output) DOUBLE PRECISION |
|
|
40 * The real and imaginary parts of the eigenvalues. If the |
|
|
41 * eigenvalues are a complex conjugate pair, RT1I > 0. |
|
|
42 * |
|
|
43 * CS (output) DOUBLE PRECISION |
|
|
44 * SN (output) DOUBLE PRECISION |
|
|
45 * Parameters of the rotation matrix. |
|
|
46 * |
|
3333
|
47 * Further Details |
|
|
48 * =============== |
|
|
49 * |
|
|
50 * Modified by V. Sima, Research Institute for Informatics, Bucharest, |
|
|
51 * Romania, to reduce the risk of cancellation errors, |
|
|
52 * when computing real eigenvalues, and to ensure, if possible, that |
|
|
53 * abs(RT1R) >= abs(RT2R). |
|
|
54 * |
|
2329
|
55 * ===================================================================== |
|
|
56 * |
|
|
57 * .. Parameters .. |
|
|
58 DOUBLE PRECISION ZERO, HALF, ONE |
|
|
59 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) |
|
3333
|
60 DOUBLE PRECISION MULTPL |
|
|
61 PARAMETER ( MULTPL = 4.0D+0 ) |
|
2329
|
62 * .. |
|
|
63 * .. Local Scalars .. |
|
3333
|
64 DOUBLE PRECISION AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB, |
|
|
65 $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z |
|
2329
|
66 * .. |
|
|
67 * .. External Functions .. |
|
3333
|
68 DOUBLE PRECISION DLAMCH, DLAPY2 |
|
|
69 EXTERNAL DLAMCH, DLAPY2 |
|
2329
|
70 * .. |
|
|
71 * .. Intrinsic Functions .. |
|
3333
|
72 INTRINSIC ABS, MAX, MIN, SIGN, SQRT |
|
2329
|
73 * .. |
|
|
74 * .. Executable Statements .. |
|
|
75 * |
|
3333
|
76 EPS = DLAMCH( 'P' ) |
|
2329
|
77 IF( C.EQ.ZERO ) THEN |
|
3333
|
78 CS = ONE |
|
|
79 SN = ZERO |
|
2329
|
80 GO TO 10 |
|
|
81 * |
|
|
82 ELSE IF( B.EQ.ZERO ) THEN |
|
|
83 * |
|
|
84 * Swap rows and columns |
|
|
85 * |
|
|
86 CS = ZERO |
|
|
87 SN = ONE |
|
|
88 TEMP = D |
|
|
89 D = A |
|
|
90 A = TEMP |
|
|
91 B = -C |
|
|
92 C = ZERO |
|
|
93 GO TO 10 |
|
3333
|
94 ELSE IF( ( A-D ).EQ.ZERO .AND. SIGN( ONE, B ).NE.SIGN( ONE, C ) ) |
|
|
95 $ THEN |
|
|
96 CS = ONE |
|
|
97 SN = ZERO |
|
2329
|
98 GO TO 10 |
|
|
99 ELSE |
|
|
100 * |
|
|
101 TEMP = A - D |
|
|
102 P = HALF*TEMP |
|
3333
|
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 |
|
2329
|
110 * |
|
3333
|
111 IF( Z.GE.MULTPL*EPS ) THEN |
|
|
112 * |
|
|
113 * Real eigenvalues. Compute A and D. |
|
2329
|
114 * |
|
3333
|
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 |
|
2329
|
120 * |
|
3333
|
121 TAU = DLAPY2( C, Z ) |
|
|
122 CS = Z / TAU |
|
|
123 SN = C / TAU |
|
|
124 B = B - C |
|
|
125 C = ZERO |
|
|
126 ELSE |
|
2329
|
127 * |
|
3333
|
128 * Complex eigenvalues, or real (almost) equal eigenvalues. |
|
|
129 * Make diagonal elements equal. |
|
2329
|
130 * |
|
3333
|
131 SIGMA = B + C |
|
|
132 TAU = DLAPY2( SIGMA, TEMP ) |
|
|
133 CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) ) |
|
|
134 SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA ) |
|
2329
|
135 * |
|
3333
|
136 * Compute [ AA BB ] = [ A B ] [ CS -SN ] |
|
|
137 * [ CC DD ] [ C D ] [ SN CS ] |
|
2329
|
138 * |
|
3333
|
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 ] |
|
2329
|
146 * |
|
3333
|
147 A = AA*CS + CC*SN |
|
|
148 B = BB*CS + DD*SN |
|
|
149 C = -AA*SN + CC*CS |
|
|
150 D = -BB*SN + DD*CS |
|
2329
|
151 * |
|
3333
|
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 |
|
2329
|
159 * |
|
3333
|
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 |
|
2329
|
178 C = ZERO |
|
3333
|
179 TEMP = CS |
|
|
180 CS = -SN |
|
|
181 SN = TEMP |
|
2329
|
182 END IF |
|
|
183 END IF |
|
|
184 END IF |
|
3333
|
185 * |
|
2329
|
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 DLANV2 |
|
|
204 * |
|
|
205 END |
