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