Return to dlae2.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    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