21 #define REPORT { static ExeCounter ExeCount(__LINE__,17); ++ExeCount; } 31 Tracer et(
"Evalue(tred2)");
39 for (i=n-1; i > 0; i--)
42 int k = i-1;
Real* zik = z + i*n;
43 while (k--) g +=
square(*zik++);
50 Z.
element(i,i-1) = f-g; f = 0.0;
55 *zji = (*zij++)/h; g = 0.0;
56 Real* zjk = z + j*n; zik = z + i*n;
57 k = j;
while (k--) g += *zjk++ * (*zik++);
60 { g += *zjk * (*zik++);
if (!(--k))
break; zjk += n; }
61 *ej++ = g/h; f += g * (*zji); zji += n;
63 Real hh = f / (h + h); zij = z + i*n; ej = E.
Store();
66 f = *zij++; g = *ej - hh * f; *ej++ = g;
67 Real* zjk = z + j*n;
Real* zik = z + i*n;
69 while (k--) *zjk++ -= ( f*(*ek++) + g*(*zik++) );
81 for (
int j=0; j<i; j++)
84 Real* zik = z + i*n;
Real* zkj = z + j;
87 { g += *zik++ * (*zkj);
if (!(--k))
break; zkj += n; }
88 Real* zki = z + i; zkj = z + j;
91 { *zkj -= g * (*zki);
if (!(--k))
break; zkj += n; zki += n; }
94 Real* zij = z + i*n;
Real* zji = z + i;
97 { *zij++ = 0.0; *zji = 0.0;
if (!(--j))
break; zji += n; }
98 D.
element(i) = *zij; *zij = 1.0;
104 Tracer et(
"Evalue(tql2)");
114 Real h = eps * ( fabs(dl) + fabs(el) );
115 if (b < h) {
REPORT b = h; }
117 for (m=l; m<n; m++)
if (fabs(E.
element(m)) <= b)
break;
121 if (m==l) {
REPORT test =
true;
break; }
123 Real g = dl;
Real p = (dl1-g) / (2.0*el);
Real r = sqrt(p*p + 1.0);
124 dl = el / (p < 0.0 ? p-r : p+r);
Real h = g - dl; f += h;
125 Real* dlx = &dl1; i = n-l-1;
while (i--) *dlx++ -= h;
128 for (i=m-1; i>=l; i--)
132 g = c * ei; h = c * p;
133 if ( fabs(p) >= fabs(ei))
136 c = ei / p; r = sqrt(c*c + 1.0);
137 ei1 = s*p*r; s = c/r; c = 1.0/r;
142 c = p / ei; r = sqrt(c*c + 1.0);
143 ei1 = s * ei * r; s = 1.0/r; c /= r;
145 p = c * di - s*g; D.
element(i+1) = h + s * (c*g + s*di);
147 Real* zki = z + i;
Real* zki1 = zki + 1;
int k = n;
151 h = *zki1; *zki1 = s*(*zki) + c*h; *zki = c*(*zki) - s*h;
157 if (fabs(el) <= b) {
REPORT; test =
true;
break; }
186 Tracer et(
"Evalue(tred3)");
192 for (
int i = n-1; i >= 0; i--)
196 while (k--) { f = *a++; *d++ = f; h +=
square(f); }
197 if (h <= tol) {
REPORT *(--ei) = 0.0; h = 0.0; }
201 Real g =
sign(-sqrt(h), f); *(--ei) = g; h -= f*g;
202 f -= g; *(d-1) = f; *(a-1) = f; f = 0.0;
204 for (j = 0; j < i; j++)
208 while (k--) g += *ak++ * *dk++;
210 if (k)
for (;;) { g += *ak * *dk++;
if (!(--k))
break; ak += ++l; }
211 g /= h; *ej++ = g; f += g * *dj++;
215 for (j = 0; j < i; j++)
217 f = *dj++; g = *ej - hh * f; *ej++ = g;
219 while (k--) { *ak++ -= (f * *ek++ + g * *dk++); }
228 Tracer et(
"Evalue(tql1)");
231 int n = D.
Nrows();
int l;
238 Real h = eps * ( fabs(dl) + fabs(el) );
241 for (m=l; m<n; m++)
if (fabs(E.
element(m)) <= b)
break;
245 if (m==l) {
REPORT test =
true;
break; }
247 Real g = dl;
Real p = (dl1-g) / (2.0*el);
Real r = sqrt(p*p + 1.0);
248 dl = el / (p < 0.0 ? p-r : p+r);
Real h = g - dl; f += h;
249 Real* dlx = &dl1; i = n-l-1;
while (i--) *dlx++ -= h;
252 for (i=m-1; i>=l; i--)
256 g = c * ei; h = c * p;
257 if ( fabs(p) >= fabs(ei))
260 c = ei / p; r = sqrt(c*c + 1.0);
261 ei1 = s*p*r; s = c/r; c = 1.0/r;
266 c = p / ei; r = sqrt(c*c + 1.0);
267 ei1 = s * ei * r; s = 1.0/r; c /= r;
269 p = c * di - s*g; D.
element(i+1) = h + s * (c*g + s*di);
272 if (fabs(el) <= b) {
REPORT test =
true;
break; }
280 else {
REPORT test =
true;
break; }
static void tred2(const SymmetricMatrix &A, DiagonalMatrix &D, DiagonalMatrix &E, Matrix &Z)
static void tql1(DiagonalMatrix &D, DiagonalMatrix &E)
void Inject(const GeneralMatrix &GM)
The usual rectangular matrix.
static void tred3(const SymmetricMatrix &X, DiagonalMatrix &D, DiagonalMatrix &E, SymmetricMatrix &A)
FloatVector FloatVector * a
void eigenvalues(const SymmetricMatrix &A, DiagonalMatrix &D, Matrix &Z)
virtual void resize(int, int)
Covergence failure exception.
static void tql2(DiagonalMatrix &D, DiagonalMatrix &E, Matrix &Z)
void SortSV(DiagonalMatrix &D, Matrix &U, bool ascending=false)