|
2814
|
1 SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) |
|
|
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, RT1, RT2 |
|
|
10 * .. |
|
|
11 * |
|
|
12 * Purpose |
|
|
13 * ======= |
|
|
14 * |
|
|
15 * DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix |
|
|
16 * [ A B ] |
|
|
17 * [ B C ]. |
|
|
18 * On return, RT1 is the eigenvalue of larger absolute value, and RT2 |
|
|
19 * is the eigenvalue of smaller absolute value. |
|
|
20 * |
|
|
21 * Arguments |
|
|
22 * ========= |
|
|
23 * |
|
|
24 * A (input) DOUBLE PRECISION |
|
|
25 * The (1,1) element of the 2-by-2 matrix. |
|
|
26 * |
|
|
27 * B (input) DOUBLE PRECISION |
|
|
28 * The (1,2) and (2,1) elements of the 2-by-2 matrix. |
|
|
29 * |
|
|
30 * C (input) DOUBLE PRECISION |
|
|
31 * The (2,2) element of the 2-by-2 matrix. |
|
|
32 * |
|
|
33 * RT1 (output) DOUBLE PRECISION |
|
|
34 * The eigenvalue of larger absolute value. |
|
|
35 * |
|
|
36 * RT2 (output) DOUBLE PRECISION |
|
|
37 * The eigenvalue of smaller absolute value. |
|
|
38 * |
|
|
39 * Further Details |
|
|
40 * =============== |
|
|
41 * |
|
|
42 * RT1 is accurate to a few ulps barring over/underflow. |
|
|
43 * |
|
|
44 * RT2 may be inaccurate if there is massive cancellation in the |
|
|
45 * determinant A*C-B*B; higher precision or correctly rounded or |
|
|
46 * correctly truncated arithmetic would be needed to compute RT2 |
|
|
47 * accurately in all cases. |
|
|
48 * |
|
|
49 * Overflow is possible only if RT1 is within a factor of 5 of overflow. |
|
|
50 * Underflow is harmless if the input data is 0 or exceeds |
|
|
51 * underflow_threshold / macheps. |
|
|
52 * |
|
|
53 * ===================================================================== |
|
|
54 * |
|
|
55 * .. Parameters .. |
|
|
56 DOUBLE PRECISION ONE |
|
|
57 PARAMETER ( ONE = 1.0D0 ) |
|
|
58 DOUBLE PRECISION TWO |
|
|
59 PARAMETER ( TWO = 2.0D0 ) |
|
|
60 DOUBLE PRECISION ZERO |
|
|
61 PARAMETER ( ZERO = 0.0D0 ) |
|
|
62 DOUBLE PRECISION HALF |
|
|
63 PARAMETER ( HALF = 0.5D0 ) |
|
|
64 * .. |
|
|
65 * .. Local Scalars .. |
|
|
66 DOUBLE PRECISION AB, ACMN, ACMX, ADF, DF, RT, SM, TB |
|
|
67 * .. |
|
|
68 * .. Intrinsic Functions .. |
|
|
69 INTRINSIC ABS, SQRT |
|
|
70 * .. |
|
|
71 * .. Executable Statements .. |
|
|
72 * |
|
|
73 * Compute the eigenvalues |
|
|
74 * |
|
|
75 SM = A + C |
|
|
76 DF = A - C |
|
|
77 ADF = ABS( DF ) |
|
|
78 TB = B + B |
|
|
79 AB = ABS( TB ) |
|
|
80 IF( ABS( A ).GT.ABS( C ) ) THEN |
|
|
81 ACMX = A |
|
|
82 ACMN = C |
|
|
83 ELSE |
|
|
84 ACMX = C |
|
|
85 ACMN = A |
|
|
86 END IF |
|
|
87 IF( ADF.GT.AB ) THEN |
|
|
88 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) |
|
|
89 ELSE IF( ADF.LT.AB ) THEN |
|
|
90 RT = AB*SQRT( ONE+( ADF / AB )**2 ) |
|
|
91 ELSE |
|
|
92 * |
|
|
93 * Includes case AB=ADF=0 |
|
|
94 * |
|
|
95 RT = AB*SQRT( TWO ) |
|
|
96 END IF |
|
|
97 IF( SM.LT.ZERO ) THEN |
|
|
98 RT1 = HALF*( SM-RT ) |
|
|
99 * |
|
|
100 * Order of execution important. |
|
|
101 * To get fully accurate smaller eigenvalue, |
|
|
102 * next line needs to be executed in higher precision. |
|
|
103 * |
|
|
104 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B |
|
|
105 ELSE IF( SM.GT.ZERO ) THEN |
|
|
106 RT1 = HALF*( SM+RT ) |
|
|
107 * |
|
|
108 * Order of execution important. |
|
|
109 * To get fully accurate smaller eigenvalue, |
|
|
110 * next line needs to be executed in higher precision. |
|
|
111 * |
|
|
112 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B |
|
|
113 ELSE |
|
|
114 * |
|
|
115 * Includes case RT1 = RT2 = 0 |
|
|
116 * |
|
|
117 RT1 = HALF*RT |
|
|
118 RT2 = -HALF*RT |
|
|
119 END IF |
|
|
120 RETURN |
|
|
121 * |
|
|
122 * End of DLAE2 |
|
|
123 * |
|
|
124 END |
