|
2814
|
1 SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) |
|
|
2 * |
|
3333
|
3 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2814
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
|
6 * October 31, 1992 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 |
|
|
10 * .. |
|
|
11 * |
|
|
12 * Purpose |
|
|
13 * ======= |
|
|
14 * |
|
|
15 * DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix |
|
|
16 * [ A B ] |
|
|
17 * [ B C ]. |
|
|
18 * On return, RT1 is the eigenvalue of larger absolute value, RT2 is the |
|
|
19 * eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right |
|
|
20 * eigenvector for RT1, giving the decomposition |
|
|
21 * |
|
|
22 * [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] |
|
|
23 * [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. |
|
|
24 * |
|
|
25 * Arguments |
|
|
26 * ========= |
|
|
27 * |
|
|
28 * A (input) DOUBLE PRECISION |
|
|
29 * The (1,1) element of the 2-by-2 matrix. |
|
|
30 * |
|
|
31 * B (input) DOUBLE PRECISION |
|
|
32 * The (1,2) element and the conjugate of the (2,1) element of |
|
|
33 * the 2-by-2 matrix. |
|
|
34 * |
|
|
35 * C (input) DOUBLE PRECISION |
|
|
36 * The (2,2) element of the 2-by-2 matrix. |
|
|
37 * |
|
|
38 * RT1 (output) DOUBLE PRECISION |
|
|
39 * The eigenvalue of larger absolute value. |
|
|
40 * |
|
|
41 * RT2 (output) DOUBLE PRECISION |
|
|
42 * The eigenvalue of smaller absolute value. |
|
|
43 * |
|
|
44 * CS1 (output) DOUBLE PRECISION |
|
|
45 * SN1 (output) DOUBLE PRECISION |
|
|
46 * The vector (CS1, SN1) is a unit right eigenvector for RT1. |
|
|
47 * |
|
|
48 * Further Details |
|
|
49 * =============== |
|
|
50 * |
|
|
51 * RT1 is accurate to a few ulps barring over/underflow. |
|
|
52 * |
|
|
53 * RT2 may be inaccurate if there is massive cancellation in the |
|
|
54 * determinant A*C-B*B; higher precision or correctly rounded or |
|
|
55 * correctly truncated arithmetic would be needed to compute RT2 |
|
|
56 * accurately in all cases. |
|
|
57 * |
|
|
58 * CS1 and SN1 are accurate to a few ulps barring over/underflow. |
|
|
59 * |
|
|
60 * Overflow is possible only if RT1 is within a factor of 5 of overflow. |
|
|
61 * Underflow is harmless if the input data is 0 or exceeds |
|
|
62 * underflow_threshold / macheps. |
|
|
63 * |
|
|
64 * ===================================================================== |
|
|
65 * |
|
|
66 * .. Parameters .. |
|
|
67 DOUBLE PRECISION ONE |
|
|
68 PARAMETER ( ONE = 1.0D0 ) |
|
|
69 DOUBLE PRECISION TWO |
|
|
70 PARAMETER ( TWO = 2.0D0 ) |
|
|
71 DOUBLE PRECISION ZERO |
|
|
72 PARAMETER ( ZERO = 0.0D0 ) |
|
|
73 DOUBLE PRECISION HALF |
|
|
74 PARAMETER ( HALF = 0.5D0 ) |
|
|
75 * .. |
|
|
76 * .. Local Scalars .. |
|
|
77 INTEGER SGN1, SGN2 |
|
|
78 DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, |
|
|
79 $ TB, TN |
|
|
80 * .. |
|
|
81 * .. Intrinsic Functions .. |
|
|
82 INTRINSIC ABS, SQRT |
|
|
83 * .. |
|
|
84 * .. Executable Statements .. |
|
|
85 * |
|
|
86 * Compute the eigenvalues |
|
|
87 * |
|
|
88 SM = A + C |
|
|
89 DF = A - C |
|
|
90 ADF = ABS( DF ) |
|
|
91 TB = B + B |
|
|
92 AB = ABS( TB ) |
|
|
93 IF( ABS( A ).GT.ABS( C ) ) THEN |
|
|
94 ACMX = A |
|
|
95 ACMN = C |
|
|
96 ELSE |
|
|
97 ACMX = C |
|
|
98 ACMN = A |
|
|
99 END IF |
|
|
100 IF( ADF.GT.AB ) THEN |
|
|
101 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) |
|
|
102 ELSE IF( ADF.LT.AB ) THEN |
|
|
103 RT = AB*SQRT( ONE+( ADF / AB )**2 ) |
|
|
104 ELSE |
|
|
105 * |
|
|
106 * Includes case AB=ADF=0 |
|
|
107 * |
|
|
108 RT = AB*SQRT( TWO ) |
|
|
109 END IF |
|
|
110 IF( SM.LT.ZERO ) THEN |
|
|
111 RT1 = HALF*( SM-RT ) |
|
|
112 SGN1 = -1 |
|
|
113 * |
|
|
114 * Order of execution important. |
|
|
115 * To get fully accurate smaller eigenvalue, |
|
|
116 * next line needs to be executed in higher precision. |
|
|
117 * |
|
|
118 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B |
|
|
119 ELSE IF( SM.GT.ZERO ) THEN |
|
|
120 RT1 = HALF*( SM+RT ) |
|
|
121 SGN1 = 1 |
|
|
122 * |
|
|
123 * Order of execution important. |
|
|
124 * To get fully accurate smaller eigenvalue, |
|
|
125 * next line needs to be executed in higher precision. |
|
|
126 * |
|
|
127 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B |
|
|
128 ELSE |
|
|
129 * |
|
|
130 * Includes case RT1 = RT2 = 0 |
|
|
131 * |
|
|
132 RT1 = HALF*RT |
|
|
133 RT2 = -HALF*RT |
|
|
134 SGN1 = 1 |
|
|
135 END IF |
|
|
136 * |
|
|
137 * Compute the eigenvector |
|
|
138 * |
|
|
139 IF( DF.GE.ZERO ) THEN |
|
|
140 CS = DF + RT |
|
|
141 SGN2 = 1 |
|
|
142 ELSE |
|
|
143 CS = DF - RT |
|
|
144 SGN2 = -1 |
|
|
145 END IF |
|
|
146 ACS = ABS( CS ) |
|
|
147 IF( ACS.GT.AB ) THEN |
|
|
148 CT = -TB / CS |
|
|
149 SN1 = ONE / SQRT( ONE+CT*CT ) |
|
|
150 CS1 = CT*SN1 |
|
|
151 ELSE |
|
|
152 IF( AB.EQ.ZERO ) THEN |
|
|
153 CS1 = ONE |
|
|
154 SN1 = ZERO |
|
|
155 ELSE |
|
|
156 TN = -CS / TB |
|
|
157 CS1 = ONE / SQRT( ONE+TN*TN ) |
|
|
158 SN1 = TN*CS1 |
|
|
159 END IF |
|
|
160 END IF |
|
|
161 IF( SGN1.EQ.SGN2 ) THEN |
|
|
162 TN = CS1 |
|
|
163 CS1 = -SN1 |
|
|
164 SN1 = TN |
|
|
165 END IF |
|
|
166 RETURN |
|
|
167 * |
|
|
168 * End of DLAEV2 |
|
|
169 * |
|
|
170 END |
