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
|
#include "f2c.h"
/* Subroutine */ int zlartg_(doublecomplex *f, doublecomplex *g, doublereal *
cs, doublecomplex *sn, doublecomplex *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
=======
ZLARTG generates a plane rotation so that
[ CS SN ] [ F ] [ R ]
[ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1.
[ -SN CS ] [ G ] [ 0 ]
This is a faster version of the BLAS1 routine ZROTG, except for
the following 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.
Arguments
=========
F (input) COMPLEX*16
The first component of vector to be rotated.
G (input) COMPLEX*16
The second component of vector to be rotated.
CS (output) DOUBLE PRECISION
The cosine of the rotation.
SN (output) COMPLEX*16
The sine of the rotation.
R (output) COMPLEX*16
The nonzero component of the rotated vector.
=====================================================================
[ 25 or 38 ops for main paths ] */
/* System generated locals */
doublereal d__1, d__2;
doublecomplex z__1, z__2, z__3;
/* Builtin functions */
void d_cnjg(doublecomplex *, doublecomplex *);
double z_abs(doublecomplex *), d_imag(doublecomplex *), sqrt(doublereal);
/* Local variables */
static doublereal d, f1, f2, g1, g2, fa, ga, di;
static doublecomplex fs, gs, ss;
if (g->r == 0. && g->i == 0.) {
*cs = 1.;
sn->r = 0., sn->i = 0.;
r->r = f->r, r->i = f->i;
} else if (f->r == 0. && f->i == 0.) {
*cs = 0.;
d_cnjg(&z__2, g);
d__1 = z_abs(g);
z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
sn->r = z__1.r, sn->i = z__1.i;
d__1 = z_abs(g);
r->r = d__1, r->i = 0.;
/* SN = ONE
R = G */
} else {
f1 = (d__1 = f->r, abs(d__1)) + (d__2 = d_imag(f), abs(d__2));
g1 = (d__1 = g->r, abs(d__1)) + (d__2 = d_imag(g), abs(d__2));
if (f1 >= g1) {
z__1.r = g->r / f1, z__1.i = g->i / f1;
gs.r = z__1.r, gs.i = z__1.i;
/* Computing 2nd power */
d__1 = gs.r;
/* Computing 2nd power */
d__2 = d_imag(&gs);
g2 = d__1 * d__1 + d__2 * d__2;
z__1.r = f->r / f1, z__1.i = f->i / f1;
fs.r = z__1.r, fs.i = z__1.i;
/* Computing 2nd power */
d__1 = fs.r;
/* Computing 2nd power */
d__2 = d_imag(&fs);
f2 = d__1 * d__1 + d__2 * d__2;
d = sqrt(g2 / f2 + 1.);
*cs = 1. / d;
d_cnjg(&z__3, &gs);
z__2.r = z__3.r * fs.r - z__3.i * fs.i, z__2.i = z__3.r * fs.i +
z__3.i * fs.r;
d__1 = *cs / f2;
z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
sn->r = z__1.r, sn->i = z__1.i;
z__1.r = d * f->r, z__1.i = d * f->i;
r->r = z__1.r, r->i = z__1.i;
} else {
z__1.r = f->r / g1, z__1.i = f->i / g1;
fs.r = z__1.r, fs.i = z__1.i;
/* Computing 2nd power */
d__1 = fs.r;
/* Computing 2nd power */
d__2 = d_imag(&fs);
f2 = d__1 * d__1 + d__2 * d__2;
fa = sqrt(f2);
z__1.r = g->r / g1, z__1.i = g->i / g1;
gs.r = z__1.r, gs.i = z__1.i;
/* Computing 2nd power */
d__1 = gs.r;
/* Computing 2nd power */
d__2 = d_imag(&gs);
g2 = d__1 * d__1 + d__2 * d__2;
ga = sqrt(g2);
d = sqrt(f2 / g2 + 1.);
di = 1. / d;
*cs = fa / ga * di;
d_cnjg(&z__3, &gs);
z__2.r = z__3.r * fs.r - z__3.i * fs.i, z__2.i = z__3.r * fs.i +
z__3.i * fs.r;
d__1 = fa * ga;
z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
ss.r = z__1.r, ss.i = z__1.i;
z__1.r = di * ss.r, z__1.i = di * ss.i;
sn->r = z__1.r, sn->i = z__1.i;
z__2.r = g->r * ss.r - g->i * ss.i, z__2.i = g->r * ss.i + g->i *
ss.r;
z__1.r = d * z__2.r, z__1.i = d * z__2.i;
r->r = z__1.r, r->i = z__1.i;
}
}
return 0;
/* End of ZLARTG */
} /* zlartg_ */
|
