SRC/claesy.f

   1 *> \brief \b CLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex 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 CLAESY + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claesy.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claesy.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claesy.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
  25 *       ..
  26 *
  27 *
  28 *> \par Purpose:
  29 *  =============
  30 *>
  31 *> \verbatim
  32 *>
  33 *> CLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
  34 *>    ( ( A, B );( B, C ) )
  35 *> provided the norm of the matrix of eigenvectors is larger than
  36 *> some threshold value.
  37 *>
  38 *> RT1 is the eigenvalue of larger absolute value, and RT2 of
  39 *> smaller absolute value.  If the eigenvectors are computed, then
  40 *> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
  41 *>
  42 *> [  CS1     SN1   ] . [ A  B ] . [ CS1    -SN1   ] = [ RT1  0  ]
  43 *> [ -SN1     CS1   ]   [ B  C ]   [ SN1     CS1   ]   [  0  RT2 ]
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] A
  50 *> \verbatim
  51 *>          A is COMPLEX
  52 *>          The ( 1, 1 ) element of input matrix.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] B
  56 *> \verbatim
  57 *>          B is COMPLEX
  58 *>          The ( 1, 2 ) element of input matrix.  The ( 2, 1 ) element
  59 *>          is also given by B, since the 2-by-2 matrix is symmetric.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] C
  63 *> \verbatim
  64 *>          C is COMPLEX
  65 *>          The ( 2, 2 ) element of input matrix.
  66 *> \endverbatim
  67 *>
  68 *> \param[out] RT1
  69 *> \verbatim
  70 *>          RT1 is COMPLEX
  71 *>          The eigenvalue of larger modulus.
  72 *> \endverbatim
  73 *>
  74 *> \param[out] RT2
  75 *> \verbatim
  76 *>          RT2 is COMPLEX
  77 *>          The eigenvalue of smaller modulus.
  78 *> \endverbatim
  79 *>
  80 *> \param[out] EVSCAL
  81 *> \verbatim
  82 *>          EVSCAL is COMPLEX
  83 *>          The complex value by which the eigenvector matrix was scaled
  84 *>          to make it orthonormal.  If EVSCAL is zero, the eigenvectors
  85 *>          were not computed.  This means one of two things:  the 2-by-2
  86 *>          matrix could not be diagonalized, or the norm of the matrix
  87 *>          of eigenvectors before scaling was larger than the threshold
  88 *>          value THRESH (set below).
  89 *> \endverbatim
  90 *>
  91 *> \param[out] CS1
  92 *> \verbatim
  93 *>          CS1 is COMPLEX
  94 *> \endverbatim
  95 *>
  96 *> \param[out] SN1
  97 *> \verbatim
  98 *>          SN1 is COMPLEX
  99 *>          If EVSCAL .NE. 0,  ( CS1, SN1 ) is the unit right eigenvector
 100 *>          for RT1.
 101 *> \endverbatim
 102 *
 103 *  Authors:
 104 *  ========
 105 *
 106 *> \author Univ. of Tennessee
 107 *> \author Univ. of California Berkeley
 108 *> \author Univ. of Colorado Denver
 109 *> \author NAG Ltd.
 110 *
 111 *> \date September 2012
 112 *
 113 *> \ingroup complexSYauxiliary
 114 *
 115 *  =====================================================================
 116       SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
 117 *
 118 *  -- LAPACK auxiliary routine (version 3.4.2) --
 119 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 120 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 121 *     September 2012
 122 *
 123 *     .. Scalar Arguments ..
 124       COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
 125 *     ..
 126 *
 127 * =====================================================================
 128 *
 129 *     .. Parameters ..
 130       REAL               ZERO
 131       PARAMETER          ( ZERO = 0.0E0 )
 132       REAL               ONE
 133       PARAMETER          ( ONE = 1.0E0 )
 134       COMPLEX            CONE
 135       PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
 136       REAL               HALF
 137       PARAMETER          ( HALF = 0.5E0 )
 138       REAL               THRESH
 139       PARAMETER          ( THRESH = 0.1E0 )
 140 *     ..
 141 *     .. Local Scalars ..
 142       REAL               BABS, EVNORM, TABS, Z
 143       COMPLEX            S, T, TMP
 144 *     ..
 145 *     .. Intrinsic Functions ..
 146       INTRINSIC          ABS, MAX, SQRT
 147 *     ..
 148 *     .. Executable Statements ..
 149 *
 150 *
 151 *     Special case:  The matrix is actually diagonal.
 152 *     To avoid divide by zero later, we treat this case separately.
 153 *
 154       IF( ABS( B ).EQ.ZERO ) THEN
 155          RT1 = A
 156          RT2 = C
 157          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
 158             TMP = RT1
 159             RT1 = RT2
 160             RT2 = TMP
 161             CS1 = ZERO
 162             SN1 = ONE
 163          ELSE
 164             CS1 = ONE
 165             SN1 = ZERO
 166          END IF
 167       ELSE
 168 *
 169 *        Compute the eigenvalues and eigenvectors.
 170 *        The characteristic equation is
 171 *           lambda **2 - (A+C) lambda + (A*C - B*B)
 172 *        and we solve it using the quadratic formula.
 173 *
 174          S = ( A+C )*HALF
 175          T = ( A-C )*HALF
 176 *
 177 *        Take the square root carefully to avoid over/under flow.
 178 *
 179          BABS = ABS( B )
 180          TABS = ABS( T )
 181          Z = MAX( BABS, TABS )
 182          IF( Z.GT.ZERO )
 183      $      T = Z*SQRT( ( T / Z )**2+( B / Z )**2 )
 184 *
 185 *        Compute the two eigenvalues.  RT1 and RT2 are exchanged
 186 *        if necessary so that RT1 will have the greater magnitude.
 187 *
 188          RT1 = S + T
 189          RT2 = S - T
 190          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
 191             TMP = RT1
 192             RT1 = RT2
 193             RT2 = TMP
 194          END IF
 195 *
 196 *        Choose CS1 = 1 and SN1 to satisfy the first equation, then
 197 *        scale the components of this eigenvector so that the matrix
 198 *        of eigenvectors X satisfies  X * X**T = I .  (No scaling is
 199 *        done if the norm of the eigenvalue matrix is less than THRESH.)
 200 *
 201          SN1 = ( RT1-A ) / B
 202          TABS = ABS( SN1 )
 203          IF( TABS.GT.ONE ) THEN
 204             T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 )
 205          ELSE
 206             T = SQRT( CONE+SN1*SN1 )
 207          END IF
 208          EVNORM = ABS( T )
 209          IF( EVNORM.GE.THRESH ) THEN
 210             EVSCAL = CONE / T
 211             CS1 = EVSCAL
 212             SN1 = SN1*EVSCAL
 213          ELSE
 214             EVSCAL = ZERO
 215          END IF
 216       END IF
 217       RETURN
 218 *
 219 *     End of CLAESY
 220 *
 221       END