jacobi.cpp File Reference

#include <ossim/matrix/include.h>
#include <ossim/matrix/newmatap.h>
#include <ossim/matrix/precisio.h>
#include <ossim/matrix/newmatrm.h>

Go to the source code of this file.

Macros
#define	WANT_MATH
#define	REPORT   {}

Functions
void	Jacobi (const SymmetricMatrix &X, DiagonalMatrix &D, SymmetricMatrix &A, Matrix &V, bool eivec)
void	Jacobi (const SymmetricMatrix &X, DiagonalMatrix &D)
void	Jacobi (const SymmetricMatrix &X, DiagonalMatrix &D, SymmetricMatrix &A)
void	Jacobi (const SymmetricMatrix &X, DiagonalMatrix &D, Matrix &V)

Macro Definition Documentation

REPORT

#define REPORT   {}

Definition at line 23 of file jacobi.cpp.

WANT_MATH

#define WANT_MATH

Definition at line 9 of file jacobi.cpp.

Function Documentation

Jacobi() [1/4]

void Jacobi	(	const SymmetricMatrix &	X,
		DiagonalMatrix &	D,
		SymmetricMatrix &	A,
		Matrix &	V,
		bool	eivec
	)

Definition at line 27 of file jacobi.cpp.

{
   Real epsilon = FloatingPointPrecision::Epsilon();
   Tracer et("Jacobi");
   REPORT
   int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.ReSize(n); A = X;
   if (eivec) { REPORT V.ReSize(n,n); D = 1.0; V = D; }
   B << A; D = B; Z = 0.0; A.Inject(Z);
   bool converged = false;
   for (int i=1; i<=50; i++)
   {
      Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
      while (p--) sm += fabs(*a++);            // have previously zeroed diags
      if (sm==0.0) { REPORT converged = true; break; }
      Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
      for (p = 0; p < n; p++)
      {
         Real* ap1 = a + (p*(p+1))/2;
         Real& zp = Z.element(p); Real& dp = D.element(p);
         for (int q = p+1; q < n; q++)
         {
            Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
            Real& zq = Z.element(q); Real& dq = D.element(q);
            Real& apq = A.element(q,p);
            Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);

            if (i>4 && g < epsilon*adp && g < epsilon*adq) { REPORT apq = 0.0; }
            else if (fabs(apq) > tresh)
            {
               REPORT
               Real t; Real h = dq - dp; Real ah = fabs(h);
               if (g < epsilon*ah) { REPORT t = apq / h; }
               else
               {
                  REPORT
                  Real theta = 0.5 * h / apq;
                  t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
                  if (theta<0.0) { REPORT t = -t; }
               }
               Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
               Real tau = s / (1.0 + c); h = t * apq;
               zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
               int j = p;
               while (j--)
               {
                  g = *ap; h = *aq;
                  *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
               }
               int ip = p+1; j = q-ip; ap += ip++; aq++;
               while (j--)
               {
                  g = *ap; h = *aq;
                  *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
                  ap += ip++;
               }
               if (q < n-1)             // last loop is non-empty
               {
                  int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
                  for (;;)
                  {
                     g = *ap; h = *aq;
                     *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
                     if (!(--j)) break;
                     ap += ip++; aq += iq++;
                  }
               }
               if (eivec)
               {
                  REPORT
                  RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
                  Rotate(VP, VQ, tau, s);
               }
            }
         }
      }
      B = B + Z; D = B; Z = 0.0;
   }
   if (!converged) Throw(ConvergenceException(X));
   if (eivec) SortSV(D, V, true);
   else SortAscending(D);
}

Jacobi() [2/4]

void Jacobi	(	const SymmetricMatrix &	X,
		DiagonalMatrix &	D
	)

Definition at line 110 of file jacobi.cpp.

{ REPORT SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }

Jacobi() [3/4]

void Jacobi	(	const SymmetricMatrix &	X,
		DiagonalMatrix &	D,
		SymmetricMatrix &	A
	)

Definition at line 113 of file jacobi.cpp.

{ REPORT Matrix V; Jacobi(X,D,A,V,false); }

Jacobi() [4/4]

void Jacobi	(	const SymmetricMatrix &	X,
		DiagonalMatrix &	D,
		Matrix &	V
	)

Definition at line 116 of file jacobi.cpp.

{ REPORT SymmetricMatrix A; Jacobi(X,D,A,V,true); }