LAPACK 3.3.0
|
00001 SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 DOUBLE PRECISION F, G, H, SSMAX, SSMIN 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * DLAS2 computes the singular values of the 2-by-2 matrix 00016 * [ F G ] 00017 * [ 0 H ]. 00018 * On return, SSMIN is the smaller singular value and SSMAX is the 00019 * larger singular value. 00020 * 00021 * Arguments 00022 * ========= 00023 * 00024 * F (input) DOUBLE PRECISION 00025 * The (1,1) element of the 2-by-2 matrix. 00026 * 00027 * G (input) DOUBLE PRECISION 00028 * The (1,2) element of the 2-by-2 matrix. 00029 * 00030 * H (input) DOUBLE PRECISION 00031 * The (2,2) element of the 2-by-2 matrix. 00032 * 00033 * SSMIN (output) DOUBLE PRECISION 00034 * The smaller singular value. 00035 * 00036 * SSMAX (output) DOUBLE PRECISION 00037 * The larger singular value. 00038 * 00039 * Further Details 00040 * =============== 00041 * 00042 * Barring over/underflow, all output quantities are correct to within 00043 * a few units in the last place (ulps), even in the absence of a guard 00044 * digit in addition/subtraction. 00045 * 00046 * In IEEE arithmetic, the code works correctly if one matrix element is 00047 * infinite. 00048 * 00049 * Overflow will not occur unless the largest singular value itself 00050 * overflows, or is within a few ulps of overflow. (On machines with 00051 * partial overflow, like the Cray, overflow may occur if the largest 00052 * singular value is within a factor of 2 of overflow.) 00053 * 00054 * Underflow is harmless if underflow is gradual. Otherwise, results 00055 * may correspond to a matrix modified by perturbations of size near 00056 * the underflow threshold. 00057 * 00058 * ==================================================================== 00059 * 00060 * .. Parameters .. 00061 DOUBLE PRECISION ZERO 00062 PARAMETER ( ZERO = 0.0D0 ) 00063 DOUBLE PRECISION ONE 00064 PARAMETER ( ONE = 1.0D0 ) 00065 DOUBLE PRECISION TWO 00066 PARAMETER ( TWO = 2.0D0 ) 00067 * .. 00068 * .. Local Scalars .. 00069 DOUBLE PRECISION AS, AT, AU, C, FA, FHMN, FHMX, GA, HA 00070 * .. 00071 * .. Intrinsic Functions .. 00072 INTRINSIC ABS, MAX, MIN, SQRT 00073 * .. 00074 * .. Executable Statements .. 00075 * 00076 FA = ABS( F ) 00077 GA = ABS( G ) 00078 HA = ABS( H ) 00079 FHMN = MIN( FA, HA ) 00080 FHMX = MAX( FA, HA ) 00081 IF( FHMN.EQ.ZERO ) THEN 00082 SSMIN = ZERO 00083 IF( FHMX.EQ.ZERO ) THEN 00084 SSMAX = GA 00085 ELSE 00086 SSMAX = MAX( FHMX, GA )*SQRT( ONE+ 00087 $ ( MIN( FHMX, GA ) / MAX( FHMX, GA ) )**2 ) 00088 END IF 00089 ELSE 00090 IF( GA.LT.FHMX ) THEN 00091 AS = ONE + FHMN / FHMX 00092 AT = ( FHMX-FHMN ) / FHMX 00093 AU = ( GA / FHMX )**2 00094 C = TWO / ( SQRT( AS*AS+AU )+SQRT( AT*AT+AU ) ) 00095 SSMIN = FHMN*C 00096 SSMAX = FHMX / C 00097 ELSE 00098 AU = FHMX / GA 00099 IF( AU.EQ.ZERO ) THEN 00100 * 00101 * Avoid possible harmful underflow if exponent range 00102 * asymmetric (true SSMIN may not underflow even if 00103 * AU underflows) 00104 * 00105 SSMIN = ( FHMN*FHMX ) / GA 00106 SSMAX = GA 00107 ELSE 00108 AS = ONE + FHMN / FHMX 00109 AT = ( FHMX-FHMN ) / FHMX 00110 C = ONE / ( SQRT( ONE+( AS*AU )**2 )+ 00111 $ SQRT( ONE+( AT*AU )**2 ) ) 00112 SSMIN = ( FHMN*C )*AU 00113 SSMIN = SSMIN + SSMIN 00114 SSMAX = GA / ( C+C ) 00115 END IF 00116 END IF 00117 END IF 00118 RETURN 00119 * 00120 * End of DLAS2 00121 * 00122 END