Annotation of rpl/lapack/lapack/dlaev2.f, revision 1.5
1.1 bertrand 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
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