30 #define MAX(a, b) ((a)>(b)?(a):(b)) 36 static double hypot2(
double x,
double y) {
42 static void tred2(
double V[n][n],
double d[n],
double e[n]) {
51 for (j = 0; j < n; j++) {
57 for (i = n-1; i > 0; i--) {
63 for (k = 0; k < i; k++) {
64 scale = scale + fabs(d[k]);
68 for (j = 0; j < i; j++) {
77 for (k = 0; k < i; k++) {
89 for (j = 0; j < i; j++) {
95 for (j = 0; j < i; j++) {
98 g = e[j] + V[j][j] * f;
99 for (k = j+1; k <= i-1; k++) {
106 for (j = 0; j < i; j++) {
111 for (j = 0; j < i; j++) {
114 for (j = 0; j < i; j++) {
117 for (k = j; k <= i-1; k++) {
118 V[k][j] -= (f * e[k] + g * d[k]);
129 for (i = 0; i < n-1; i++) {
134 for (k = 0; k <= i; k++) {
135 d[k] = V[k][i+1] / h;
137 for (j = 0; j <= i; j++) {
139 for (k = 0; k <= i; k++) {
140 g += V[k][i+1] * V[k][j];
142 for (k = 0; k <= i; k++) {
147 for (k = 0; k <= i; k++) {
151 for (j = 0; j < n; j++) {
161 static void tql2(
double V[n][n],
double d[n],
double e[n]) {
169 double g,p,r,dl1,h,f,tst1,eps;
170 double c,c2,c3,el1,s,s2;
172 for (i = 1; i < n; i++) {
179 eps = pow(2.0,-52.0);
180 for (l = 0; l < n; l++) {
184 tst1 = MAX(tst1,fabs(d[l]) + fabs(e[l]));
187 if (fabs(e[m]) <= eps*tst1) {
204 p = (d[l+1] - g) / (2.0 * e[l]);
209 d[l] = e[l] / (p + r);
210 d[l+1] = e[l] * (p + r);
213 for (i = l+2; i < n; i++) {
227 for (i = m-1; i >= l; i--) {
237 p = c * d[i] - s * g;
238 d[i+1] = h + s * (c * g + s * d[i]);
242 for (k = 0; k < n; k++) {
244 V[k][i+1] = s * V[k][i] + c * h;
245 V[k][i] = c * V[k][i] - s * h;
248 p = -s * s2 * c3 * el1 * e[l] / dl1;
254 }
while (fabs(e[l]) > eps*tst1);
262 for (i = 0; i < n-1; i++) {
265 for (j = i+1; j < n; j++) {
274 for (j = 0; j < n; j++) {
283 void eigen_decomposition(
double A[n][n],
double V[n][n],
double d[n]) {
286 for (i = 0; i < n; i++) {
287 for (j = 0; j < n; j++) {