26 #define REPORT { static ExeCounter ExeCount(__LINE__,18); ++ExeCount; } 40 B << A; D = B; Z = 0.0; A.
Inject(Z);
41 bool converged =
false;
42 for (
int i=1; i<=50; i++)
44 Real sm=0.0;
Real*
a = A.Store();
int p = A.Storage();
45 while (p--) sm += fabs(*a++);
46 if (sm==0.0) {
REPORT converged =
true;
break; }
47 Real tresh = (i<4) ? 0.2 * sm /
square(n) : 0.0; a = A.Store();
48 for (p = 0; p < n; p++)
50 Real* ap1 = a + (p*(p+1))/2;
52 for (
int q = p+1; q < n; q++)
54 Real* ap = ap1;
Real* aq = a + (q*(q+1))/2;
56 Real& apq = A.element(q,p);
57 Real g = 100 * fabs(apq);
Real adp = fabs(dp);
Real adq = fabs(dq);
59 if (i>4 && g < epsilon*adp && g < epsilon*adq) {
REPORT apq = 0.0; }
60 else if (fabs(apq) > tresh)
64 if (g < epsilon*ah) {
REPORT t = apq / h; }
68 Real theta = 0.5 * h / apq;
69 t = 1.0 / ( fabs(theta) + sqrt(1.0 +
square(theta)) );
70 if (theta<0.0) {
REPORT t = -t; }
73 Real tau = s / (1.0 + c); h = t * apq;
74 zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
79 *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
81 int ip = p+1; j = q-ip; ap += ip++; aq++;
85 *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
90 int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
94 *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
96 ap += ip++; aq += iq++;
108 B = B + Z; D = B; Z = 0.0;
111 if (eivec)
SortSV(D, V,
true);
void Jacobi(const SymmetricMatrix &X, DiagonalMatrix &D, SymmetricMatrix &A, Matrix &V, bool eivec)
void Rotate(RectMatrixCol &U, RectMatrixCol &V, Real tau, Real s)
void Inject(const GeneralMatrix &GM)
void SortAscending(GeneralMatrix &gm)
The usual rectangular matrix.
FloatVector FloatVector * a
virtual void resize(int, int)
Covergence failure exception.
void SortSV(DiagonalMatrix &D, Matrix &U, bool ascending=false)