github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dsterftest/dlae2.f (about)

     1  *> \brief \b DLAE2 computes the eigenvalues of a 2-by-2 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 DLAE2 + dependencies 
    10  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlae2.f"> 
    11  *> [TGZ]</a> 
    12  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlae2.f"> 
    13  *> [ZIP]</a> 
    14  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlae2.f"> 
    15  *> [TXT]</a>
    16  *> \endhtmlonly 
    17  *
    18  *  Definition:
    19  *  ===========
    20  *
    21  *       SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
    22  * 
    23  *       .. Scalar Arguments ..
    24  *       DOUBLE PRECISION   A, B, C, RT1, RT2
    25  *       ..
    26  *  
    27  *
    28  *> \par Purpose:
    29  *  =============
    30  *>
    31  *> \verbatim
    32  *>
    33  *> DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
    34  *>    [  A   B  ]
    35  *>    [  B   C  ].
    36  *> On return, RT1 is the eigenvalue of larger absolute value, and RT2
    37  *> is the eigenvalue of smaller absolute value.
    38  *> \endverbatim
    39  *
    40  *  Arguments:
    41  *  ==========
    42  *
    43  *> \param[in] A
    44  *> \verbatim
    45  *>          A is DOUBLE PRECISION
    46  *>          The (1,1) element of the 2-by-2 matrix.
    47  *> \endverbatim
    48  *>
    49  *> \param[in] B
    50  *> \verbatim
    51  *>          B is DOUBLE PRECISION
    52  *>          The (1,2) and (2,1) elements of the 2-by-2 matrix.
    53  *> \endverbatim
    54  *>
    55  *> \param[in] C
    56  *> \verbatim
    57  *>          C is DOUBLE PRECISION
    58  *>          The (2,2) element of the 2-by-2 matrix.
    59  *> \endverbatim
    60  *>
    61  *> \param[out] RT1
    62  *> \verbatim
    63  *>          RT1 is DOUBLE PRECISION
    64  *>          The eigenvalue of larger absolute value.
    65  *> \endverbatim
    66  *>
    67  *> \param[out] RT2
    68  *> \verbatim
    69  *>          RT2 is DOUBLE PRECISION
    70  *>          The eigenvalue of smaller absolute value.
    71  *> \endverbatim
    72  *
    73  *  Authors:
    74  *  ========
    75  *
    76  *> \author Univ. of Tennessee 
    77  *> \author Univ. of California Berkeley 
    78  *> \author Univ. of Colorado Denver 
    79  *> \author NAG Ltd. 
    80  *
    81  *> \date September 2012
    82  *
    83  *> \ingroup auxOTHERauxiliary
    84  *
    85  *> \par Further Details:
    86  *  =====================
    87  *>
    88  *> \verbatim
    89  *>
    90  *>  RT1 is accurate to a few ulps barring over/underflow.
    91  *>
    92  *>  RT2 may be inaccurate if there is massive cancellation in the
    93  *>  determinant A*C-B*B; higher precision or correctly rounded or
    94  *>  correctly truncated arithmetic would be needed to compute RT2
    95  *>  accurately in all cases.
    96  *>
    97  *>  Overflow is possible only if RT1 is within a factor of 5 of overflow.
    98  *>  Underflow is harmless if the input data is 0 or exceeds
    99  *>     underflow_threshold / macheps.
   100  *> \endverbatim
   101  *>
   102  *  =====================================================================
   103        SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
   104  *
   105  *  -- LAPACK auxiliary routine (version 3.4.2) --
   106  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   107  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   108  *     September 2012
   109  *
   110  *     .. Scalar Arguments ..
   111        DOUBLE PRECISION   A, B, C, RT1, RT2
   112  *     ..
   113  *
   114  * =====================================================================
   115  *
   116  *     .. Parameters ..
   117        DOUBLE PRECISION   ONE
   118        PARAMETER          ( ONE = 1.0D0 )
   119        DOUBLE PRECISION   TWO
   120        PARAMETER          ( TWO = 2.0D0 )
   121        DOUBLE PRECISION   ZERO
   122        PARAMETER          ( ZERO = 0.0D0 )
   123        DOUBLE PRECISION   HALF
   124        PARAMETER          ( HALF = 0.5D0 )
   125  *     ..
   126  *     .. Local Scalars ..
   127        DOUBLE PRECISION   AB, ACMN, ACMX, ADF, DF, RT, SM, TB
   128  *     ..
   129  *     .. Intrinsic Functions ..
   130        INTRINSIC          ABS, SQRT
   131  *     ..
   132  *     .. Executable Statements ..
   133  *
   134  *     Compute the eigenvalues
   135  *
   136        SM = A + C
   137        DF = A - C
   138        ADF = ABS( DF )
   139        TB = B + B
   140        AB = ABS( TB )
   141        IF( ABS( A ).GT.ABS( C ) ) THEN
   142           ACMX = A
   143           ACMN = C
   144        ELSE
   145           ACMX = C
   146           ACMN = A
   147        END IF
   148        IF( ADF.GT.AB ) THEN
   149           RT = ADF*SQRT( ONE+( AB / ADF )**2 )
   150        ELSE IF( ADF.LT.AB ) THEN
   151           RT = AB*SQRT( ONE+( ADF / AB )**2 )
   152        ELSE
   153  *
   154  *        Includes case AB=ADF=0
   155  *
   156           RT = AB*SQRT( TWO )
   157        END IF
   158        IF( SM.LT.ZERO ) THEN
   159           RT1 = HALF*( SM-RT )
   160  *
   161  *        Order of execution important.
   162  *        To get fully accurate smaller eigenvalue,
   163  *        next line needs to be executed in higher precision.
   164  *
   165           RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
   166        ELSE IF( SM.GT.ZERO ) THEN
   167           RT1 = HALF*( SM+RT )
   168  *
   169  *        Order of execution important.
   170  *        To get fully accurate smaller eigenvalue,
   171  *        next line needs to be executed in higher precision.
   172  *
   173           RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
   174        ELSE
   175  *
   176  *        Includes case RT1 = RT2 = 0
   177  *
   178           RT1 = HALF*RT
   179           RT2 = -HALF*RT
   180        END IF
   181        RETURN
   182  *
   183  *     End of DLAE2
   184  *
   185        END