1 SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
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, 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
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, 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