vendor/gonum.org/v1/gonum/lapack/internal/testdata/dsterftest/dlae2.f

   1 *> \brief \b DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAE2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlae2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlae2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlae2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       DOUBLE PRECISION   A, B, C, RT1, RT2
  25 *       ..
  26 *
  27 *
  28 *> \par Purpose:
  29 *  =============
  30 *>
  31 *> \verbatim
  32 *>
  33 *> DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
  34 *>    [  A   B  ]
  35 *>    [  B   C  ].
  36 *> On return, RT1 is the eigenvalue of larger absolute value, and RT2
  37 *> is the eigenvalue of smaller absolute value.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] A
  44 *> \verbatim
  45 *>          A is DOUBLE PRECISION
  46 *>          The (1,1) element of the 2-by-2 matrix.
  47 *> \endverbatim
  48 *>
  49 *> \param[in] B
  50 *> \verbatim
  51 *>          B is DOUBLE PRECISION
  52 *>          The (1,2) and (2,1) elements of the 2-by-2 matrix.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] C
  56 *> \verbatim
  57 *>          C is DOUBLE PRECISION
  58 *>          The (2,2) element of the 2-by-2 matrix.
  59 *> \endverbatim
  60 *>
  61 *> \param[out] RT1
  62 *> \verbatim
  63 *>          RT1 is DOUBLE PRECISION
  64 *>          The eigenvalue of larger absolute value.
  65 *> \endverbatim
  66 *>
  67 *> \param[out] RT2
  68 *> \verbatim
  69 *>          RT2 is DOUBLE PRECISION
  70 *>          The eigenvalue of smaller absolute value.
  71 *> \endverbatim
  72 *
  73 *  Authors:
  74 *  ========
  75 *
  76 *> \author Univ. of Tennessee
  77 *> \author Univ. of California Berkeley
  78 *> \author Univ. of Colorado Denver
  79 *> \author NAG Ltd.
  80 *
  81 *> \date September 2012
  82 *
  83 *> \ingroup auxOTHERauxiliary
  84 *
  85 *> \par Further Details:
  86 *  =====================
  87 *>
  88 *> \verbatim
  89 *>
  90 *>  RT1 is accurate to a few ulps barring over/underflow.
  91 *>
  92 *>  RT2 may be inaccurate if there is massive cancellation in the
  93 *>  determinant A*C-B*B; higher precision or correctly rounded or
  94 *>  correctly truncated arithmetic would be needed to compute RT2
  95 *>  accurately in all cases.
  96 *>
  97 *>  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  98 *>  Underflow is harmless if the input data is 0 or exceeds
  99 *>     underflow_threshold / macheps.
 100 *> \endverbatim
 101 *>
 102 *  =====================================================================
 103       SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
 104 *
 105 *  -- LAPACK auxiliary routine (version 3.4.2) --
 106 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 107 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 108 *     September 2012
 109 *
 110 *     .. Scalar Arguments ..
 111       DOUBLE PRECISION   A, B, C, RT1, RT2
 112 *     ..
 113 *
 114 * =====================================================================
 115 *
 116 *     .. Parameters ..
 117       DOUBLE PRECISION   ONE
 118       PARAMETER          ( ONE = 1.0D0 )
 119       DOUBLE PRECISION   TWO
 120       PARAMETER          ( TWO = 2.0D0 )
 121       DOUBLE PRECISION   ZERO
 122       PARAMETER          ( ZERO = 0.0D0 )
 123       DOUBLE PRECISION   HALF
 124       PARAMETER          ( HALF = 0.5D0 )
 125 *     ..
 126 *     .. Local Scalars ..
 127       DOUBLE PRECISION   AB, ACMN, ACMX, ADF, DF, RT, SM, TB
 128 *     ..
 129 *     .. Intrinsic Functions ..
 130       INTRINSIC          ABS, SQRT
 131 *     ..
 132 *     .. Executable Statements ..
 133 *
 134 *     Compute the eigenvalues
 135 *
 136       SM = A + C
 137       DF = A - C
 138       ADF = ABS( DF )
 139       TB = B + B
 140       AB = ABS( TB )
 141       IF( ABS( A ).GT.ABS( C ) ) THEN
 142          ACMX = A
 143          ACMN = C
 144       ELSE
 145          ACMX = C
 146          ACMN = A
 147       END IF
 148       IF( ADF.GT.AB ) THEN
 149          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
 150       ELSE IF( ADF.LT.AB ) THEN
 151          RT = AB*SQRT( ONE+( ADF / AB )**2 )
 152       ELSE
 153 *
 154 *        Includes case AB=ADF=0
 155 *
 156          RT = AB*SQRT( TWO )
 157       END IF
 158       IF( SM.LT.ZERO ) THEN
 159          RT1 = HALF*( SM-RT )
 160 *
 161 *        Order of execution important.
 162 *        To get fully accurate smaller eigenvalue,
 163 *        next line needs to be executed in higher precision.
 164 *
 165          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 166       ELSE IF( SM.GT.ZERO ) THEN
 167          RT1 = HALF*( SM+RT )
 168 *
 169 *        Order of execution important.
 170 *        To get fully accurate smaller eigenvalue,
 171 *        next line needs to be executed in higher precision.
 172 *
 173          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 174       ELSE
 175 *
 176 *        Includes case RT1 = RT2 = 0
 177 *
 178          RT1 = HALF*RT
 179          RT2 = -HALF*RT
 180       END IF
 181       RETURN
 182 *
 183 *     End of DLAE2
 184 *
 185       END