001:       SUBROUTINE SLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       REAL               A, B, C, CS1, RT1, RT2, SN1
010: *     ..
011: *
012: *  Purpose
013: *  =======
014: *
015: *  SLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
016: *     [  A   B  ]
017: *     [  B   C  ].
018: *  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
019: *  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
020: *  eigenvector for RT1, giving the decomposition
021: *
022: *     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
023: *     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
024: *
025: *  Arguments
026: *  =========
027: *
028: *  A       (input) REAL
029: *          The (1,1) element of the 2-by-2 matrix.
030: *
031: *  B       (input) REAL
032: *          The (1,2) element and the conjugate of the (2,1) element of
033: *          the 2-by-2 matrix.
034: *
035: *  C       (input) REAL
036: *          The (2,2) element of the 2-by-2 matrix.
037: *
038: *  RT1     (output) REAL
039: *          The eigenvalue of larger absolute value.
040: *
041: *  RT2     (output) REAL
042: *          The eigenvalue of smaller absolute value.
043: *
044: *  CS1     (output) REAL
045: *  SN1     (output) REAL
046: *          The vector (CS1, SN1) is a unit right eigenvector for RT1.
047: *
048: *  Further Details
049: *  ===============
050: *
051: *  RT1 is accurate to a few ulps barring over/underflow.
052: *
053: *  RT2 may be inaccurate if there is massive cancellation in the
054: *  determinant A*C-B*B; higher precision or correctly rounded or
055: *  correctly truncated arithmetic would be needed to compute RT2
056: *  accurately in all cases.
057: *
058: *  CS1 and SN1 are accurate to a few ulps barring over/underflow.
059: *
060: *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
061: *  Underflow is harmless if the input data is 0 or exceeds
062: *     underflow_threshold / macheps.
063: *
064: * =====================================================================
065: *
066: *     .. Parameters ..
067:       REAL               ONE
068:       PARAMETER          ( ONE = 1.0E0 )
069:       REAL               TWO
070:       PARAMETER          ( TWO = 2.0E0 )
071:       REAL               ZERO
072:       PARAMETER          ( ZERO = 0.0E0 )
073:       REAL               HALF
074:       PARAMETER          ( HALF = 0.5E0 )
075: *     ..
076: *     .. Local Scalars ..
077:       INTEGER            SGN1, SGN2
078:       REAL               AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
079:      $                   TB, TN
080: *     ..
081: *     .. Intrinsic Functions ..
082:       INTRINSIC          ABS, SQRT
083: *     ..
084: *     .. Executable Statements ..
085: *
086: *     Compute the eigenvalues
087: *
088:       SM = A + C
089:       DF = A - C
090:       ADF = ABS( DF )
091:       TB = B + B
092:       AB = ABS( TB )
093:       IF( ABS( A ).GT.ABS( C ) ) THEN
094:          ACMX = A
095:          ACMN = C
096:       ELSE
097:          ACMX = C
098:          ACMN = A
099:       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 SLAEV2
169: *
170:       END
171: