001:       SUBROUTINE SLAE2( A, B, C, RT1, RT2 )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       REAL               A, B, C, RT1, RT2
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  SLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
015: *     [  A   B  ]
016: *     [  B   C  ].
017: *  On return, RT1 is the eigenvalue of larger absolute value, and RT2
018: *  is the eigenvalue of smaller absolute value.
019: *
020: *  Arguments
021: *  =========
022: *
023: *  A       (input) REAL
024: *          The (1,1) element of the 2-by-2 matrix.
025: *
026: *  B       (input) REAL
027: *          The (1,2) and (2,1) elements of the 2-by-2 matrix.
028: *
029: *  C       (input) REAL
030: *          The (2,2) element of the 2-by-2 matrix.
031: *
032: *  RT1     (output) REAL
033: *          The eigenvalue of larger absolute value.
034: *
035: *  RT2     (output) REAL
036: *          The eigenvalue of smaller absolute value.
037: *
038: *  Further Details
039: *  ===============
040: *
041: *  RT1 is accurate to a few ulps barring over/underflow.
042: *
043: *  RT2 may be inaccurate if there is massive cancellation in the
044: *  determinant A*C-B*B; higher precision or correctly rounded or
045: *  correctly truncated arithmetic would be needed to compute RT2
046: *  accurately in all cases.
047: *
048: *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
049: *  Underflow is harmless if the input data is 0 or exceeds
050: *     underflow_threshold / macheps.
051: *
052: * =====================================================================
053: *
054: *     .. Parameters ..
055:       REAL               ONE
056:       PARAMETER          ( ONE = 1.0E0 )
057:       REAL               TWO
058:       PARAMETER          ( TWO = 2.0E0 )
059:       REAL               ZERO
060:       PARAMETER          ( ZERO = 0.0E0 )
061:       REAL               HALF
062:       PARAMETER          ( HALF = 0.5E0 )
063: *     ..
064: *     .. Local Scalars ..
065:       REAL               AB, ACMN, ACMX, ADF, DF, RT, SM, TB
066: *     ..
067: *     .. Intrinsic Functions ..
068:       INTRINSIC          ABS, SQRT
069: *     ..
070: *     .. Executable Statements ..
071: *
072: *     Compute the eigenvalues
073: *
074:       SM = A + C
075:       DF = A - C
076:       ADF = ABS( DF )
077:       TB = B + B
078:       AB = ABS( TB )
079:       IF( ABS( A ).GT.ABS( C ) ) THEN
080:          ACMX = A
081:          ACMN = C
082:       ELSE
083:          ACMX = C
084:          ACMN = A
085:       END IF
086:       IF( ADF.GT.AB ) THEN
087:          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
088:       ELSE IF( ADF.LT.AB ) THEN
089:          RT = AB*SQRT( ONE+( ADF / AB )**2 )
090:       ELSE
091: *
092: *        Includes case AB=ADF=0
093: *
094:          RT = AB*SQRT( TWO )
095:       END IF
096:       IF( SM.LT.ZERO ) THEN
097:          RT1 = HALF*( SM-RT )
098: *
099: *        Order of execution important.
100: *        To get fully accurate smaller eigenvalue,
101: *        next line needs to be executed in higher precision.
102: *
103:          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
104:       ELSE IF( SM.GT.ZERO ) THEN
105:          RT1 = HALF*( SM+RT )
106: *
107: *        Order of execution important.
108: *        To get fully accurate smaller eigenvalue,
109: *        next line needs to be executed in higher precision.
110: *
111:          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
112:       ELSE
113: *
114: *        Includes case RT1 = RT2 = 0
115: *
116:          RT1 = HALF*RT
117:          RT2 = -HALF*RT
118:       END IF
119:       RETURN
120: *
121: *     End of SLAE2
122: *
123:       END
124: