9 void qrecvc(
double *ev,
Cpx * evec,
double *dp,
int n)
11 double cc, sc = 0.0, d,
x = 0.0, y, h = 0.0, tzr = 1.e-15;
13 int i, j, k, m, nqr = 50 * n;
17 for (j = 0, m = n - 1; j < nqr; ++j) {
22 if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
25 x = (ev[k] - ev[m]) / 2.;
26 h = sqrt(x * x + dp[k] * dp[k]);
27 if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
29 if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
30 sc = dp[k] / (2. * cc * h);
36 for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
38 p[0].
re = cc * h + sc * p[n].
re;
39 p[n].
re = cc * p[n].
re - sc * h;
41 p[0].
im = cc * h + sc * p[n].
im;
42 p[n].
im = cc * p[n].
im - sc * h;
53 for (k = 0; k < m; ++k) {
56 h = sqrt(x * x + dp[k] * dp[k]);
64 ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
65 for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
67 p[0].
re = cc * h + sc * p[n].
re;
68 p[n].
re = cc * p[n].
re - sc * h;
70 p[0].
im = cc * h + sc * p[n].
im;
71 p[n].
im = cc * p[n].
im - sc * h;
74 ev[k] = ev[k] * cc - y;
75 dp[k - 1] = ev[k] * sc;
76 ev[k] = ev[k] * cc + d;
void qrecvc(double *ev, Cpx *evec, double *dp, int n)