|
2329
|
1 SUBROUTINE DLARTG( F, G, CS, SN, R ) |
|
|
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 CS, F, G, R, SN |
|
|
9 * .. |
|
|
10 * |
|
|
11 * Purpose |
|
|
12 * ======= |
|
|
13 * |
|
|
14 * DLARTG generate a plane rotation so that |
|
|
15 * |
|
|
16 * [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. |
|
|
17 * [ -SN CS ] [ G ] [ 0 ] |
|
|
18 * |
|
|
19 * This is a slower, more accurate version of the BLAS1 routine DROTG, |
|
|
20 * with the following other differences: |
|
|
21 * F and G are unchanged on return. |
|
|
22 * If G=0, then CS=1 and SN=0. |
|
|
23 * If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any |
|
|
24 * floating point operations (saves work in DBDSQR when |
|
|
25 * there are zeros on the diagonal). |
|
|
26 * |
|
|
27 * If F exceeds G in magnitude, CS will be positive. |
|
|
28 * |
|
|
29 * Arguments |
|
|
30 * ========= |
|
|
31 * |
|
|
32 * F (input) DOUBLE PRECISION |
|
|
33 * The first component of vector to be rotated. |
|
|
34 * |
|
|
35 * G (input) DOUBLE PRECISION |
|
|
36 * The second component of vector to be rotated. |
|
|
37 * |
|
|
38 * CS (output) DOUBLE PRECISION |
|
|
39 * The cosine of the rotation. |
|
|
40 * |
|
|
41 * SN (output) DOUBLE PRECISION |
|
|
42 * The sine of the rotation. |
|
|
43 * |
|
|
44 * R (output) DOUBLE PRECISION |
|
|
45 * The nonzero component of the rotated vector. |
|
|
46 * |
|
7034
|
47 * This version has a few statements commented out for thread safety |
|
|
48 * (machine parameters are computed on each entry). 10 feb 03, SJH. |
|
|
49 * |
|
2329
|
50 * ===================================================================== |
|
|
51 * |
|
|
52 * .. Parameters .. |
|
|
53 DOUBLE PRECISION ZERO |
|
|
54 PARAMETER ( ZERO = 0.0D0 ) |
|
|
55 DOUBLE PRECISION ONE |
|
|
56 PARAMETER ( ONE = 1.0D0 ) |
|
|
57 DOUBLE PRECISION TWO |
|
|
58 PARAMETER ( TWO = 2.0D0 ) |
|
|
59 * .. |
|
|
60 * .. Local Scalars .. |
|
7034
|
61 * LOGICAL FIRST |
|
2329
|
62 INTEGER COUNT, I |
|
|
63 DOUBLE PRECISION EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE |
|
|
64 * .. |
|
|
65 * .. External Functions .. |
|
|
66 DOUBLE PRECISION DLAMCH |
|
|
67 EXTERNAL DLAMCH |
|
|
68 * .. |
|
|
69 * .. Intrinsic Functions .. |
|
|
70 INTRINSIC ABS, INT, LOG, MAX, SQRT |
|
|
71 * .. |
|
|
72 * .. Save statement .. |
|
7034
|
73 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 |
|
2329
|
74 * .. |
|
|
75 * .. Data statements .. |
|
7034
|
76 * DATA FIRST / .TRUE. / |
|
2329
|
77 * .. |
|
|
78 * .. Executable Statements .. |
|
|
79 * |
|
7034
|
80 * IF( FIRST ) THEN |
|
2329
|
81 SAFMIN = DLAMCH( 'S' ) |
|
|
82 EPS = DLAMCH( 'E' ) |
|
|
83 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / |
|
|
84 $ LOG( DLAMCH( 'B' ) ) / TWO ) |
|
|
85 SAFMX2 = ONE / SAFMN2 |
|
7034
|
86 * FIRST = .FALSE. |
|
|
87 * END IF |
|
2329
|
88 IF( G.EQ.ZERO ) THEN |
|
|
89 CS = ONE |
|
|
90 SN = ZERO |
|
|
91 R = F |
|
|
92 ELSE IF( F.EQ.ZERO ) THEN |
|
|
93 CS = ZERO |
|
|
94 SN = ONE |
|
|
95 R = G |
|
|
96 ELSE |
|
|
97 F1 = F |
|
|
98 G1 = G |
|
|
99 SCALE = MAX( ABS( F1 ), ABS( G1 ) ) |
|
|
100 IF( SCALE.GE.SAFMX2 ) THEN |
|
|
101 COUNT = 0 |
|
|
102 10 CONTINUE |
|
|
103 COUNT = COUNT + 1 |
|
|
104 F1 = F1*SAFMN2 |
|
|
105 G1 = G1*SAFMN2 |
|
|
106 SCALE = MAX( ABS( F1 ), ABS( G1 ) ) |
|
|
107 IF( SCALE.GE.SAFMX2 ) |
|
|
108 $ GO TO 10 |
|
|
109 R = SQRT( F1**2+G1**2 ) |
|
|
110 CS = F1 / R |
|
|
111 SN = G1 / R |
|
|
112 DO 20 I = 1, COUNT |
|
|
113 R = R*SAFMX2 |
|
|
114 20 CONTINUE |
|
|
115 ELSE IF( SCALE.LE.SAFMN2 ) THEN |
|
|
116 COUNT = 0 |
|
|
117 30 CONTINUE |
|
|
118 COUNT = COUNT + 1 |
|
|
119 F1 = F1*SAFMX2 |
|
|
120 G1 = G1*SAFMX2 |
|
|
121 SCALE = MAX( ABS( F1 ), ABS( G1 ) ) |
|
|
122 IF( SCALE.LE.SAFMN2 ) |
|
|
123 $ GO TO 30 |
|
|
124 R = SQRT( F1**2+G1**2 ) |
|
|
125 CS = F1 / R |
|
|
126 SN = G1 / R |
|
|
127 DO 40 I = 1, COUNT |
|
|
128 R = R*SAFMN2 |
|
|
129 40 CONTINUE |
|
|
130 ELSE |
|
|
131 R = SQRT( F1**2+G1**2 ) |
|
|
132 CS = F1 / R |
|
|
133 SN = G1 / R |
|
|
134 END IF |
|
|
135 IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN |
|
|
136 CS = -CS |
|
|
137 SN = -SN |
|
|
138 R = -R |
|
|
139 END IF |
|
|
140 END IF |
|
|
141 RETURN |
|
|
142 * |
|
|
143 * End of DLARTG |
|
|
144 * |
|
|
145 END |
