1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162 |
- #include "ABACUS.h"
- using namespace std;
-
- void ABACUS::tred2 (SQMat_DP& a, Vect_DP& d, Vect_DP& e)
- {
- int l, k, j, i;
- DP scale, hh, h, g, f;
-
- int n = a.size();
- for (i = n-1; i > 0; i--) {
- l = i - 1;
- h = scale = 0.0;
- if (l > 0) {
- for (k = 0; k < l + 1; k++) scale += fabs(a[i][k]);
- if (scale == 0.0) e[i] = a[i][l];
- else {
- for (k = 0; k < l + 1; k++) {
- a[i][k] /= scale;
- h += a[i][k] * a[i][k];
- }
- f = a[i][l];
- g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
- e[i] = scale * g;
- h -= f * g;
- a[i][l] = f - g;
- f = 0.0;
- for (j = 0; j < l + 1; j++) {
- a[j][i] = a[i][j]/h;
- g = 0.0;
- for (k = 0; k < j + 1; k++) g += a[j][k] * a[i][k];
- for (k = j + 1; k < l + 1; k++) g += a[k][j] * a[i][k];
- e[j] = g/h;
- f += e[j] * a[i][j];
- }
- hh = f/(h +h);
- for (j = 0; j < l + 1; j++) {
- f = a[i][j];
- e[j] = g = e[j] - hh * f;
- for (k = 0; k < j + 1; k++) a[j][k] -= (f * e[k] + g * a[i][k]);
- }
- }
- } else
- e[i] = a[i][l];
- d[i] = h;
- }
- d[0] = 0.0;
- e[0] = 0.0;
-
- for (i = 0; i < n; i++) {
- l = i;
- if (d[i] != 0.0) {
- for (j = 0; j < l; j++) {
- g = 0.0;
- for (k = 0; k < l; k++) g += a[i][k] * a[k][j];
- for (k = 0; k < l; k++) a[k][j] -= g * a[k][i];
- }
- }
- d[i] = a[i][i];
- a[i][i] = 1.0;
- for (j = 0; j < l; j++) a[j][i] = a[i][j] = 0.0;
- }
- } // void tred2
|