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