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