var/external/lapack/dlae2.inc

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