LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaev2.f
Go to the documentation of this file.
00001 *> \brief \b DLAEV2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAEV2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
00025 *       ..
00026 *  
00027 *
00028 *> \par Purpose:
00029 *  =============
00030 *>
00031 *> \verbatim
00032 *>
00033 *> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
00034 *>    [  A   B  ]
00035 *>    [  B   C  ].
00036 *> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
00037 *> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
00038 *> eigenvector for RT1, giving the decomposition
00039 *>
00040 *>    [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
00041 *>    [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] A
00048 *> \verbatim
00049 *>          A is DOUBLE PRECISION
00050 *>          The (1,1) element of the 2-by-2 matrix.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] B
00054 *> \verbatim
00055 *>          B is DOUBLE PRECISION
00056 *>          The (1,2) element and the conjugate of the (2,1) element of
00057 *>          the 2-by-2 matrix.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] C
00061 *> \verbatim
00062 *>          C is DOUBLE PRECISION
00063 *>          The (2,2) element of the 2-by-2 matrix.
00064 *> \endverbatim
00065 *>
00066 *> \param[out] RT1
00067 *> \verbatim
00068 *>          RT1 is DOUBLE PRECISION
00069 *>          The eigenvalue of larger absolute value.
00070 *> \endverbatim
00071 *>
00072 *> \param[out] RT2
00073 *> \verbatim
00074 *>          RT2 is DOUBLE PRECISION
00075 *>          The eigenvalue of smaller absolute value.
00076 *> \endverbatim
00077 *>
00078 *> \param[out] CS1
00079 *> \verbatim
00080 *>          CS1 is DOUBLE PRECISION
00081 *> \endverbatim
00082 *>
00083 *> \param[out] SN1
00084 *> \verbatim
00085 *>          SN1 is DOUBLE PRECISION
00086 *>          The vector (CS1, SN1) is a unit right eigenvector for RT1.
00087 *> \endverbatim
00088 *
00089 *  Authors:
00090 *  ========
00091 *
00092 *> \author Univ. of Tennessee 
00093 *> \author Univ. of California Berkeley 
00094 *> \author Univ. of Colorado Denver 
00095 *> \author NAG Ltd. 
00096 *
00097 *> \date November 2011
00098 *
00099 *> \ingroup auxOTHERauxiliary
00100 *
00101 *> \par Further Details:
00102 *  =====================
00103 *>
00104 *> \verbatim
00105 *>
00106 *>  RT1 is accurate to a few ulps barring over/underflow.
00107 *>
00108 *>  RT2 may be inaccurate if there is massive cancellation in the
00109 *>  determinant A*C-B*B; higher precision or correctly rounded or
00110 *>  correctly truncated arithmetic would be needed to compute RT2
00111 *>  accurately in all cases.
00112 *>
00113 *>  CS1 and SN1 are accurate to a few ulps barring over/underflow.
00114 *>
00115 *>  Overflow is possible only if RT1 is within a factor of 5 of overflow.
00116 *>  Underflow is harmless if the input data is 0 or exceeds
00117 *>     underflow_threshold / macheps.
00118 *> \endverbatim
00119 *>
00120 *  =====================================================================
00121       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
00122 *
00123 *  -- LAPACK auxiliary routine (version 3.4.0) --
00124 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00125 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00126 *     November 2011
00127 *
00128 *     .. Scalar Arguments ..
00129       DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
00130 *     ..
00131 *
00132 * =====================================================================
00133 *
00134 *     .. Parameters ..
00135       DOUBLE PRECISION   ONE
00136       PARAMETER          ( ONE = 1.0D0 )
00137       DOUBLE PRECISION   TWO
00138       PARAMETER          ( TWO = 2.0D0 )
00139       DOUBLE PRECISION   ZERO
00140       PARAMETER          ( ZERO = 0.0D0 )
00141       DOUBLE PRECISION   HALF
00142       PARAMETER          ( HALF = 0.5D0 )
00143 *     ..
00144 *     .. Local Scalars ..
00145       INTEGER            SGN1, SGN2
00146       DOUBLE PRECISION   AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
00147      $                   TB, TN
00148 *     ..
00149 *     .. Intrinsic Functions ..
00150       INTRINSIC          ABS, SQRT
00151 *     ..
00152 *     .. Executable Statements ..
00153 *
00154 *     Compute the eigenvalues
00155 *
00156       SM = A + C
00157       DF = A - C
00158       ADF = ABS( DF )
00159       TB = B + B
00160       AB = ABS( TB )
00161       IF( ABS( A ).GT.ABS( C ) ) THEN
00162          ACMX = A
00163          ACMN = C
00164       ELSE
00165          ACMX = C
00166          ACMN = A
00167       END IF
00168       IF( ADF.GT.AB ) THEN
00169          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
00170       ELSE IF( ADF.LT.AB ) THEN
00171          RT = AB*SQRT( ONE+( ADF / AB )**2 )
00172       ELSE
00173 *
00174 *        Includes case AB=ADF=0
00175 *
00176          RT = AB*SQRT( TWO )
00177       END IF
00178       IF( SM.LT.ZERO ) THEN
00179          RT1 = HALF*( SM-RT )
00180          SGN1 = -1
00181 *
00182 *        Order of execution important.
00183 *        To get fully accurate smaller eigenvalue,
00184 *        next line needs to be executed in higher precision.
00185 *
00186          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
00187       ELSE IF( SM.GT.ZERO ) THEN
00188          RT1 = HALF*( SM+RT )
00189          SGN1 = 1
00190 *
00191 *        Order of execution important.
00192 *        To get fully accurate smaller eigenvalue,
00193 *        next line needs to be executed in higher precision.
00194 *
00195          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
00196       ELSE
00197 *
00198 *        Includes case RT1 = RT2 = 0
00199 *
00200          RT1 = HALF*RT
00201          RT2 = -HALF*RT
00202          SGN1 = 1
00203       END IF
00204 *
00205 *     Compute the eigenvector
00206 *
00207       IF( DF.GE.ZERO ) THEN
00208          CS = DF + RT
00209          SGN2 = 1
00210       ELSE
00211          CS = DF - RT
00212          SGN2 = -1
00213       END IF
00214       ACS = ABS( CS )
00215       IF( ACS.GT.AB ) THEN
00216          CT = -TB / CS
00217          SN1 = ONE / SQRT( ONE+CT*CT )
00218          CS1 = CT*SN1
00219       ELSE
00220          IF( AB.EQ.ZERO ) THEN
00221             CS1 = ONE
00222             SN1 = ZERO
00223          ELSE
00224             TN = -CS / TB
00225             CS1 = ONE / SQRT( ONE+TN*TN )
00226             SN1 = TN*CS1
00227          END IF
00228       END IF
00229       IF( SGN1.EQ.SGN2 ) THEN
00230          TN = CS1
00231          CS1 = -SN1
00232          SN1 = TN
00233       END IF
00234       RETURN
00235 *
00236 *     End of DLAEV2
00237 *
00238       END
 All Files Functions