KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlaev2.f
Go to the documentation of this file.
1  SUBROUTINE dlaev2( A, B, C, RT1, RT2, CS1, SN1 )
2 *
3 * -- LAPACK auxiliary routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * October 31, 1992
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
subroutine dlaev2(A, B, C, RT1, RT2, CS1, SN1)
Definition: dlaev2.f:2