diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/dlas2.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/dlas2.f | 121 |
1 files changed, 121 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dlas2.f b/2.3-1/src/fortran/lapack/dlas2.f new file mode 100644 index 00000000..e100a4d8 --- /dev/null +++ b/2.3-1/src/fortran/lapack/dlas2.f @@ -0,0 +1,121 @@ + SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. 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 |