LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
claesy.f
Go to the documentation of this file.
00001 *> \brief \b CLAESY
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLAESY + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claesy.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claesy.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claesy.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
00025 *       ..
00026 *  
00027 *
00028 *> \par Purpose:
00029 *  =============
00030 *>
00031 *> \verbatim
00032 *>
00033 *> CLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
00034 *>    ( ( A, B );( B, C ) )
00035 *> provided the norm of the matrix of eigenvectors is larger than
00036 *> some threshold value.
00037 *>
00038 *> RT1 is the eigenvalue of larger absolute value, and RT2 of
00039 *> smaller absolute value.  If the eigenvectors are computed, then
00040 *> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
00041 *>
00042 *> [  CS1     SN1   ] . [ A  B ] . [ CS1    -SN1   ] = [ RT1  0  ]
00043 *> [ -SN1     CS1   ]   [ B  C ]   [ SN1     CS1   ]   [  0  RT2 ]
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] A
00050 *> \verbatim
00051 *>          A is COMPLEX
00052 *>          The ( 1, 1 ) element of input matrix.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] B
00056 *> \verbatim
00057 *>          B is COMPLEX
00058 *>          The ( 1, 2 ) element of input matrix.  The ( 2, 1 ) element
00059 *>          is also given by B, since the 2-by-2 matrix is symmetric.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] C
00063 *> \verbatim
00064 *>          C is COMPLEX
00065 *>          The ( 2, 2 ) element of input matrix.
00066 *> \endverbatim
00067 *>
00068 *> \param[out] RT1
00069 *> \verbatim
00070 *>          RT1 is COMPLEX
00071 *>          The eigenvalue of larger modulus.
00072 *> \endverbatim
00073 *>
00074 *> \param[out] RT2
00075 *> \verbatim
00076 *>          RT2 is COMPLEX
00077 *>          The eigenvalue of smaller modulus.
00078 *> \endverbatim
00079 *>
00080 *> \param[out] EVSCAL
00081 *> \verbatim
00082 *>          EVSCAL is COMPLEX
00083 *>          The complex value by which the eigenvector matrix was scaled
00084 *>          to make it orthonormal.  If EVSCAL is zero, the eigenvectors
00085 *>          were not computed.  This means one of two things:  the 2-by-2
00086 *>          matrix could not be diagonalized, or the norm of the matrix
00087 *>          of eigenvectors before scaling was larger than the threshold
00088 *>          value THRESH (set below).
00089 *> \endverbatim
00090 *>
00091 *> \param[out] CS1
00092 *> \verbatim
00093 *>          CS1 is COMPLEX
00094 *> \endverbatim
00095 *>
00096 *> \param[out] SN1
00097 *> \verbatim
00098 *>          SN1 is COMPLEX
00099 *>          If EVSCAL .NE. 0,  ( CS1, SN1 ) is the unit right eigenvector
00100 *>          for RT1.
00101 *> \endverbatim
00102 *
00103 *  Authors:
00104 *  ========
00105 *
00106 *> \author Univ. of Tennessee 
00107 *> \author Univ. of California Berkeley 
00108 *> \author Univ. of Colorado Denver 
00109 *> \author NAG Ltd. 
00110 *
00111 *> \date November 2011
00112 *
00113 *> \ingroup complexSYauxiliary
00114 *
00115 *  =====================================================================
00116       SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
00117 *
00118 *  -- LAPACK auxiliary routine (version 3.4.0) --
00119 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00120 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00121 *     November 2011
00122 *
00123 *     .. Scalar Arguments ..
00124       COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
00125 *     ..
00126 *
00127 * =====================================================================
00128 *
00129 *     .. Parameters ..
00130       REAL               ZERO
00131       PARAMETER          ( ZERO = 0.0E0 )
00132       REAL               ONE
00133       PARAMETER          ( ONE = 1.0E0 )
00134       COMPLEX            CONE
00135       PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
00136       REAL               HALF
00137       PARAMETER          ( HALF = 0.5E0 )
00138       REAL               THRESH
00139       PARAMETER          ( THRESH = 0.1E0 )
00140 *     ..
00141 *     .. Local Scalars ..
00142       REAL               BABS, EVNORM, TABS, Z
00143       COMPLEX            S, T, TMP
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC          ABS, MAX, SQRT
00147 *     ..
00148 *     .. Executable Statements ..
00149 *
00150 *
00151 *     Special case:  The matrix is actually diagonal.
00152 *     To avoid divide by zero later, we treat this case separately.
00153 *
00154       IF( ABS( B ).EQ.ZERO ) THEN
00155          RT1 = A
00156          RT2 = C
00157          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
00158             TMP = RT1
00159             RT1 = RT2
00160             RT2 = TMP
00161             CS1 = ZERO
00162             SN1 = ONE
00163          ELSE
00164             CS1 = ONE
00165             SN1 = ZERO
00166          END IF
00167       ELSE
00168 *
00169 *        Compute the eigenvalues and eigenvectors.
00170 *        The characteristic equation is
00171 *           lambda **2 - (A+C) lambda + (A*C - B*B)
00172 *        and we solve it using the quadratic formula.
00173 *
00174          S = ( A+C )*HALF
00175          T = ( A-C )*HALF
00176 *
00177 *        Take the square root carefully to avoid over/under flow.
00178 *
00179          BABS = ABS( B )
00180          TABS = ABS( T )
00181          Z = MAX( BABS, TABS )
00182          IF( Z.GT.ZERO )
00183      $      T = Z*SQRT( ( T / Z )**2+( B / Z )**2 )
00184 *
00185 *        Compute the two eigenvalues.  RT1 and RT2 are exchanged
00186 *        if necessary so that RT1 will have the greater magnitude.
00187 *
00188          RT1 = S + T
00189          RT2 = S - T
00190          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
00191             TMP = RT1
00192             RT1 = RT2
00193             RT2 = TMP
00194          END IF
00195 *
00196 *        Choose CS1 = 1 and SN1 to satisfy the first equation, then
00197 *        scale the components of this eigenvector so that the matrix
00198 *        of eigenvectors X satisfies  X * X**T = I .  (No scaling is
00199 *        done if the norm of the eigenvalue matrix is less than THRESH.)
00200 *
00201          SN1 = ( RT1-A ) / B
00202          TABS = ABS( SN1 )
00203          IF( TABS.GT.ONE ) THEN
00204             T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 )
00205          ELSE
00206             T = SQRT( CONE+SN1*SN1 )
00207          END IF
00208          EVNORM = ABS( T )
00209          IF( EVNORM.GE.THRESH ) THEN
00210             EVSCAL = CONE / T
00211             CS1 = EVSCAL
00212             SN1 = SN1*EVSCAL
00213          ELSE
00214             EVSCAL = ZERO
00215          END IF
00216       END IF
00217       RETURN
00218 *
00219 *     End of CLAESY
00220 *
00221       END
 All Files Functions