SRC/slaev2.f

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