1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
#include "f2c.h"
/* Subroutine */ int dlartg_(doublereal *f, doublereal *g, doublereal *cs,
doublereal *sn, doublereal *r)
{
/* -- LAPACK auxiliary routine (version 2.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
September 30, 1994
Purpose
=======
DLARTG generate a plane rotation so that
[ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1.
[ -SN CS ] [ G ] [ 0 ]
This is a slower, more accurate version of the BLAS1 routine DROTG,
with the following other differences:
F and G are unchanged on return.
If G=0, then CS=1 and SN=0.
If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
floating point operations (saves work in DBDSQR when
there are zeros on the diagonal).
If F exceeds G in magnitude, CS will be positive.
Arguments
=========
F (input) DOUBLE PRECISION
The first component of vector to be rotated.
G (input) DOUBLE PRECISION
The second component of vector to be rotated.
CS (output) DOUBLE PRECISION
The cosine of the rotation.
SN (output) DOUBLE PRECISION
The sine of the rotation.
R (output) DOUBLE PRECISION
The nonzero component of the rotated vector.
=====================================================================
*/
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer i__1;
doublereal d__1, d__2;
/* Builtin functions */
double log(doublereal), pow_di(doublereal *, integer *), sqrt(doublereal);
/* Local variables */
static integer i;
static doublereal scale;
static integer count;
static doublereal f1, g1, safmn2, safmx2;
extern doublereal dlamch_(char *);
static doublereal safmin, eps;
if (first) {
first = FALSE_;
safmin = dlamch_("S");
eps = dlamch_("E");
d__1 = dlamch_("B");
i__1 = (integer) (log(safmin / eps) / log(dlamch_("B")) / 2.);
safmn2 = pow_di(&d__1, &i__1);
safmx2 = 1. / safmn2;
}
if (*g == 0.) {
*cs = 1.;
*sn = 0.;
*r = *f;
} else if (*f == 0.) {
*cs = 0.;
*sn = 1.;
*r = *g;
} else {
f1 = *f;
g1 = *g;
/* Computing MAX */
d__1 = abs(f1), d__2 = abs(g1);
scale = max(d__1,d__2);
if (scale >= safmx2) {
count = 0;
L10:
++count;
f1 *= safmn2;
g1 *= safmn2;
/* Computing MAX */
d__1 = abs(f1), d__2 = abs(g1);
scale = max(d__1,d__2);
if (scale >= safmx2) {
goto L10;
}
/* Computing 2nd power */
d__1 = f1;
/* Computing 2nd power */
d__2 = g1;
*r = sqrt(d__1 * d__1 + d__2 * d__2);
*cs = f1 / *r;
*sn = g1 / *r;
i__1 = count;
for (i = 1; i <= count; ++i) {
*r *= safmx2;
/* L20: */
}
} else if (scale <= safmn2) {
count = 0;
L30:
++count;
f1 *= safmx2;
g1 *= safmx2;
/* Computing MAX */
d__1 = abs(f1), d__2 = abs(g1);
scale = max(d__1,d__2);
if (scale <= safmn2) {
goto L30;
}
/* Computing 2nd power */
d__1 = f1;
/* Computing 2nd power */
d__2 = g1;
*r = sqrt(d__1 * d__1 + d__2 * d__2);
*cs = f1 / *r;
*sn = g1 / *r;
i__1 = count;
for (i = 1; i <= count; ++i) {
*r *= safmn2;
/* L40: */
}
} else {
/* Computing 2nd power */
d__1 = f1;
/* Computing 2nd power */
d__2 = g1;
*r = sqrt(d__1 * d__1 + d__2 * d__2);
*cs = f1 / *r;
*sn = g1 / *r;
}
if (abs(*f) > abs(*g) && *cs < 0.) {
*cs = -(*cs);
*sn = -(*sn);
*r = -(*r);
}
}
return 0;
/* End of DLARTG */
} /* dlartg_ */
|
