var/external/lapack/dlaev2.inc

   1       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
   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, CS1, RT1, RT2, SN1
   9 !     ..
  10 !
  11 !  Purpose
  12 !  =======
  13 !
  14 !  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
  15 !     [  A   B  ]
  16 !     [  B   C  ].
  17 !  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
  18 !  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
  19 !  eigenvector for RT1, giving the decomposition
  20 !
  21 !     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
  22 !     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
  23 !
  24 !  Arguments
  25 !  =========
  26 !
  27 !  A       (input) DOUBLE PRECISION
  28 !          The (1,1) element of the 2-by-2 matrix.
  29 !
  30 !  B       (input) DOUBLE PRECISION
  31 !          The (1,2) element and the conjugate of the (2,1) element of
  32 !          the 2-by-2 matrix.
  33 !
  34 !  C       (input) DOUBLE PRECISION
  35 !          The (2,2) element of the 2-by-2 matrix.
  36 !
  37 !  RT1     (output) DOUBLE PRECISION
  38 !          The eigenvalue of larger absolute value.
  39 !
  40 !  RT2     (output) DOUBLE PRECISION
  41 !          The eigenvalue of smaller absolute value.
  42 !
  43 !  CS1     (output) DOUBLE PRECISION
  44 !  SN1     (output) DOUBLE PRECISION
  45 !          The vector (CS1, SN1) is a unit right eigenvector for RT1.
  46 !
  47 !  Further Details
  48 !  ===============
  49 !
  50 !  RT1 is accurate to a few ulps barring over/underflow.
  51 !
  52 !  RT2 may be inaccurate if there is massive cancellation in the
  53 !  determinant A*C-B*B; higher precision or correctly rounded or
  54 !  correctly truncated arithmetic would be needed to compute RT2
  55 !  accurately in all cases.
  56 !
  57 !  CS1 and SN1 are accurate to a few ulps barring over/underflow.
  58 !
  59 !  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  60 !  Underflow is harmless if the input data is 0 or exceeds
  61 !     underflow_threshold / macheps.
  62 !
  63 ! =====================================================================
  64 !
  65 !     .. Parameters ..
  66       DOUBLE PRECISION   ONE
  67       PARAMETER          ( ONE = 1.0D0 )
  68       DOUBLE PRECISION   TWO
  69       PARAMETER          ( TWO = 2.0D0 )
  70       DOUBLE PRECISION   ZERO
  71       PARAMETER          ( ZERO = 0.0D0 )
  72       DOUBLE PRECISION   HALF
  73       PARAMETER          ( HALF = 0.5D0 )
  74 !     ..
  75 !     .. Local Scalars ..
  76       INTEGER            SGN1, SGN2
  77       DOUBLE PRECISION   AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, &
  78                          TB, TN
  79 !     ..
  80 !     .. Intrinsic Functions ..
  81       INTRINSIC          ABS, SQRT
  82 !     ..
  83 !     .. Executable Statements ..
  84 !
  85 !     Compute the eigenvalues
  86 !
  87       SM = A + C
  88       DF = A - C
  89       ADF = ABS( DF )
  90       TB = B + B
  91       AB = ABS( TB )
  92       IF( ABS( A ).GT.ABS( C ) ) THEN
  93          ACMX = A
  94          ACMN = C
  95       ELSE
  96          ACMX = C
  97          ACMN = A
  98       END IF
  99       IF( ADF.GT.AB ) THEN
 100          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
 101       ELSE IF( ADF.LT.AB ) THEN
 102          RT = AB*SQRT( ONE+( ADF / AB )**2 )
 103       ELSE
 104 !
 105 !        Includes case AB=ADF=0
 106 !
 107          RT = AB*SQRT( TWO )
 108       END IF
 109       IF( SM.LT.ZERO ) THEN
 110          RT1 = HALF*( SM-RT )
 111          SGN1 = -1
 112 !
 113 !        Order of execution important.
 114 !        To get fully accurate smaller eigenvalue,
 115 !        next line needs to be executed in higher precision.
 116 !
 117          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 118       ELSE IF( SM.GT.ZERO ) THEN
 119          RT1 = HALF*( SM+RT )
 120          SGN1 = 1
 121 !
 122 !        Order of execution important.
 123 !        To get fully accurate smaller eigenvalue,
 124 !        next line needs to be executed in higher precision.
 125 !
 126          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 127       ELSE
 128 !
 129 !        Includes case RT1 = RT2 = 0
 130 !
 131          RT1 = HALF*RT
 132          RT2 = -HALF*RT
 133          SGN1 = 1
 134       END IF
 135 !
 136 !     Compute the eigenvector
 137 !
 138       IF( DF.GE.ZERO ) THEN
 139          CS = DF + RT
 140          SGN2 = 1
 141       ELSE
 142          CS = DF - RT
 143          SGN2 = -1
 144       END IF
 145       ACS = ABS( CS )
 146       IF( ACS.GT.AB ) THEN
 147          CT = -TB / CS
 148          SN1 = ONE / SQRT( ONE+CT*CT )
 149          CS1 = CT*SN1
 150       ELSE
 151          IF( AB.EQ.ZERO ) THEN
 152             CS1 = ONE
 153             SN1 = ZERO
 154          ELSE
 155             TN = -CS / TB
 156             CS1 = ONE / SQRT( ONE+TN*TN )
 157             SN1 = TN*CS1
 158          END IF
 159       END IF
 160       IF( SGN1.EQ.SGN2 ) THEN
 161          TN = CS1
 162          CS1 = -SN1
 163          SN1 = TN
 164       END IF
 165       RETURN
 166 !
 167 !     End of DLAEV2
 168 !
 169       END SUBROUTINE DLAEV2