d4/d37/slae2_8f_source.html

00001       SUBROUTINE SLAE2( A, B, C, RT1, RT2 )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       REAL               A, B, C, RT1, RT2
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  SLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
00016 *     [  A   B  ]
00017 *     [  B   C  ].
00018 *  On return, RT1 is the eigenvalue of larger absolute value, and RT2
00019 *  is the eigenvalue of smaller absolute value.
00020 *
00021 *  Arguments
00022 *  =========
00023 *
00024 *  A       (input) REAL
00025 *          The (1,1) element of the 2-by-2 matrix.
00026 *
00027 *  B       (input) REAL
00028 *          The (1,2) and (2,1) elements of the 2-by-2 matrix.
00029 *
00030 *  C       (input) REAL
00031 *          The (2,2) element of the 2-by-2 matrix.
00032 *
00033 *  RT1     (output) REAL
00034 *          The eigenvalue of larger absolute value.
00035 *
00036 *  RT2     (output) REAL
00037 *          The eigenvalue of smaller absolute value.
00038 *
00039 *  Further Details
00040 *  ===============
00041 *
00042 *  RT1 is accurate to a few ulps barring over/underflow.
00043 *
00044 *  RT2 may be inaccurate if there is massive cancellation in the
00045 *  determinant A*C-B*B; higher precision or correctly rounded or
00046 *  correctly truncated arithmetic would be needed to compute RT2
00047 *  accurately in all cases.
00048 *
00049 *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
00050 *  Underflow is harmless if the input data is 0 or exceeds
00051 *     underflow_threshold / macheps.
00052 *
00053 * =====================================================================
00054 *
00055 *     .. Parameters ..
00056       REAL               ONE
00057       PARAMETER          ( ONE = 1.0E0 )
00058       REAL               TWO
00059       PARAMETER          ( TWO = 2.0E0 )
00060       REAL               ZERO
00061       PARAMETER          ( ZERO = 0.0E0 )
00062       REAL               HALF
00063       PARAMETER          ( HALF = 0.5E0 )
00064 *     ..
00065 *     .. Local Scalars ..
00066       REAL               AB, ACMN, ACMX, ADF, DF, RT, SM, TB
00067 *     ..
00068 *     .. Intrinsic Functions ..
00069       INTRINSIC          ABS, SQRT
00070 *     ..
00071 *     .. Executable Statements ..
00072 *
00073 *     Compute the eigenvalues
00074 *
00075       SM = A + C
00076       DF = A - C
00077       ADF = ABS( DF )
00078       TB = B + B
00079       AB = ABS( TB )
00080       IF( ABS( A ).GT.ABS( C ) ) THEN
00081          ACMX = A
00082          ACMN = C
00083       ELSE
00084          ACMX = C
00085          ACMN = A
00086       END IF
00087       IF( ADF.GT.AB ) THEN
00088          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
00089       ELSE IF( ADF.LT.AB ) THEN
00090          RT = AB*SQRT( ONE+( ADF / AB )**2 )
00091       ELSE
00092 *
00093 *        Includes case AB=ADF=0
00094 *
00095          RT = AB*SQRT( TWO )
00096       END IF
00097       IF( SM.LT.ZERO ) THEN
00098          RT1 = HALF*( SM-RT )
00099 *
00100 *        Order of execution important.
00101 *        To get fully accurate smaller eigenvalue,
00102 *        next line needs to be executed in higher precision.
00103 *
00104          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
00105       ELSE IF( SM.GT.ZERO ) THEN
00106          RT1 = HALF*( SM+RT )
00107 *
00108 *        Order of execution important.
00109 *        To get fully accurate smaller eigenvalue,
00110 *        next line needs to be executed in higher precision.
00111 *
00112          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
00113       ELSE
00114 *
00115 *        Includes case RT1 = RT2 = 0
00116 *
00117          RT1 = HALF*RT
00118          RT2 = -HALF*RT
00119       END IF
00120       RETURN
00121 *
00122 *     End of SLAE2
00123 *
00124       END