1: *> \brief \b DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAEV2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
22: *
23: * .. Scalar Arguments ..
24: * DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1
25: * ..
26: *
27: *
28: *> \par Purpose:
29: * =============
30: *>
31: *> \verbatim
32: *>
33: *> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
34: *> [ A B ]
35: *> [ B C ].
36: *> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
37: *> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
38: *> eigenvector for RT1, giving the decomposition
39: *>
40: *> [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
41: *> [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] A
48: *> \verbatim
49: *> A is DOUBLE PRECISION
50: *> The (1,1) element of the 2-by-2 matrix.
51: *> \endverbatim
52: *>
53: *> \param[in] B
54: *> \verbatim
55: *> B is DOUBLE PRECISION
56: *> The (1,2) element and the conjugate of the (2,1) element of
57: *> the 2-by-2 matrix.
58: *> \endverbatim
59: *>
60: *> \param[in] C
61: *> \verbatim
62: *> C is DOUBLE PRECISION
63: *> The (2,2) element of the 2-by-2 matrix.
64: *> \endverbatim
65: *>
66: *> \param[out] RT1
67: *> \verbatim
68: *> RT1 is DOUBLE PRECISION
69: *> The eigenvalue of larger absolute value.
70: *> \endverbatim
71: *>
72: *> \param[out] RT2
73: *> \verbatim
74: *> RT2 is DOUBLE PRECISION
75: *> The eigenvalue of smaller absolute value.
76: *> \endverbatim
77: *>
78: *> \param[out] CS1
79: *> \verbatim
80: *> CS1 is DOUBLE PRECISION
81: *> \endverbatim
82: *>
83: *> \param[out] SN1
84: *> \verbatim
85: *> SN1 is DOUBLE PRECISION
86: *> The vector (CS1, SN1) is a unit right eigenvector for RT1.
87: *> \endverbatim
88: *
89: * Authors:
90: * ========
91: *
92: *> \author Univ. of Tennessee
93: *> \author Univ. of California Berkeley
94: *> \author Univ. of Colorado Denver
95: *> \author NAG Ltd.
96: *
97: *> \ingroup OTHERauxiliary
98: *
99: *> \par Further Details:
100: * =====================
101: *>
102: *> \verbatim
103: *>
104: *> RT1 is accurate to a few ulps barring over/underflow.
105: *>
106: *> RT2 may be inaccurate if there is massive cancellation in the
107: *> determinant A*C-B*B; higher precision or correctly rounded or
108: *> correctly truncated arithmetic would be needed to compute RT2
109: *> accurately in all cases.
110: *>
111: *> CS1 and SN1 are accurate to a few ulps barring over/underflow.
112: *>
113: *> Overflow is possible only if RT1 is within a factor of 5 of overflow.
114: *> Underflow is harmless if the input data is 0 or exceeds
115: *> underflow_threshold / macheps.
116: *> \endverbatim
117: *>
118: * =====================================================================
119: SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
120: *
121: * -- LAPACK auxiliary routine --
122: * -- LAPACK is a software package provided by Univ. of Tennessee, --
123: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124: *
125: * .. Scalar Arguments ..
126: DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1
127: * ..
128: *
129: * =====================================================================
130: *
131: * .. Parameters ..
132: DOUBLE PRECISION ONE
133: PARAMETER ( ONE = 1.0D0 )
134: DOUBLE PRECISION TWO
135: PARAMETER ( TWO = 2.0D0 )
136: DOUBLE PRECISION ZERO
137: PARAMETER ( ZERO = 0.0D0 )
138: DOUBLE PRECISION HALF
139: PARAMETER ( HALF = 0.5D0 )
140: * ..
141: * .. Local Scalars ..
142: INTEGER SGN1, SGN2
143: DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
144: $ TB, TN
145: * ..
146: * .. Intrinsic Functions ..
147: INTRINSIC ABS, SQRT
148: * ..
149: * .. Executable Statements ..
150: *
151: * Compute the eigenvalues
152: *
153: SM = A + C
154: DF = A - C
155: ADF = ABS( DF )
156: TB = B + B
157: AB = ABS( TB )
158: IF( ABS( A ).GT.ABS( C ) ) THEN
159: ACMX = A
160: ACMN = C
161: ELSE
162: ACMX = C
163: ACMN = A
164: END IF
165: IF( ADF.GT.AB ) THEN
166: RT = ADF*SQRT( ONE+( AB / ADF )**2 )
167: ELSE IF( ADF.LT.AB ) THEN
168: RT = AB*SQRT( ONE+( ADF / AB )**2 )
169: ELSE
170: *
171: * Includes case AB=ADF=0
172: *
173: RT = AB*SQRT( TWO )
174: END IF
175: IF( SM.LT.ZERO ) THEN
176: RT1 = HALF*( SM-RT )
177: SGN1 = -1
178: *
179: * Order of execution important.
180: * To get fully accurate smaller eigenvalue,
181: * next line needs to be executed in higher precision.
182: *
183: RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
184: ELSE IF( SM.GT.ZERO ) THEN
185: RT1 = HALF*( SM+RT )
186: SGN1 = 1
187: *
188: * Order of execution important.
189: * To get fully accurate smaller eigenvalue,
190: * next line needs to be executed in higher precision.
191: *
192: RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
193: ELSE
194: *
195: * Includes case RT1 = RT2 = 0
196: *
197: RT1 = HALF*RT
198: RT2 = -HALF*RT
199: SGN1 = 1
200: END IF
201: *
202: * Compute the eigenvector
203: *
204: IF( DF.GE.ZERO ) THEN
205: CS = DF + RT
206: SGN2 = 1
207: ELSE
208: CS = DF - RT
209: SGN2 = -1
210: END IF
211: ACS = ABS( CS )
212: IF( ACS.GT.AB ) THEN
213: CT = -TB / CS
214: SN1 = ONE / SQRT( ONE+CT*CT )
215: CS1 = CT*SN1
216: ELSE
217: IF( AB.EQ.ZERO ) THEN
218: CS1 = ONE
219: SN1 = ZERO
220: ELSE
221: TN = -CS / TB
222: CS1 = ONE / SQRT( ONE+TN*TN )
223: SN1 = TN*CS1
224: END IF
225: END IF
226: IF( SGN1.EQ.SGN2 ) THEN
227: TN = CS1
228: CS1 = -SN1
229: SN1 = TN
230: END IF
231: RETURN
232: *
233: * End of DLAEV2
234: *
235: END
CVSweb interface <joel.bertrand@systella.fr>