hholder.cpp File Reference

#include <cmath>
#include <ossim/matrix/include.h>
#include <ossim/matrix/newmatap.h>

Go to the source code of this file.

Macros
#define	WANT_MATH
#define	REPORT   {}

Functions
Real	square (Real x)
void	QRZT (Matrix &X, LowerTriangularMatrix &L)
void	QRZT (const Matrix &X, Matrix &Y, Matrix &M)
void	QRZ (Matrix &X, UpperTriangularMatrix &U)
void	QRZ (const Matrix &X, Matrix &Y, Matrix &M)
void	UpdateQRZT (Matrix &X, LowerTriangularMatrix &L)
void	UpdateQRZ (Matrix &X, UpperTriangularMatrix &U)

Macro Definition Documentation

REPORT

#define REPORT   {}

Definition at line 20 of file hholder.cpp.

Referenced by QRZ(), QRZT(), UpdateQRZ(), and UpdateQRZT().

WANT_MATH

#define WANT_MATH

Definition at line 5 of file hholder.cpp.

Function Documentation

QRZ() [1/2]

void QRZ	(	Matrix &	X,
		UpperTriangularMatrix &	U
	)

Definition at line 111 of file hholder.cpp.

References n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, UpperTriangularMatrix::ReSize(), and GeneralMatrix::Store().

Referenced by NonLinearLeastSquares::NextPoint().

{
   REPORT
   Tracer et("QRZ(1)");
   int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s); U = 0.0;
   Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
   int j, k; int J = s; int i = s;
   while (i--)
   {
      Real* xj0 = xi0; Real* xi = xi0; k = n;
      if (k) for (;;)
      {
         u = u0; Real Xi = *xi; Real* xj = xj0;
         j = J; while(j--) *u++ += Xi * *xj++;
         if (!(--k)) break;
         xi += s; xj0 += s;
      }

      Real sum = std::sqrt(*u0); *u0 = sum; u = u0+1;
      if (sum == 0.0)
      {
         REPORT
         j = J - 1; while(j--) *u++ = 0.0;

         xj0 = xi0++; k = n;
         if (k) for (;;)
         {
            *xj0 = 0.0;
            if (!(--k)) break;
              xj0 += s;
         }
         u0 += J--;
      }
      else
      {
         int J1 = J-1; j = J1; while(j--) *u++ /= sum;

         xj0 = xi0; xi = xi0++; k = n;
         if (k) for (;;)
         {
            u = u0+1; Real Xi = *xi; Real* xj = xj0;
            Xi /= sum; *xj++ = Xi;
            j = J1; while(j--) *xj++ -= *u++ * Xi;
            if (!(--k)) break;
              xi += s; xj0 += s;
         }
         u0 += J--;
      }
   }
}

QRZ() [2/2]

void QRZ	(	const Matrix &	X,
		Matrix &	Y,
		Matrix &	M
	)

Definition at line 162 of file hholder.cpp.

References n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, Matrix::ReSize(), and GeneralMatrix::Store().

{
   REPORT
   Tracer et("QRZ(2)");
   int n = X.Nrows(); int s = X.Ncols(); int t = Y.Ncols();
   if (Y.Nrows() != n)
      { Throw(ProgramException("Unequal column lengths",X,Y)); }
   M.ReSize(s,t); M = 0;Real* m0 = M.Store(); Real* m;
   Real* xi0 = X.Store();
   int j, k; int i = s;
   while (i--)
   {
      Real* xj0 = Y.Store(); Real* xi = xi0; k = n;
      if (k) for (;;)
      {
         m = m0; Real Xi = *xi; Real* xj = xj0;
         j = t; while(j--) *m++ += Xi * *xj++;
         if (!(--k)) break;
         xi += s; xj0 += t;
      }

      xj0 = Y.Store(); xi = xi0++; k = n;
      if (k) for (;;)
      {
         m = m0; Real Xi = *xi; Real* xj = xj0;
         j = t; while(j--) *xj++ -= *m++ * Xi;
         if (!(--k)) break;
         xi += s; xj0 += t;
      }
      m0 += t;
   }
}

QRZT() [1/2]

void QRZT	(	Matrix &	X,
		LowerTriangularMatrix &	L
	)

Definition at line 28 of file hholder.cpp.

