dlas2.f

Go to the documentation of this file.

      SUBROUTINE dlas2( F, G, H, SSMIN, SSMAX )
*
*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     September 30, 1994
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   F, G, H, SSMAX, SSMIN
*     ..
*
*  Purpose
*  =======
*
*  DLAS2  computes the singular values of the 2-by-2 matrix
*     [  F   G  ]
*     [  0   H  ].
*  On return, SSMIN is the smaller singular value and SSMAX is the
*  larger singular value.
*
*  Arguments
*  =========
*
*  F       (input) DOUBLE PRECISION
*          The (1,1) element of the 2-by-2 matrix.
*
*  G       (input) DOUBLE PRECISION
*          The (1,2) element of the 2-by-2 matrix.
*
*  H       (input) DOUBLE PRECISION
*          The (2,2) element of the 2-by-2 matrix.
*
*  SSMIN   (output) DOUBLE PRECISION
*          The smaller singular value.
*
*  SSMAX   (output) DOUBLE PRECISION
*          The larger singular value.
*
*  Further Details
*  ===============
*
*  Barring over/underflow, all output quantities are correct to within
*  a few units in the last place (ulps), even in the absence of a guard
*  digit in addition/subtraction.
*
*  In IEEE arithmetic, the code works correctly if one matrix element is
*  infinite.
*
*  Overflow will not occur unless the largest singular value itself
*  overflows, or is within a few ulps of overflow. (On machines with
*  partial overflow, like the Cray, overflow may occur if the largest
*  singular value is within a factor of 2 of overflow.)
*
*  Underflow is harmless if underflow is gradual. Otherwise, results
*  may correspond to a matrix modified by perturbations of size near
*  the underflow threshold.
*
*  ====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d0 )
      DOUBLE PRECISION   ONE
      parameter( one = 1.0d0 )
      DOUBLE PRECISION   TWO
      parameter( two = 2.0d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   AS, AT, AU, C, FA, FHMN, FHMX, GA, HA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt
*     ..
*     .. Executable Statements ..
*
      fa = abs( f )
      ga = abs( g )
      ha = abs( h )
      fhmn = min( fa, ha )
      fhmx = max( fa, ha )
      IF( fhmn.EQ.zero ) THEN
         ssmin = zero
         IF( fhmx.EQ.zero ) THEN
            ssmax = ga
         ELSE
            ssmax = max( fhmx, ga )*sqrt( one+
     $              ( min( fhmx, ga ) / max( fhmx, ga ) )**2 )
         END IF
      ELSE
         IF( ga.LT.fhmx ) THEN
            as = one + fhmn / fhmx
            at = ( fhmx-fhmn ) / fhmx
            au = ( ga / fhmx )**2
            c = two / ( sqrt( as*as+au )+sqrt( at*at+au ) )
            ssmin = fhmn*c
            ssmax = fhmx / c
         ELSE
            au = fhmx / ga
            IF( au.EQ.zero ) THEN
*
*              Avoid possible harmful underflow if exponent range
*              asymmetric (true SSMIN may not underflow even if
*              AU underflows)
*
               ssmin = ( fhmn*fhmx ) / ga
               ssmax = ga
            ELSE
               as = one + fhmn / fhmx
               at = ( fhmx-fhmn ) / fhmx
               c = one / ( sqrt( one+( as*au )**2 )+
     $             sqrt( one+( at*au )**2 ) )
               ssmin = ( fhmn*c )*au
               ssmin = ssmin + ssmin
               ssmax = ga / ( c+c )
            END IF
         END IF
      END IF
      RETURN
*
*     End of DLAS2
*
      END