6 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
9 int G_tqli(
double d[],
double e[],
int n,
double **z)
12 double s,
r, p,
g, f, dd, c,
b;
14 for (i = 1; i < n; i++)
17 for (l = 0; l < n; l++) {
21 for (m = l; m < n - 1; m++) {
22 dd = fabs(d[m]) + fabs(d[m + 1]);
23 if (fabs(e[m]) + dd == dd)
30 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
31 r = sqrt((g * g) + 1.0);
32 g = d[m] - d[l] + e[l] / (g +
SIGN(r, g));
36 for (i = m - 1; i >= l; i--) {
40 if (fabs(f) >= fabs(g)) {
42 r = sqrt((c * c) + 1.0);
48 r = sqrt((s * s) + 1.0);
54 r = (d[i] -
g) * s + 2.0 * c * b;
60 for (k = 0; k < n; k++) {
62 z[k][i + 1] = s * z[k][i] + c * f;
63 z[k][i] = c * z[k][i] - s * f;
77 void G_tred2(
double **a,
int n,
double d[],
double e[])
80 double scale, hh,
h,
g, f;
82 for (i = n - 1; i >= 1; i--) {
87 for (k = 0; k <= l; k++)
88 scale += fabs(a[i][k]);
93 for (k = 0; k <= l; k++) {
95 h += a[i][k] * a[i][k];
99 g = f > 0 ? -sqrt(h) : sqrt(h);
105 for (j = 0; j <= l; j++) {
107 a[j][i] = a[i][j] /
h;
109 for (k = 0; k <= j; k++)
110 g += a[j][k] * a[i][k];
111 for (k = j + 1; k <= l; k++)
112 g += a[k][j] * a[i][k];
118 for (j = 0; j <= l; j++) {
120 e[j] = g = e[j] - hh * f;
122 for (k = 0; k <= j; k++)
123 a[j][k] -= (f * e[k] + g * a[i][k]);
138 for (i = 0; i < n; i++) {
142 for (j = 0; j <= l; j++) {
144 for (k = 0; k <= l; k++)
145 g += a[i][k] * a[k][j];
146 for (k = 0; k <= l; k++)
147 a[k][j] -= g * a[k][i];
153 for (j = 0; j <= l; j++)
154 a[j][i] = a[i][j] = 0.0;