References LowerTriangularMatrix::element(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, LowerTriangularMatrix::ReSize(), square(), and GeneralMatrix::Store().

Referenced by HHDecompose().

{
   REPORT
     Tracer et("QRZT(1)");
   int n = X.Ncols(); int s = X.Nrows(); L.ReSize(s);
   Real* xi = X.Store(); int k;
   for (int i=0; i<s; i++)
   {
      Real sum = 0.0;
      Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
      sum = std::sqrt(sum);
      if (sum == 0.0)
      {
         REPORT
         k=n; while(k--) { *xi0++ = 0.0; }
         for (int j=i; j<s; j++) L.element(j,i) = 0.0;
      }
      else
      {
         L.element(i,i) = sum;
         Real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
         for (int j=i+1; j<s; j++)
         {
            sum=0.0;
            xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
            xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
            L.element(j,i) = sum;
         }
      }
   }
}

QRZT() [2/2]

void QRZT	(	const Matrix &	X,
		Matrix &	Y,
		Matrix &	M
	)

Definition at line 60 of file hholder.cpp.

References Matrix::element(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, Matrix::ReSize(), and GeneralMatrix::Store().

{
   REPORT
   Tracer et("QRZT(2)");
   int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
   if (Y.Ncols() != n)
      { Throw(ProgramException("Unequal row lengths",X,Y)); }
   M.ReSize(t,s);
   Real* xi = X.Store(); int k;
   for (int i=0; i<s; i++)
   {
      Real* xj0 = Y.Store(); Real* xi0 = xi;
      for (int j=0; j<t; j++)
      {
         Real sum=0.0;
         xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
         xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
         M.element(j,i) = sum;
      }
   }
}

square()

Real square ( Real  x )

inline

Definition at line 26 of file hholder.cpp.

References x.

Referenced by QRZT(), UpdateQRZ(), and UpdateQRZT().

{ return x*x; }

UpdateQRZ()

void UpdateQRZ	(	Matrix &	X,
		UpperTriangularMatrix &	U
	)

Definition at line 261 of file hholder.cpp.

References n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, square(), and GeneralMatrix::Store().

{
   REPORT
   Tracer et("UpdateQRZ");
   int n = X.Nrows(); int s = X.Ncols();
   if (s != U.Ncols())
      Throw(ProgramException("Incompatible dimensions",X,U)); 
   Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
   RowVector V(s); Real* v0 = V.Store(); Real* v; V = 0.0;
   int j, k; int J = s; int i = s;
   while (i--)
   {
      Real* xj0 = xi0; Real* xi = xi0; k = n;
      if (k) for (;;)
      {
         v = v0; Real Xi = *xi; Real* xj = xj0;
         j = J; while(j--) *v++ += Xi * *xj++;
         if (!(--k)) break;
         xi += s; xj0 += s;
      }

      Real r = *u0;
      Real sum = std::sqrt(*v0 + square(r));
      
      if (sum == 0.0)
      {
         REPORT
         u = u0; v = v0;
         j = J; while(j--) { *u++ = 0.0; *v++ = 0.0; }
         xj0 = xi0++; k = n;
         if (k) for (;;)
         {
            *xj0 = 0.0;
            if (!(--k)) break;
              xj0 += s;
         }
         u0 += J--;
      }
      else
      {
         Real frs = std::fabs(r) + sum;
         Real a0 = std::sqrt(frs / sum); Real alpha = a0 / frs;
         if (r <= 0) { REPORT alpha = -alpha; *u0 = sum; }
         else { REPORT *u0 = -sum; }
      
         j = J - 1; v = v0 + 1; u = u0 + 1;     
         while (j--)
            { *v = a0 * *u + alpha * *v; *u -= a0 * *v; ++v; ++u; }

         xj0 = xi0; xi = xi0++; k = n;
         if (k) for (;;)
         {
            v = v0 + 1; Real Xi = *xi; Real* xj = xj0;
            Xi *= alpha; *xj++ = Xi;
            j = J - 1; while(j--) *xj++ -= *v++ * Xi;
            if (!(--k)) break;
              xi += s; xj0 += s;
         }
         
         j = J; v = v0;
         while (j--) *v++ = 0.0;
         
         u0 += J--;
      }
   }
}

UpdateQRZT()

void UpdateQRZT	(	Matrix &	X,
		LowerTriangularMatrix &	L
	)

Definition at line 222 of file hholder.cpp.

References LowerTriangularMatrix::element(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, square(), and GeneralMatrix::Store().

{
   REPORT
     Tracer et("UpdateQRZT");
   int n = X.Ncols(); int s = X.Nrows();
   if (s != L.Nrows())
      Throw(ProgramException("Incompatible dimensions",X,L)); 
   Real* xi = X.Store(); int k;
   for (int i=0; i<s; i++)
   {
      Real r = L.element(i,i); 
      Real sum = 0.0;
      Real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
      sum = std::sqrt(sum + square(r));
      if (sum == 0.0)
      {
         REPORT
         k=n; while(k--) { *xi0++ = 0.0; }
         for (int j=i; j<s; j++) L.element(j,i) = 0.0;
      }
      else
      {
         Real frs = std::fabs(r) + sum;
         Real a0 = std::sqrt(frs / sum); Real alpha = a0 / frs;
         if (r <= 0) { REPORT L.element(i,i) = sum; alpha = -alpha; }
         else { REPORT L.element(i,i) = -sum; }
         Real* xj0=xi0; k=n; while(k--) { *xj0++ *= alpha; }
         for (int j=i+1; j<s; j++)
         {
            sum = 0.0;
            xi=xi0; Real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
            sum += a0 * L.element(j,i);
            xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
            L.element(j,i) -= sum * a0;
         }
      }
   }
}