11 void qrecvc(
double *ev,
Cpx *evec,
double *dp,
int n)
13 double cc, sc = 0.0, d,
x = 0.0, y, h = 0.0, tzr = 1.e-15;
15 int i, j, k, m, nqr = 50 * n;
19 for (j = 0, m = n - 1; j < nqr; ++j) {
24 if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
27 x = (ev[k] - ev[m]) / 2.;
28 h = sqrt(
x *
x + dp[k] * dp[k]);
29 if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
31 if ((cc = sqrt((1. +
x / h) / 2.)) != 0.)
32 sc = dp[k] / (2. * cc * h);
38 for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
40 p[0].
re = cc * h + sc * p[n].
re;
41 p[n].
re = cc * p[n].
re - sc * h;
43 p[0].
im = cc * h + sc * p[n].
im;
44 p[n].
im = cc * p[n].
im - sc * h;
55 for (k = 0; k < m; ++k) {
58 h = sqrt(
x *
x + dp[k] * dp[k]);
66 ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
67 for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
69 p[0].
re = cc * h + sc * p[n].
re;
70 p[n].
re = cc * p[n].
re - sc * h;
72 p[0].
im = cc * h + sc * p[n].
im;
73 p[n].
im = cc * p[n].
im - sc * h;
76 ev[k] = ev[k] * cc - y;
77 dp[k - 1] = ev[k] * sc;
78 ev[k] = ev[k] * cc + d;
void qrecvc(double *ev, Cpx *evec, double *dp, int n)