Annotation of rpl/lapack/lapack/dlaev2.f, revision 1.14
1.11 bertrand 1: *> \brief \b DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
1.8 bertrand 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: *
1.11 bertrand 97: *> \date September 2012
1.8 bertrand 98: *
99: *> \ingroup auxOTHERauxiliary
100: *
101: *> \par Further Details:
102: * =====================
103: *>
104: *> \verbatim
105: *>
106: *> RT1 is accurate to a few ulps barring over/underflow.
107: *>
108: *> RT2 may be inaccurate if there is massive cancellation in the
109: *> determinant A*C-B*B; higher precision or correctly rounded or
110: *> correctly truncated arithmetic would be needed to compute RT2
111: *> accurately in all cases.
112: *>
113: *> CS1 and SN1 are accurate to a few ulps barring over/underflow.
114: *>
115: *> Overflow is possible only if RT1 is within a factor of 5 of overflow.
116: *> Underflow is harmless if the input data is 0 or exceeds
117: *> underflow_threshold / macheps.
118: *> \endverbatim
119: *>
120: * =====================================================================
1.1 bertrand 121: SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
122: *
1.11 bertrand 123: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 124: * -- LAPACK is a software package provided by Univ. of Tennessee, --
125: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 126: * September 2012
1.1 bertrand 127: *
128: * .. Scalar Arguments ..
129: DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1
130: * ..
131: *
132: * =====================================================================
133: *
134: * .. Parameters ..
135: DOUBLE PRECISION ONE
136: PARAMETER ( ONE = 1.0D0 )
137: DOUBLE PRECISION TWO
138: PARAMETER ( TWO = 2.0D0 )
139: DOUBLE PRECISION ZERO
140: PARAMETER ( ZERO = 0.0D0 )
141: DOUBLE PRECISION HALF
142: PARAMETER ( HALF = 0.5D0 )
143: * ..
144: * .. Local Scalars ..
145: INTEGER SGN1, SGN2
146: DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
147: $ TB, TN
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC ABS, SQRT
151: * ..
152: * .. Executable Statements ..
153: *
154: * Compute the eigenvalues
155: *
156: SM = A + C
157: DF = A - C
158: ADF = ABS( DF )
159: TB = B + B
160: AB = ABS( TB )
161: IF( ABS( A ).GT.ABS( C ) ) THEN
162: ACMX = A
163: ACMN = C
164: ELSE
165: ACMX = C
166: ACMN = A
167: END IF
168: IF( ADF.GT.AB ) THEN
169: RT = ADF*SQRT( ONE+( AB / ADF )**2 )
170: ELSE IF( ADF.LT.AB ) THEN
171: RT = AB*SQRT( ONE+( ADF / AB )**2 )
172: ELSE
173: *
174: * Includes case AB=ADF=0
175: *
176: RT = AB*SQRT( TWO )
177: END IF
178: IF( SM.LT.ZERO ) THEN
179: RT1 = HALF*( SM-RT )
180: SGN1 = -1
181: *
182: * Order of execution important.
183: * To get fully accurate smaller eigenvalue,
184: * next line needs to be executed in higher precision.
185: *
186: RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
187: ELSE IF( SM.GT.ZERO ) THEN
188: RT1 = HALF*( SM+RT )
189: SGN1 = 1
190: *
191: * Order of execution important.
192: * To get fully accurate smaller eigenvalue,
193: * next line needs to be executed in higher precision.
194: *
195: RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
196: ELSE
197: *
198: * Includes case RT1 = RT2 = 0
199: *
200: RT1 = HALF*RT
201: RT2 = -HALF*RT
202: SGN1 = 1
203: END IF
204: *
205: * Compute the eigenvector
206: *
207: IF( DF.GE.ZERO ) THEN
208: CS = DF + RT
209: SGN2 = 1
210: ELSE
211: CS = DF - RT
212: SGN2 = -1
213: END IF
214: ACS = ABS( CS )
215: IF( ACS.GT.AB ) THEN
216: CT = -TB / CS
217: SN1 = ONE / SQRT( ONE+CT*CT )
218: CS1 = CT*SN1
219: ELSE
220: IF( AB.EQ.ZERO ) THEN
221: CS1 = ONE
222: SN1 = ZERO
223: ELSE
224: TN = -CS / TB
225: CS1 = ONE / SQRT( ONE+TN*TN )
226: SN1 = TN*CS1
227: END IF
228: END IF
229: IF( SGN1.EQ.SGN2 ) THEN
230: TN = CS1
231: CS1 = -SN1
232: SN1 = TN
233: END IF
234: RETURN
235: *
236: * End of DLAEV2
237: *
238: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaev2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack