1 SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
2 *
3 * -- LAPACK auxiliary 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, 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
2 *
3 * -- LAPACK auxiliary 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, 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