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
171