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

     1  *> \brief \b DLAS2 computes singular values of a 2-by-2 triangular matrix.
     2  *
     3  *  =========== DOCUMENTATION ===========
     4  *
     5  * Online html documentation available at 
     6  *            http://www.netlib.org/lapack/explore-html/ 
     7  *
     8  *> \htmlonly
     9  *> Download DLAS2 + dependencies 
    10  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlas2.f"> 
    11  *> [TGZ]</a> 
    12  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlas2.f"> 
    13  *> [ZIP]</a> 
    14  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlas2.f"> 
    15  *> [TXT]</a>
    16  *> \endhtmlonly 
    17  *
    18  *  Definition:
    19  *  ===========
    20  *
    21  *       SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX )
    22  * 
    23  *       .. Scalar Arguments ..
    24  *       DOUBLE PRECISION   F, G, H, SSMAX, SSMIN
    25  *       ..
    26  *  
    27  *
    28  *> \par Purpose:
    29  *  =============
    30  *>
    31  *> \verbatim
    32  *>
    33  *> DLAS2  computes the singular values of the 2-by-2 matrix
    34  *>    [  F   G  ]
    35  *>    [  0   H  ].
    36  *> On return, SSMIN is the smaller singular value and SSMAX is the
    37  *> larger singular value.
    38  *> \endverbatim
    39  *
    40  *  Arguments:
    41  *  ==========
    42  *
    43  *> \param[in] F
    44  *> \verbatim
    45  *>          F is DOUBLE PRECISION
    46  *>          The (1,1) element of the 2-by-2 matrix.
    47  *> \endverbatim
    48  *>
    49  *> \param[in] G
    50  *> \verbatim
    51  *>          G is DOUBLE PRECISION
    52  *>          The (1,2) element of the 2-by-2 matrix.
    53  *> \endverbatim
    54  *>
    55  *> \param[in] H
    56  *> \verbatim
    57  *>          H is DOUBLE PRECISION
    58  *>          The (2,2) element of the 2-by-2 matrix.
    59  *> \endverbatim
    60  *>
    61  *> \param[out] SSMIN
    62  *> \verbatim
    63  *>          SSMIN is DOUBLE PRECISION
    64  *>          The smaller singular value.
    65  *> \endverbatim
    66  *>
    67  *> \param[out] SSMAX
    68  *> \verbatim
    69  *>          SSMAX is DOUBLE PRECISION
    70  *>          The larger singular 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  *>  Barring over/underflow, all output quantities are correct to within
    91  *>  a few units in the last place (ulps), even in the absence of a guard
    92  *>  digit in addition/subtraction.
    93  *>
    94  *>  In IEEE arithmetic, the code works correctly if one matrix element is
    95  *>  infinite.
    96  *>
    97  *>  Overflow will not occur unless the largest singular value itself
    98  *>  overflows, or is within a few ulps of overflow. (On machines with
    99  *>  partial overflow, like the Cray, overflow may occur if the largest
   100  *>  singular value is within a factor of 2 of overflow.)
   101  *>
   102  *>  Underflow is harmless if underflow is gradual. Otherwise, results
   103  *>  may correspond to a matrix modified by perturbations of size near
   104  *>  the underflow threshold.
   105  *> \endverbatim
   106  *>
   107  *  =====================================================================
   108        SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX )
   109  *
   110  *  -- LAPACK auxiliary routine (version 3.4.2) --
   111  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   112  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   113  *     September 2012
   114  *
   115  *     .. Scalar Arguments ..
   116        DOUBLE PRECISION   F, G, H, SSMAX, SSMIN
   117  *     ..
   118  *
   119  *  ====================================================================
   120  *
   121  *     .. Parameters ..
   122        DOUBLE PRECISION   ZERO
   123        PARAMETER          ( ZERO = 0.0D0 )
   124        DOUBLE PRECISION   ONE
   125        PARAMETER          ( ONE = 1.0D0 )
   126        DOUBLE PRECISION   TWO
   127        PARAMETER          ( TWO = 2.0D0 )
   128  *     ..
   129  *     .. Local Scalars ..
   130        DOUBLE PRECISION   AS, AT, AU, C, FA, FHMN, FHMX, GA, HA
   131  *     ..
   132  *     .. Intrinsic Functions ..
   133        INTRINSIC          ABS, MAX, MIN, SQRT
   134  *     ..
   135  *     .. Executable Statements ..
   136  *
   137        FA = ABS( F )
   138        GA = ABS( G )
   139        HA = ABS( H )
   140        FHMN = MIN( FA, HA )
   141        FHMX = MAX( FA, HA )
   142        IF( FHMN.EQ.ZERO ) THEN
   143           SSMIN = ZERO
   144           IF( FHMX.EQ.ZERO ) THEN
   145              SSMAX = GA
   146           ELSE
   147              SSMAX = MAX( FHMX, GA )*SQRT( ONE+
   148       $              ( MIN( FHMX, GA ) / MAX( FHMX, GA ) )**2 )
   149           END IF
   150        ELSE
   151           IF( GA.LT.FHMX ) THEN
   152              AS = ONE + FHMN / FHMX
   153              AT = ( FHMX-FHMN ) / FHMX
   154              AU = ( GA / FHMX )**2
   155              C = TWO / ( SQRT( AS*AS+AU )+SQRT( AT*AT+AU ) )
   156              SSMIN = FHMN*C
   157              SSMAX = FHMX / C
   158           ELSE
   159              AU = FHMX / GA
   160              IF( AU.EQ.ZERO ) THEN
   161  *
   162  *              Avoid possible harmful underflow if exponent range
   163  *              asymmetric (true SSMIN may not underflow even if
   164  *              AU underflows)
   165  *
   166                 SSMIN = ( FHMN*FHMX ) / GA
   167                 SSMAX = GA
   168              ELSE
   169                 AS = ONE + FHMN / FHMX
   170                 AT = ( FHMX-FHMN ) / FHMX
   171                 C = ONE / ( SQRT( ONE+( AS*AU )**2 )+
   172       $             SQRT( ONE+( AT*AU )**2 ) )
   173                 SSMIN = ( FHMN*C )*AU
   174                 SSMIN = SSMIN + SSMIN
   175                 SSMAX = GA / ( C+C )
   176              END IF
   177           END IF
   178        END IF
   179        RETURN
   180  *
   181  *     End of DLAS2
   182  *
   183        END