vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlasqtest/dlas2.f

   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