Return to dlae2.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dlae2.f, revision 1.5

1.1       bertrand    1:       SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
                      2: *
                      3: *  -- LAPACK auxiliary routine (version 3.2) --
                      4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      6: *     November 2006
                      7: *
                      8: *     .. Scalar Arguments ..
                      9:       DOUBLE PRECISION   A, B, C, RT1, RT2
                     10: *     ..
                     11: *
                     12: *  Purpose
                     13: *  =======
                     14: *
                     15: *  DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
                     16: *     [  A   B  ]
                     17: *     [  B   C  ].
                     18: *  On return, RT1 is the eigenvalue of larger absolute value, and RT2
                     19: *  is the eigenvalue of smaller absolute value.
                     20: *
                     21: *  Arguments
                     22: *  =========
                     23: *
                     24: *  A       (input) DOUBLE PRECISION
                     25: *          The (1,1) element of the 2-by-2 matrix.
                     26: *
                     27: *  B       (input) DOUBLE PRECISION
                     28: *          The (1,2) and (2,1) elements of the 2-by-2 matrix.
                     29: *
                     30: *  C       (input) DOUBLE PRECISION
                     31: *          The (2,2) element of the 2-by-2 matrix.
                     32: *
                     33: *  RT1     (output) DOUBLE PRECISION
                     34: *          The eigenvalue of larger absolute value.
                     35: *
                     36: *  RT2     (output) DOUBLE PRECISION
                     37: *          The eigenvalue of smaller absolute value.
                     38: *
                     39: *  Further Details
                     40: *  ===============
                     41: *
                     42: *  RT1 is accurate to a few ulps barring over/underflow.
                     43: *
                     44: *  RT2 may be inaccurate if there is massive cancellation in the
                     45: *  determinant A*C-B*B; higher precision or correctly rounded or
                     46: *  correctly truncated arithmetic would be needed to compute RT2
                     47: *  accurately in all cases.
                     48: *
                     49: *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
                     50: *  Underflow is harmless if the input data is 0 or exceeds
                     51: *     underflow_threshold / macheps.
                     52: *
                     53: * =====================================================================
                     54: *
                     55: *     .. Parameters ..
                     56:       DOUBLE PRECISION   ONE
                     57:       PARAMETER          ( ONE = 1.0D0 )
                     58:       DOUBLE PRECISION   TWO
                     59:       PARAMETER          ( TWO = 2.0D0 )
                     60:       DOUBLE PRECISION   ZERO
                     61:       PARAMETER          ( ZERO = 0.0D0 )
                     62:       DOUBLE PRECISION   HALF
                     63:       PARAMETER          ( HALF = 0.5D0 )
                     64: *     ..
                     65: *     .. Local Scalars ..
                     66:       DOUBLE PRECISION   AB, ACMN, ACMX, ADF, DF, RT, SM, TB
                     67: *     ..
                     68: *     .. Intrinsic Functions ..
                     69:       INTRINSIC          ABS, SQRT
                     70: *     ..
                     71: *     .. Executable Statements ..
                     72: *
                     73: *     Compute the eigenvalues
                     74: *
                     75:       SM = A + C
                     76:       DF = A - C
                     77:       ADF = ABS( DF )
                     78:       TB = B + B
                     79:       AB = ABS( TB )
                     80:       IF( ABS( A ).GT.ABS( C ) ) THEN
                     81:          ACMX = A
                     82:          ACMN = C
                     83:       ELSE
                     84:          ACMX = C
                     85:          ACMN = A
                     86:       END IF
                     87:       IF( ADF.GT.AB ) THEN
                     88:          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
                     89:       ELSE IF( ADF.LT.AB ) THEN
                     90:          RT = AB*SQRT( ONE+( ADF / AB )**2 )
                     91:       ELSE
                     92: *
                     93: *        Includes case AB=ADF=0
                     94: *
                     95:          RT = AB*SQRT( TWO )
                     96:       END IF
                     97:       IF( SM.LT.ZERO ) THEN
                     98:          RT1 = HALF*( SM-RT )
                     99: *
                    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 DLAE2
                    123: *
                    124:       END