Return to dlas2.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dlas2.f, revision 1.4

1.1       bertrand    1:       SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX )
                      2: *
                      3: *  -- LAPACK auxiliary routine (version 3.2) --
                      4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      6: *     November 2006
                      7: *
                      8: *     .. Scalar Arguments ..
                      9:       DOUBLE PRECISION   F, G, H, SSMAX, SSMIN
                     10: *     ..
                     11: *
                     12: *  Purpose
                     13: *  =======
                     14: *
                     15: *  DLAS2  computes the singular values of the 2-by-2 matrix
                     16: *     [  F   G  ]
                     17: *     [  0   H  ].
                     18: *  On return, SSMIN is the smaller singular value and SSMAX is the
                     19: *  larger singular value.
                     20: *
                     21: *  Arguments
                     22: *  =========
                     23: *
                     24: *  F       (input) DOUBLE PRECISION
                     25: *          The (1,1) element of the 2-by-2 matrix.
                     26: *
                     27: *  G       (input) DOUBLE PRECISION
                     28: *          The (1,2) element of the 2-by-2 matrix.
                     29: *
                     30: *  H       (input) DOUBLE PRECISION
                     31: *          The (2,2) element of the 2-by-2 matrix.
                     32: *
                     33: *  SSMIN   (output) DOUBLE PRECISION
                     34: *          The smaller singular value.
                     35: *
                     36: *  SSMAX   (output) DOUBLE PRECISION
                     37: *          The larger singular value.
                     38: *
                     39: *  Further Details
                     40: *  ===============
                     41: *
                     42: *  Barring over/underflow, all output quantities are correct to within
                     43: *  a few units in the last place (ulps), even in the absence of a guard
                     44: *  digit in addition/subtraction.
                     45: *
                     46: *  In IEEE arithmetic, the code works correctly if one matrix element is
                     47: *  infinite.
                     48: *
                     49: *  Overflow will not occur unless the largest singular value itself
                     50: *  overflows, or is within a few ulps of overflow. (On machines with
                     51: *  partial overflow, like the Cray, overflow may occur if the largest
                     52: *  singular value is within a factor of 2 of overflow.)
                     53: *
                     54: *  Underflow is harmless if underflow is gradual. Otherwise, results
                     55: *  may correspond to a matrix modified by perturbations of size near
                     56: *  the underflow threshold.
                     57: *
                     58: *  ====================================================================
                     59: *
                     60: *     .. Parameters ..
                     61:       DOUBLE PRECISION   ZERO
                     62:       PARAMETER          ( ZERO = 0.0D0 )
                     63:       DOUBLE PRECISION   ONE
                     64:       PARAMETER          ( ONE = 1.0D0 )
                     65:       DOUBLE PRECISION   TWO
                     66:       PARAMETER          ( TWO = 2.0D0 )
                     67: *     ..
                     68: *     .. Local Scalars ..
                     69:       DOUBLE PRECISION   AS, AT, AU, C, FA, FHMN, FHMX, GA, HA
                     70: *     ..
                     71: *     .. Intrinsic Functions ..
                     72:       INTRINSIC          ABS, MAX, MIN, SQRT
                     73: *     ..
                     74: *     .. Executable Statements ..
                     75: *
                     76:       FA = ABS( F )
                     77:       GA = ABS( G )
                     78:       HA = ABS( H )
                     79:       FHMN = MIN( FA, HA )
                     80:       FHMX = MAX( FA, HA )
                     81:       IF( FHMN.EQ.ZERO ) THEN
                     82:          SSMIN = ZERO
                     83:          IF( FHMX.EQ.ZERO ) THEN
                     84:             SSMAX = GA
                     85:          ELSE
                     86:             SSMAX = MAX( FHMX, GA )*SQRT( ONE+
                     87:      $              ( MIN( FHMX, GA ) / MAX( FHMX, GA ) )**2 )
                     88:          END IF
                     89:       ELSE
                     90:          IF( GA.LT.FHMX ) THEN
                     91:             AS = ONE + FHMN / FHMX
                     92:             AT = ( FHMX-FHMN ) / FHMX
                     93:             AU = ( GA / FHMX )**2
                     94:             C = TWO / ( SQRT( AS*AS+AU )+SQRT( AT*AT+AU ) )
                     95:             SSMIN = FHMN*C
                     96:             SSMAX = FHMX / C
                     97:          ELSE
                     98:             AU = FHMX / GA
                     99:             IF( AU.EQ.ZERO ) THEN
                    100: *
                    101: *              Avoid possible harmful underflow if exponent range
                    102: *              asymmetric (true SSMIN may not underflow even if
                    103: *              AU underflows)
                    104: *
                    105:                SSMIN = ( FHMN*FHMX ) / GA
                    106:                SSMAX = GA
                    107:             ELSE
                    108:                AS = ONE + FHMN / FHMX
                    109:                AT = ( FHMX-FHMN ) / FHMX
                    110:                C = ONE / ( SQRT( ONE+( AS*AU )**2 )+
                    111:      $             SQRT( ONE+( AT*AU )**2 ) )
                    112:                SSMIN = ( FHMN*C )*AU
                    113:                SSMIN = SSMIN + SSMIN
                    114:                SSMAX = GA / ( C+C )
                    115:             END IF
                    116:          END IF
                    117:       END IF
                    118:       RETURN
                    119: *
                    120: *     End of DLAS2
                    121: *
                    122:       END