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*> \date September 2012
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*  =====================================================================
121      SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
122*
123*  -- LAPACK auxiliary routine (version 3.4.2) --
124*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
125*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*     September 2012
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
239c $Id$
240