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