1 /* ./src_f77/zlartg.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
zlartg_(doublecomplex * f,doublecomplex * g,doublereal * cs,doublecomplex * sn,doublecomplex * r__)8 /* Subroutine */ int zlartg_(doublecomplex *f, doublecomplex *g, doublereal *
9 cs, doublecomplex *sn, doublecomplex *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, d__3, d__4, d__5, d__6, d__7, d__8, d__9, d__10;
18 doublecomplex z__1, z__2, z__3;
19
20 /* Builtin functions */
21 double log(doublereal), pow_di(doublereal *, integer *), d_imag(
22 doublecomplex *), sqrt(doublereal);
23 void d_cnjg(doublecomplex *, doublecomplex *);
24
25 /* Local variables */
26 static doublereal d__;
27 static integer i__;
28 static doublereal f2, g2;
29 static doublecomplex ff;
30 static doublereal di, dr;
31 static doublecomplex fs, gs;
32 static doublereal f2s, g2s, eps, scale;
33 static integer count;
34 static doublereal safmn2;
35 extern doublereal dlapy2_(doublereal *, doublereal *);
36 static doublereal safmx2;
37 extern doublereal dlamch_(char *, ftnlen);
38 static doublereal safmin;
39
40
41 /* -- LAPACK auxiliary routine (version 3.0) -- */
42 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
43 /* Courant Institute, Argonne National Lab, and Rice University */
44 /* June 30, 1999 */
45
46 /* .. Scalar Arguments .. */
47 /* .. */
48
49 /* Purpose */
50 /* ======= */
51
52 /* ZLARTG generates a plane rotation so that */
53
54 /* [ CS SN ] [ F ] [ R ] */
55 /* [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. */
56 /* [ -SN CS ] [ G ] [ 0 ] */
57
58 /* This is a faster version of the BLAS1 routine ZROTG, except for */
59 /* the following differences: */
60 /* F and G are unchanged on return. */
61 /* If G=0, then CS=1 and SN=0. */
62 /* If F=0, then CS=0 and SN is chosen so that R is real. */
63
64 /* Arguments */
65 /* ========= */
66
67 /* F (input) COMPLEX*16 */
68 /* The first component of vector to be rotated. */
69
70 /* G (input) COMPLEX*16 */
71 /* The second component of vector to be rotated. */
72
73 /* CS (output) DOUBLE PRECISION */
74 /* The cosine of the rotation. */
75
76 /* SN (output) COMPLEX*16 */
77 /* The sine of the rotation. */
78
79 /* R (output) COMPLEX*16 */
80 /* The nonzero component of the rotated vector. */
81
82 /* Further Details */
83 /* ======= ======= */
84
85 /* 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel */
86
87 /* ===================================================================== */
88
89 /* .. Parameters .. */
90 /* .. */
91 /* .. Local Scalars .. */
92 /* .. */
93 /* .. External Functions .. */
94 /* .. */
95 /* .. Intrinsic Functions .. */
96 /* .. */
97 /* .. Statement Functions .. */
98 /* .. */
99 /* .. Save statement .. */
100 /* .. */
101 /* .. Data statements .. */
102 /* .. */
103 /* .. Statement Function definitions .. */
104 /* .. */
105 /* .. Executable Statements .. */
106
107 if (first) {
108 first = FALSE_;
109 safmin = dlamch_("S", (ftnlen)1);
110 eps = dlamch_("E", (ftnlen)1);
111 d__1 = dlamch_("B", (ftnlen)1);
112 i__1 = (integer) (log(safmin / eps) / log(dlamch_("B", (ftnlen)1)) /
113 2.);
114 safmn2 = pow_di(&d__1, &i__1);
115 safmx2 = 1. / safmn2;
116 }
117 /* Computing MAX */
118 /* Computing MAX */
119 d__7 = (d__1 = f->r, abs(d__1)), d__8 = (d__2 = d_imag(f), abs(d__2));
120 /* Computing MAX */
121 d__9 = (d__3 = g->r, abs(d__3)), d__10 = (d__4 = d_imag(g), abs(d__4));
122 d__5 = max(d__7,d__8), d__6 = max(d__9,d__10);
123 scale = max(d__5,d__6);
124 fs.r = f->r, fs.i = f->i;
125 gs.r = g->r, gs.i = g->i;
126 count = 0;
127 if (scale >= safmx2) {
128 L10:
129 ++count;
130 z__1.r = safmn2 * fs.r, z__1.i = safmn2 * fs.i;
131 fs.r = z__1.r, fs.i = z__1.i;
132 z__1.r = safmn2 * gs.r, z__1.i = safmn2 * gs.i;
133 gs.r = z__1.r, gs.i = z__1.i;
134 scale *= safmn2;
135 if (scale >= safmx2) {
136 goto L10;
137 }
138 } else if (scale <= safmn2) {
139 if (g->r == 0. && g->i == 0.) {
140 *cs = 1.;
141 sn->r = 0., sn->i = 0.;
142 r__->r = f->r, r__->i = f->i;
143 return 0;
144 }
145 L20:
146 --count;
147 z__1.r = safmx2 * fs.r, z__1.i = safmx2 * fs.i;
148 fs.r = z__1.r, fs.i = z__1.i;
149 z__1.r = safmx2 * gs.r, z__1.i = safmx2 * gs.i;
150 gs.r = z__1.r, gs.i = z__1.i;
151 scale *= safmx2;
152 if (scale <= safmn2) {
153 goto L20;
154 }
155 }
156 /* Computing 2nd power */
157 d__1 = fs.r;
158 /* Computing 2nd power */
159 d__2 = d_imag(&fs);
160 f2 = d__1 * d__1 + d__2 * d__2;
161 /* Computing 2nd power */
162 d__1 = gs.r;
163 /* Computing 2nd power */
164 d__2 = d_imag(&gs);
165 g2 = d__1 * d__1 + d__2 * d__2;
166 if (f2 <= max(g2,1.) * safmin) {
167
168 /* This is a rare case: F is very small. */
169
170 if (f->r == 0. && f->i == 0.) {
171 *cs = 0.;
172 d__2 = g->r;
173 d__3 = d_imag(g);
174 d__1 = dlapy2_(&d__2, &d__3);
175 r__->r = d__1, r__->i = 0.;
176 /* Do complex/real division explicitly with two real divisions */
177 d__1 = gs.r;
178 d__2 = d_imag(&gs);
179 d__ = dlapy2_(&d__1, &d__2);
180 d__1 = gs.r / d__;
181 d__2 = -d_imag(&gs) / d__;
182 z__1.r = d__1, z__1.i = d__2;
183 sn->r = z__1.r, sn->i = z__1.i;
184 return 0;
185 }
186 d__1 = fs.r;
187 d__2 = d_imag(&fs);
188 f2s = dlapy2_(&d__1, &d__2);
189 /* G2 and G2S are accurate */
190 /* G2 is at least SAFMIN, and G2S is at least SAFMN2 */
191 g2s = sqrt(g2);
192 /* Error in CS from underflow in F2S is at most */
193 /* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS */
194 /* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, */
195 /* and so CS .lt. sqrt(SAFMIN) */
196 /* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN */
197 /* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) */
198 /* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S */
199 *cs = f2s / g2s;
200 /* Make sure abs(FF) = 1 */
201 /* Do complex/real division explicitly with 2 real divisions */
202 /* Computing MAX */
203 d__3 = (d__1 = f->r, abs(d__1)), d__4 = (d__2 = d_imag(f), abs(d__2));
204 if (max(d__3,d__4) > 1.) {
205 d__1 = f->r;
206 d__2 = d_imag(f);
207 d__ = dlapy2_(&d__1, &d__2);
208 d__1 = f->r / d__;
209 d__2 = d_imag(f) / d__;
210 z__1.r = d__1, z__1.i = d__2;
211 ff.r = z__1.r, ff.i = z__1.i;
212 } else {
213 dr = safmx2 * f->r;
214 di = safmx2 * d_imag(f);
215 d__ = dlapy2_(&dr, &di);
216 d__1 = dr / d__;
217 d__2 = di / d__;
218 z__1.r = d__1, z__1.i = d__2;
219 ff.r = z__1.r, ff.i = z__1.i;
220 }
221 d__1 = gs.r / g2s;
222 d__2 = -d_imag(&gs) / g2s;
223 z__2.r = d__1, z__2.i = d__2;
224 z__1.r = ff.r * z__2.r - ff.i * z__2.i, z__1.i = ff.r * z__2.i + ff.i
225 * z__2.r;
226 sn->r = z__1.r, sn->i = z__1.i;
227 z__2.r = *cs * f->r, z__2.i = *cs * f->i;
228 z__3.r = sn->r * g->r - sn->i * g->i, z__3.i = sn->r * g->i + sn->i *
229 g->r;
230 z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
231 r__->r = z__1.r, r__->i = z__1.i;
232 } else {
233
234 /* This is the most common case. */
235 /* Neither F2 nor F2/G2 are less than SAFMIN */
236 /* F2S cannot overflow, and it is accurate */
237
238 f2s = sqrt(g2 / f2 + 1.);
239 /* Do the F2S(real)*FS(complex) multiply with two real multiplies */
240 d__1 = f2s * fs.r;
241 d__2 = f2s * d_imag(&fs);
242 z__1.r = d__1, z__1.i = d__2;
243 r__->r = z__1.r, r__->i = z__1.i;
244 *cs = 1. / f2s;
245 d__ = f2 + g2;
246 /* Do complex/real division explicitly with two real divisions */
247 d__1 = r__->r / d__;
248 d__2 = d_imag(r__) / d__;
249 z__1.r = d__1, z__1.i = d__2;
250 sn->r = z__1.r, sn->i = z__1.i;
251 d_cnjg(&z__2, &gs);
252 z__1.r = sn->r * z__2.r - sn->i * z__2.i, z__1.i = sn->r * z__2.i +
253 sn->i * z__2.r;
254 sn->r = z__1.r, sn->i = z__1.i;
255 if (count != 0) {
256 if (count > 0) {
257 i__1 = count;
258 for (i__ = 1; i__ <= i__1; ++i__) {
259 z__1.r = safmx2 * r__->r, z__1.i = safmx2 * r__->i;
260 r__->r = z__1.r, r__->i = z__1.i;
261 /* L30: */
262 }
263 } else {
264 i__1 = -count;
265 for (i__ = 1; i__ <= i__1; ++i__) {
266 z__1.r = safmn2 * r__->r, z__1.i = safmn2 * r__->i;
267 r__->r = z__1.r, r__->i = z__1.i;
268 /* L40: */
269 }
270 }
271 }
272 }
273 return 0;
274
275 /* End of ZLARTG */
276
277 } /* zlartg_ */
278
279