LAPACK 3.3.1
Linear Algebra PACKage

claesy.f

Go to the documentation of this file.
00001       SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  CLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
00016 *     ( ( A, B );( B, C ) )
00017 *  provided the norm of the matrix of eigenvectors is larger than
00018 *  some threshold value.
00019 *
00020 *  RT1 is the eigenvalue of larger absolute value, and RT2 of
00021 *  smaller absolute value.  If the eigenvectors are computed, then
00022 *  on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
00023 *
00024 *  [  CS1     SN1   ] . [ A  B ] . [ CS1    -SN1   ] = [ RT1  0  ]
00025 *  [ -SN1     CS1   ]   [ B  C ]   [ SN1     CS1   ]   [  0  RT2 ]
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  A       (input) COMPLEX
00031 *          The ( 1, 1 ) element of input matrix.
00032 *
00033 *  B       (input) COMPLEX
00034 *          The ( 1, 2 ) element of input matrix.  The ( 2, 1 ) element
00035 *          is also given by B, since the 2-by-2 matrix is symmetric.
00036 *
00037 *  C       (input) COMPLEX
00038 *          The ( 2, 2 ) element of input matrix.
00039 *
00040 *  RT1     (output) COMPLEX
00041 *          The eigenvalue of larger modulus.
00042 *
00043 *  RT2     (output) COMPLEX
00044 *          The eigenvalue of smaller modulus.
00045 *
00046 *  EVSCAL  (output) COMPLEX
00047 *          The complex value by which the eigenvector matrix was scaled
00048 *          to make it orthonormal.  If EVSCAL is zero, the eigenvectors
00049 *          were not computed.  This means one of two things:  the 2-by-2
00050 *          matrix could not be diagonalized, or the norm of the matrix
00051 *          of eigenvectors before scaling was larger than the threshold
00052 *          value THRESH (set below).
00053 *
00054 *  CS1     (output) COMPLEX
00055 *  SN1     (output) COMPLEX
00056 *          If EVSCAL .NE. 0,  ( CS1, SN1 ) is the unit right eigenvector
00057 *          for RT1.
00058 *
00059 * =====================================================================
00060 *
00061 *     .. Parameters ..
00062       REAL               ZERO
00063       PARAMETER          ( ZERO = 0.0E0 )
00064       REAL               ONE
00065       PARAMETER          ( ONE = 1.0E0 )
00066       COMPLEX            CONE
00067       PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
00068       REAL               HALF
00069       PARAMETER          ( HALF = 0.5E0 )
00070       REAL               THRESH
00071       PARAMETER          ( THRESH = 0.1E0 )
00072 *     ..
00073 *     .. Local Scalars ..
00074       REAL               BABS, EVNORM, TABS, Z
00075       COMPLEX            S, T, TMP
00076 *     ..
00077 *     .. Intrinsic Functions ..
00078       INTRINSIC          ABS, MAX, SQRT
00079 *     ..
00080 *     .. Executable Statements ..
00081 *
00082 *
00083 *     Special case:  The matrix is actually diagonal.
00084 *     To avoid divide by zero later, we treat this case separately.
00085 *
00086       IF( ABS( B ).EQ.ZERO ) THEN
00087          RT1 = A
00088          RT2 = C
00089          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
00090             TMP = RT1
00091             RT1 = RT2
00092             RT2 = TMP
00093             CS1 = ZERO
00094             SN1 = ONE
00095          ELSE
00096             CS1 = ONE
00097             SN1 = ZERO
00098          END IF
00099       ELSE
00100 *
00101 *        Compute the eigenvalues and eigenvectors.
00102 *        The characteristic equation is
00103 *           lambda **2 - (A+C) lambda + (A*C - B*B)
00104 *        and we solve it using the quadratic formula.
00105 *
00106          S = ( A+C )*HALF
00107          T = ( A-C )*HALF
00108 *
00109 *        Take the square root carefully to avoid over/under flow.
00110 *
00111          BABS = ABS( B )
00112          TABS = ABS( T )
00113          Z = MAX( BABS, TABS )
00114          IF( Z.GT.ZERO )
00115      $      T = Z*SQRT( ( T / Z )**2+( B / Z )**2 )
00116 *
00117 *        Compute the two eigenvalues.  RT1 and RT2 are exchanged
00118 *        if necessary so that RT1 will have the greater magnitude.
00119 *
00120          RT1 = S + T
00121          RT2 = S - T
00122          IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
00123             TMP = RT1
00124             RT1 = RT2
00125             RT2 = TMP
00126          END IF
00127 *
00128 *        Choose CS1 = 1 and SN1 to satisfy the first equation, then
00129 *        scale the components of this eigenvector so that the matrix
00130 *        of eigenvectors X satisfies  X * X**T = I .  (No scaling is
00131 *        done if the norm of the eigenvalue matrix is less than THRESH.)
00132 *
00133          SN1 = ( RT1-A ) / B
00134          TABS = ABS( SN1 )
00135          IF( TABS.GT.ONE ) THEN
00136             T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 )
00137          ELSE
00138             T = SQRT( CONE+SN1*SN1 )
00139          END IF
00140          EVNORM = ABS( T )
00141          IF( EVNORM.GE.THRESH ) THEN
00142             EVSCAL = CONE / T
00143             CS1 = EVSCAL
00144             SN1 = SN1*EVSCAL
00145          ELSE
00146             EVSCAL = ZERO
00147          END IF
00148       END IF
00149       RETURN
00150 *
00151 *     End of CLAESY
00152 *
00153       END
 All Files Functions