newmatap.h File Reference

#include <ossim/matrix/newmat.h>

Go to the source code of this file.

Classes
class	SymmetricEigenAnalysis
class	FFT_Controller
class	MultiRadixCounter

Macros
#define	NEWMATAP_LIB   0

Functions
void	QRZT (Matrix &, LowerTriangularMatrix &)
void	QRZT (const Matrix &, Matrix &, Matrix &)
void	QRZ (Matrix &, UpperTriangularMatrix &)
void	QRZ (const Matrix &, Matrix &, Matrix &)
void	HHDecompose (Matrix &X, LowerTriangularMatrix &L)
void	HHDecompose (const Matrix &X, Matrix &Y, Matrix &M)
void	UpdateQRZT (Matrix &X, LowerTriangularMatrix &L)
void	UpdateQRZ (Matrix &X, UpperTriangularMatrix &U)
ReturnMatrix	Cholesky (const SymmetricMatrix &)
ReturnMatrix	Cholesky (const SymmetricBandMatrix &)
void	UpdateCholesky (UpperTriangularMatrix &chol, RowVector r1Modification)
void	DowndateCholesky (UpperTriangularMatrix &chol, RowVector x)
void	RightCircularUpdateCholesky (UpperTriangularMatrix &chol, int k, int l)
void	LeftCircularUpdateCholesky (UpperTriangularMatrix &chol, int k, int l)
void	SVD (const Matrix &, DiagonalMatrix &, Matrix &, Matrix &, bool=true, bool=true)
void	SVD (const Matrix &, DiagonalMatrix &)
void	SVD (const Matrix &A, DiagonalMatrix &D, Matrix &U, bool withU=true)
void	SortSV (DiagonalMatrix &D, Matrix &U, bool ascending=false)
void	SortSV (DiagonalMatrix &D, Matrix &U, Matrix &V, bool ascending=false)
void	Jacobi (const SymmetricMatrix &, DiagonalMatrix &)
void	Jacobi (const SymmetricMatrix &, DiagonalMatrix &, SymmetricMatrix &)
void	Jacobi (const SymmetricMatrix &, DiagonalMatrix &, Matrix &)
void	Jacobi (const SymmetricMatrix &, DiagonalMatrix &, SymmetricMatrix &, Matrix &, bool=true)
void	EigenValues (const SymmetricMatrix &, DiagonalMatrix &)
void	EigenValues (const SymmetricMatrix &, DiagonalMatrix &, SymmetricMatrix &)
void	EigenValues (const SymmetricMatrix &, DiagonalMatrix &, Matrix &)
void	SortAscending (GeneralMatrix &)
void	SortDescending (GeneralMatrix &)
void	FFT (const ColumnVector &, const ColumnVector &, ColumnVector &, ColumnVector &)
void	FFTI (const ColumnVector &, const ColumnVector &, ColumnVector &, ColumnVector &)
void	RealFFT (const ColumnVector &, ColumnVector &, ColumnVector &)
void	RealFFTI (const ColumnVector &, const ColumnVector &, ColumnVector &)
void	DCT_II (const ColumnVector &, ColumnVector &)
void	DCT_II_inverse (const ColumnVector &, ColumnVector &)
void	DST_II (const ColumnVector &, ColumnVector &)
void	DST_II_inverse (const ColumnVector &, ColumnVector &)
void	DCT (const ColumnVector &, ColumnVector &)
void	DCT_inverse (const ColumnVector &, ColumnVector &)
void	DST (const ColumnVector &, ColumnVector &)
void	DST_inverse (const ColumnVector &, ColumnVector &)
void	FFT2 (const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
void	FFT2I (const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)

Macro Definition Documentation

NEWMATAP_LIB

#define NEWMATAP_LIB   0

Definition at line 6 of file newmatap.h.

Function Documentation

Cholesky() [1/2]

ReturnMatrix Cholesky

(

const SymmetricMatrix &

)

Definition at line 30 of file cholesky.cpp.

References GeneralMatrix::ForReturn(), GeneralMatrix::Nrows(), GeneralMatrix::Release(), REPORT, ossim::square(), and GeneralMatrix::Store().

Referenced by MLE_D_FI::NextPoint().

{
   REPORT
   Tracer trace("Cholesky");
   int nr = S.Nrows();
   LowerTriangularMatrix T(nr);
   Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
   for (int i=0; i<nr; i++)
   {
      Real* tj = t; Real sum; int k;
      for (int j=0; j<i; j++)
      {
         Real* tk = ti; sum = 0.0; k = j;
         while (k--) { sum += *tj++ * *tk++; }
         *tk = (*s++ - sum) / *tj++;
      }
      sum = 0.0; k = i;
      while (k--) { sum += square(*ti++); }
      Real d = *s++ - sum;
      if (d<=0.0)  Throw(NPDException(S));
      *ti++ = std::sqrt(d);
   }
   T.Release(); return T.ForReturn();
}

Cholesky() [2/2]

ReturnMatrix Cholesky

(

const SymmetricBandMatrix &

)

Definition at line 55 of file cholesky.cpp.

References GeneralMatrix::ForReturn(), SymmetricBandMatrix::lower, GeneralMatrix::Nrows(), GeneralMatrix::Release(), REPORT, ossim::square(), and GeneralMatrix::Store().

{
   REPORT
   Tracer trace("Band-Cholesky");
   int nr = S.Nrows(); int m = S.lower;
   LowerBandMatrix T(nr,m);
   Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;

   for (int i=0; i<nr; i++)
   {
      Real* tj = t; Real sum; int l;
      if (i<m) { REPORT l = m-i; s += l; ti += l; l = i; }
      else { REPORT t += (m+1); l = m; }

      for (int j=0; j<l; j++)
      {
         Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
         while (k--) { sum += *tj++ * *tk++; }
         *tk = (*s++ - sum) / *tj++;
      }
      sum = 0.0;
      while (l--) { sum += square(*ti++); }
      Real d = *s++ - sum;
      if (d<=0.0)  Throw(NPDException(S));
      *ti++ = std::sqrt(d);
   }

   T.Release(); return T.ForReturn();
}

DCT()

void DCT	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 393 of file fft.cpp.

References DCT_inverse(), GeneralMatrix::Nrows(), and REPORT.

{
   // Discrete cosine transform, type I
   Tracer trace("DCT");
   REPORT
   DCT_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

DCT_II()

void DCT_II	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 249 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), RealFFT(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

{
   // Discrete cosine transform, type II, of a real series
   Tracer trace("DCT_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = *u++; }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *v = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
   }
}

DCT_II_inverse()

void DCT_II_inverse	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 277 of file fft.cpp.

References n, GeneralMatrix::Nrows(), REPORT, GeneralMatrix::Store(), x, and y.

{
   // Inverse of discrete cosine transform, type II
   Tracer trace("DCT_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *v; *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *(++v); Real wi = *(--w);
      *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = *(--b); }
}

DCT_inverse()

void DCT_inverse	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 359 of file fft.cpp.

References A, ColumnVector::CleanUp(), n, GeneralMatrix::Nrows(), RealFFTI(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

Referenced by DCT().

{
   // Inverse of discrete cosine transform, type I
   Tracer trace("DCT_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *v++; *x++ = vi; *y++ = 0.0;
   Real sum1 = vi / 2.0; Real sum2 = sum1; vi = *v++;
   int i = n2-1;
   while (i--)
   {
      Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
      *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
   }
   sum1 += vi; sum2 -= vi;
   vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real bz = (ai - bi) / 4 / s; Real az = (ai + bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

DowndateCholesky()

void DowndateCholesky	(	UpperTriangularMatrix &	chol,
		RowVector	x
	)

Definition at line 123 of file cholesky.cpp.

References GivensRotation(), GeneralMatrix::Nrows(), pythag(), GeneralMatrix::SumSquare(), BaseMatrix::t(), and x.

{
   int nRC = chol.Nrows();
    
   // solve R^T a = x
   LowerTriangularMatrix L = chol.t();
   ColumnVector a(nRC); a = 0.0;
   int i, j;
    
   for (i = 1; i <= nRC; ++i)
   {
      // accumulate subtr sum
      Real subtrsum = 0.0;
      for(int k = 1; k < i; ++k) subtrsum += a(k) * L(i,k);

      a(i) = (x(i) - subtrsum) / L(i,i);
   }

   // test that l2 norm of a is < 1
   Real squareNormA = a.SumSquare();
   if (squareNormA >= 1.0)
      Throw(ProgramException("DowndateCholesky() fails", chol));

   Real alpha = std::sqrt(1.0 - squareNormA);

   // compute and apply Givens rotations to the vector a
   ColumnVector cGivens(nRC);  cGivens = 0.0;
   ColumnVector sGivens(nRC);  sGivens = 0.0;
   for(i = nRC; i >= 1; i--)
      alpha = pythag(alpha, a(i), cGivens(i), sGivens(i));

   // apply Givens rotations to the jth column of chol
   ColumnVector xtilde(nRC); xtilde = 0.0;
   for(j = nRC; j >= 1; j--)
   {
      // only the first j rotations have an affect on chol,0
      for(int k = j; k >= 1; k--)
         GivensRotation(cGivens(k), -sGivens(k), chol(k,j), xtilde(j));
   }
}

DST()

void DST	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 430 of file fft.cpp.

References DST_inverse(), GeneralMatrix::Nrows(), and REPORT.

{
   // Discrete sine transform, type I
   Tracer trace("DST");
   REPORT
   DST_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

DST_II()

void DST_II	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 304 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), RealFFT(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

{
   // Discrete sine transform, type II, of a real series
   Tracer trace("DST_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = -(*u++); }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *(--w) = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
   }
}

DST_II_inverse()

void DST_II_inverse	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 332 of file fft.cpp.

References n, GeneralMatrix::Nrows(), REPORT, GeneralMatrix::Store(), x, and y.

{
   // Inverse of discrete sine transform, type II
   Tracer trace("DST_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *(--w); *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *v++; Real wi = *(--w);
      *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
}

DST_inverse()

void DST_inverse	(	const ColumnVector &	,
		ColumnVector &
	)

Definition at line 402 of file fft.cpp.

References A, ColumnVector::CleanUp(), n, GeneralMatrix::Nrows(), RealFFTI(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

Referenced by DST().

{
   // Inverse of discrete sine transform, type I
   Tracer trace("DST_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
   int i = n2-1;
   while (i--) { *y++ = *(++v); Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
   *x = -2 * vi; *y = 0.0;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = 0.0; *v-- = 0.0;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real az = (ai + bi) / 4 / s; Real bz = (ai - bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

EigenValues() [1/3]

void EigenValues	(	const SymmetricMatrix &	,
		DiagonalMatrix &
	)

Definition at line 286 of file evalue.cpp.

Referenced by ossimMatrix3x3::getEigenValues(), and ossimMatrix4x4::getEigenValues().

{ REPORT DiagonalMatrix E; SymmetricMatrix A; tred3(X,D,E,A); tql1(D,E); }

EigenValues() [2/3]

void EigenValues	(	const SymmetricMatrix &	,
		DiagonalMatrix &	,
		SymmetricMatrix &
	)

Definition at line 289 of file evalue.cpp.

{ REPORT DiagonalMatrix E; tred3(X,D,E,A); tql1(D,E); }

EigenValues() [3/3]

void EigenValues	(	const SymmetricMatrix &	,
		DiagonalMatrix &	,
		Matrix &
	)

Definition at line 283 of file evalue.cpp.

{ REPORT DiagonalMatrix E; tred2(A, D, E, Z); tql2(D, E, Z); SortSV(D,Z,true); }

FFT()

void FFT	(	const ColumnVector &	,
		const ColumnVector &	,
		ColumnVector &	,
		ColumnVector &
	)

Definition at line 197 of file fft.cpp.

References FFT_Controller::ar_1d_ft(), FFT_Controller::CanFactor(), n, GeneralMatrix::Nrows(), FFT_Controller::OnlyOldFFT, REPORT, and GeneralMatrix::Store().

Referenced by FFT2(), and FFTI().

{
   // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
   // but first try Sande and Gentleman
   Tracer trace("FFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   if (n != V.Nrows() || n == 0)
      Throw(ProgramException("Vector lengths unequal or zero", U, V));
   if (n == 1) { REPORT X = U; Y = V; return; }

   // see if we can use the newfft routine
   if (!FFT_Controller::OnlyOldFFT && FFT_Controller::CanFactor(n))
   {
      REPORT
      X = U; Y = V;
      if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
   }

   ColumnVector B = V;
   ColumnVector A = U;
   X.ReSize(n); Y.ReSize(n);
   const int nextmx = 8;
   int prime[8] = { 2,3,5,7,11,13,17,19 };
   int after = 1; int before = n; int next = 0; bool inzee = true;
   int now = 0; int b1;             // initialised to keep gnu happy

   do
   {
      for (;;)
      {
     if (next < nextmx) { REPORT now = prime[next]; }
     b1 = before / now;  if (b1 * now == before) { REPORT break; }
     next++; now += 2;
      }
      before = b1;

      if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
      else { REPORT fftstep(X, Y, A, B, after, now, before); }

      inzee = !inzee; after *= now;
   }
   while (before != 1);

   if (inzee) { REPORT A.Release(); X = A; B.Release(); Y = B; }
}

FFT2()

void FFT2	(	const Matrix &	U,
		const Matrix &	V,
		Matrix &	X,
		Matrix &	Y
	)

Definition at line 440 of file fft.cpp.

References BaseMatrix::Column(), FFT(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, BaseMatrix::Row(), and BaseMatrix::t().

Referenced by FFT2I(), and ossimFftFilter::runFft().

{
   Tracer trace("FFT2");
   REPORT
   int m = U.Nrows(); int n = U.Ncols();
   if (m != V.Nrows() || n != V.Ncols() || m == 0 || n == 0)
      Throw(ProgramException("Matrix dimensions unequal or zero", U, V));
   X = U; Y = V;
   int i; ColumnVector CVR; ColumnVector CVI;
   for (i = 1; i <= m; ++i)
   {
      FFT(X.Row(i).t(), Y.Row(i).t(), CVR, CVI);
      X.Row(i) = CVR.t(); Y.Row(i) = CVI.t();
   }
   for (i = 1; i <= n; ++i)
   {
      FFT(X.Column(i), Y.Column(i), CVR, CVI);
      X.Column(i) = CVR; Y.Column(i) = CVI;
   }
}

FFT2I()

void FFT2I	(	const Matrix &	U,
		const Matrix &	V,
		Matrix &	X,
		Matrix &	Y
	)

Definition at line 461 of file fft.cpp.

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

Referenced by ossimFftFilter::runFft().

{
   // Inverse transform
   Tracer trace("FFT2I");
   REPORT
   FFT2(U,-V,X,Y);
   const Real n = X.Nrows() * X.Ncols(); X /= n; Y /= (-n);
}

FFTI()

void FFTI	(	const ColumnVector &	,
		const ColumnVector &	,
		ColumnVector &	,
		ColumnVector &
	)

Definition at line 116 of file fft.cpp.

References FFT(), n, GeneralMatrix::Nrows(), and REPORT.

{
   // Inverse transform
   Tracer trace("FFTI");
   REPORT
   FFT(U,-V,X,Y);
   const Real n = X.Nrows(); X /= n; Y /= (-n);
}

HHDecompose() [1/2]

void HHDecompose ( Matrix &  X, LowerTriangularMatrix &  L  )

inline

Definition at line 26 of file newmatap.h.

References QRZT().

{ QRZT(X,L); }

HHDecompose() [2/2]

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

inline

Definition at line 29 of file newmatap.h.

References QRZT().

{ QRZT(X, Y, M); }

Jacobi() [1/4]

void Jacobi	(	const SymmetricMatrix &	,
		DiagonalMatrix &
	)

Definition at line 110 of file jacobi.cpp.

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

Jacobi() [2/4]

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

Definition at line 113 of file jacobi.cpp.

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

Jacobi() [3/4]

void Jacobi	(	const SymmetricMatrix &	,
		DiagonalMatrix &	,
		Matrix &
	)

Definition at line 116 of file jacobi.cpp.

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

Jacobi() [4/4]

void Jacobi	(	const SymmetricMatrix &	,
		DiagonalMatrix &	,
		SymmetricMatrix &	,
		Matrix &	,
		bool	= true
	)

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);
}

LeftCircularUpdateCholesky()

void LeftCircularUpdateCholesky	(	UpperTriangularMatrix &	chol,
		int	k,
		int	l
	)

Definition at line 225 of file cholesky.cpp.

References BaseMatrix::Column(), GivensRotationR(), GeneralMatrix::Nrows(), and pythag().

{
   int nRC = chol.Nrows();
   int i, j;

   // I. compute shift of column k to the lth position
   Matrix cholCopy = chol;
   // a. grab column k
   ColumnVector columnK = cholCopy.Column(k);
   // b. shift columns k+1,...l to the LEFT
   for(j = k+1; j <= l; ++j)
      cholCopy.Column(j-1) = cholCopy.Column(j);
   // c. copy the elements of columnK into the lth column of cholCopy
   cholCopy.Column(l) = 0.0;
   for(i = 1; i <= k; ++i)
      cholCopy(i,l) = columnK(i);

   // II. apply and compute Given's rotations
   int nGivens = l-k;
   ColumnVector cGivens(nGivens); cGivens = 0.0;
   ColumnVector sGivens(nGivens); sGivens = 0.0;
   for(j = k; j <= nRC; ++j)
   {
      ColumnVector columnJ = cholCopy.Column(j);

      // apply the previous Givens rotations to columnJ
      int imax = j - k; if (imax > nGivens) imax = nGivens;
      for(int i = 1; i <= imax; ++i)
      {
         int gIndex = i;
         int topRowIndex = k + i - 1;
         GivensRotationR(cGivens(gIndex), sGivens(gIndex),
            columnJ(topRowIndex), columnJ(topRowIndex+1));
      }

      // compute a new Given's rotation when j < l
      if(j < l)
      {
         int gIndex = j-k+1;
         columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex), sGivens(gIndex));
         columnJ(j+1) = 0.0;
      }

      cholCopy.Column(j) = columnJ;
   }

   chol << cholCopy;
    
}

QRZ() [1/2]

void QRZ	(	Matrix &	,
		UpperTriangularMatrix &
	)

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 &	,
		Matrix &	,
		Matrix &
	)

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 &	,
		LowerTriangularMatrix &
	)

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 &	,
		Matrix &	,
		Matrix &
	)

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;
      }
   }
}

RealFFT()

void RealFFT	(	const ColumnVector &	,
		ColumnVector &	,
		ColumnVector &
	)

Definition at line 126 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), and REPORT.

Referenced by DCT_II(), and DST_II().

{
   // Fourier transform of a real series
   Tracer trace("RealFFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n2), B(n2);
   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
   while (i--) { *a++ = *u++; *b++ = *u++; }
   FFT(A,B,A,B);
   int n21 = n2 + 1;
   X.ReSize(n21); Y.ReSize(n21);
   i = n2 - 1;
   a = A.Store(); b = B.Store();              // first els of A and B
   Real* an = a + i; Real* bn = b + i;        // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y

   *x++ = *a + *b; *y++ = 0.0;                // first complex element
   *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element

   int j = -1; i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
      *x++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);
      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);
   }
}

RealFFTI()

void RealFFTI	(	const ColumnVector &	,
		const ColumnVector &	,
		ColumnVector &
	)

Definition at line 162 of file fft.cpp.

Referenced by DCT_inverse(), and DST_inverse().

{
   // inverse of a Fourier transform of a real series
   Tracer trace("RealFFTI");
   REPORT
   const int n21 = A.Nrows();                     // length of arrays
   if (n21 != B.Nrows() || n21 == 0)
      Throw(ProgramException("Vector lengths unequal or zero", A, B));
   const int n2 = n21 - 1;  const int n = 2 * n2;  int i = n2 - 1;

   ColumnVector X(n2), Y(n2);
   Real* a = A.Store(); Real* b = B.Store();  // first els of A and B
   Real* an = a + n2;   Real* bn = b + n2;    // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + i;    Real* yn = y + i;     // last els of X and Y

   Real hn = 0.5 / n2;
   *x++  = hn * (*a + *an);  *y++  = - hn * (*a - *an);
   a++; an--; b++; bn--;
   int j = -1;  i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
      *x++  =  hn * ( ap + samcbp); *y++  =  - hn * ( bm + sbpcam);
      *xn-- =  hn * ( ap - samcbp); *yn-- =  - hn * (-bm + sbpcam);
   }
   FFT(X,Y,X,Y);             // have done inverting elsewhere
   U.ReSize(n); i = n2;
   x = X.Store(); y = Y.Store(); Real* u = U.Store();
   while (i--) { *u++ = *x++; *u++ = - *y++; }
}

RightCircularUpdateCholesky()

void RightCircularUpdateCholesky	(	UpperTriangularMatrix &	chol,
		int	k,
		int	l
	)

Definition at line 170 of file cholesky.cpp.

References BaseMatrix::Column(), GivensRotationR(), GeneralMatrix::Nrows(), and pythag().

{
   int nRC = chol.Nrows();
   int i, j;
    
   // I. compute shift of column l to the kth position
   Matrix cholCopy = chol;
   // a. grab column l
   ColumnVector columnL = cholCopy.Column(l);
   // b. shift columns k,...l-1 to the RIGHT
   for(j = l-1; j >= k; --j)
      cholCopy.Column(j+1) = cholCopy.Column(j);
   // c. copy the top k-1 elements of columnL into the kth column of cholCopy
   cholCopy.Column(k) = 0.0;
   for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);

    // II. determine the l-k Given's rotations
   int nGivens = l-k;
   ColumnVector cGivens(nGivens); cGivens = 0.0;
   ColumnVector sGivens(nGivens); sGivens = 0.0;
   for(i = l; i > k; i--)
   {
      int givensIndex = l-i+1;
      columnL(i-1) = pythag(columnL(i-1), columnL(i),
         cGivens(givensIndex), sGivens(givensIndex));
      columnL(i) = 0.0;
   }
   // the kth entry of columnL is the new diagonal element in column k of cholCopy
   cholCopy(k,k) = columnL(k);
    
   // III. apply these Given's rotations to subsequent columns
   // for columns k+1,...,l-1 we only need to apply the last nGivens-(j-k) rotations
   for(j = k+1; j <= nRC; ++j)
   {
      ColumnVector columnJ = cholCopy.Column(j);
      int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
      for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
      {
         // apply gIndex Given's rotation
         int topRowIndex = k + nGivens - gIndex;
         GivensRotationR(cGivens(gIndex), sGivens(gIndex),
            columnJ(topRowIndex), columnJ(topRowIndex+1));
      }
      cholCopy.Column(j) = columnJ;
   }

   chol << cholCopy;
}

SortAscending()

void SortAscending

(

GeneralMatrix &

)

Definition at line 121 of file sort.cpp.

References DoSimpleSort, max, REPORT, GeneralMatrix::Storage(), and GeneralMatrix::Store().

{
   REPORT
   Tracer et("QuickSortAscending");

   Real* data = GM.Store(); int max = GM.Storage();

   if (max > DoSimpleSort) MyQuickSortAscending(data, data + max - 1, 0);
   InsertionSortAscending(data, max, DoSimpleSort);

}

SortDescending()

void SortDescending

(

GeneralMatrix &

)

Definition at line 45 of file sort.cpp.

References DoSimpleSort, max, REPORT, GeneralMatrix::Storage(), and GeneralMatrix::Store().

Referenced by SVD().

{
   REPORT
   Tracer et("QuickSortDescending");

   Real* data = GM.Store(); int max = GM.Storage();

   if (max > DoSimpleSort) MyQuickSortDescending(data, data + max - 1, 0);
   InsertionSortDescending(data, max, DoSimpleSort);

}

SortSV() [1/2]

void SortSV	(	DiagonalMatrix &	D,
		Matrix &	U,
		bool	ascending = false
	)

Definition at line 190 of file sort.cpp.

Referenced by SVD().

{
   REPORT
   Tracer trace("SortSV_DU");
   int m = U.Nrows(); int n = U.Ncols();
   if (n != D.Nrows()) Throw(IncompatibleDimensionsException(D,U));
   Real* u = U.Store();
   for (int i=0; i<n; i++)
   {
      int k = i; Real p = D.element(i);
      if (ascending)
      {
         for (int j=i+1; j<n; j++)
            { if (D.element(j) < p) { k = j; p = D.element(j); } }
      }
      else
      {
         for (int j=i+1; j<n; j++)
         { if (D.element(j) > p) { k = j; p = D.element(j); } }
      }
      if (k != i)
      {
         D.element(k) = D.element(i); D.element(i) = p; int j = m;
         Real* uji = u + i; Real* ujk = u + k;
         if (j) for(;;)
         {
            p = *uji; *uji = *ujk; *ujk = p;
            if (!(--j)) break;
            uji += n; ujk += n;
         }
      }
   }
}

SortSV() [2/2]

void SortSV	(	DiagonalMatrix &	D,
		Matrix &	U,
		Matrix &	V,
		bool	ascending = false
	)

Definition at line 224 of file sort.cpp.

{
   REPORT
   Tracer trace("SortSV_DUV");
   int mu = U.Nrows(); int mv = V.Nrows(); int n = D.Nrows();
   if (n != U.Ncols()) Throw(IncompatibleDimensionsException(D,U));
   if (n != V.Ncols()) Throw(IncompatibleDimensionsException(D,V));
   Real* u = U.Store(); Real* v = V.Store();
   for (int i=0; i<n; i++)
   {
      int k = i; Real p = D.element(i);
      if (ascending)
      {
         for (int j=i+1; j<n; j++)
            { if (D.element(j) < p) { k = j; p = D.element(j); } }
      }
      else
      {
         for (int j=i+1; j<n; j++)
         { if (D.element(j) > p) { k = j; p = D.element(j); } }
      }
      if (k != i)
      {
         D.element(k) = D.element(i); D.element(i) = p;
         Real* uji = u + i; Real* ujk = u + k;
         Real* vji = v + i; Real* vjk = v + k;
         int j = mu;
         if (j) for(;;)
         {
            p = *uji; *uji = *ujk; *ujk = p; if (!(--j)) break;
            uji += n; ujk += n;
         }
         j = mv;
         if (j) for(;;)
         {
            p = *vji; *vji = *vjk; *vjk = p; if (!(--j)) break;
            vji += n; vjk += n;
         }
      }
   }
}

SVD() [1/3]

void SVD	(	const Matrix &	,
		DiagonalMatrix &	,
		Matrix &	,
		Matrix &	,
		bool	= true,
		bool	= true
	)

Definition at line 26 of file svd.cpp.

References A, RectMatrixRowCol::AddScaled(), ComplexScale(), RectMatrixRowCol::Divide(), RectMatrixRow::Down(), RectMatrixCol::Down(), RectMatrixRowCol::DownDiag(), DiagonalMatrix::element(), RowVector::element(), RectMatrixRowCol::First(), RectMatrixCol::Left(), Maximum(), Minimum(), n, RectMatrixRowCol::Negate(), pythag(), REPORT, RectMatrixRow::Reset(), Matrix::ReSize(), DiagonalMatrix::ReSize(), RectMatrixRow::Right(), RectMatrixCol::Right(), sign(), SortDescending(), SortSV(), ossim::square(), RectMatrixRowCol::SumSquare(), RectMatrixCol::Up(), RectMatrixRowCol::UpDiag(), x, y, and RectMatrixRowCol::Zero().

Referenced by ossimRpcSolver::invert(), ossimSensorModelTuple::invert(), ossimRpcProjection::invert(), ossimSensorModel::invert(), and ossimPolynom< ossim_float64, 3 >::invert().

{
   REPORT
   Tracer trace("SVD");
   Real eps = FloatingPointPrecision::Epsilon();
   Real tol = FloatingPointPrecision::Minimum()/eps;

   int m = A.Nrows(); int n = A.Ncols();
   if (m<n)
      Throw(ProgramException("Want no. Rows >= no. Cols", A));
   if (withV && &U == &V)
      Throw(ProgramException("Need different matrices for U and V", U, V));
   U = A; Real g = 0.0; Real f,h; Real x = 0.0; int i;
   RowVector E(n); RectMatrixRow EI(E,0); Q.ReSize(n);
   RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1);

   if (n) for (i=0;;)
   {
      EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare();
      if (s<tol) { REPORT Q.element(i) = 0.0; }
      else
      {
         REPORT
         f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g;
         Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i;
         while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); }
      }

      s = URI.SumSquare();
      if (s<tol) { REPORT g = 0.0; }
      else
      {
         REPORT
         f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g;
         EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i;
         while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); }
      }

      Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) { REPORT x = y; }
      if (++i == n) { REPORT break; }
      UCI.DownDiag(); URI.DownDiag();
   }

   if (withV)
   {
      REPORT
      V.ReSize(n,n); V = 0.0; RectMatrixCol VCI(V,n-1,n-1,1);
      if (n) { VCI.First() = 1.0; g=E.element(n-1); if (n!=1) URI.UpDiag(); }
      for (i=n-2; i>=0; i--)
      {
         VCI.Left();
         if (g!=0.0)
         {
            VCI.Divide(URI, URI.First()*g); int j = n-i;
            RectMatrixCol VCJ = VCI;
            while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); }
         }
         VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i);
         if (i==0) break;
         URI.UpDiag();
      }
   }

   if (withU)
   {
      REPORT
      for (i=n-1; i>=0; i--)
      {
         g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero();
         if (g!=0.0)
         {
            h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI;
            while (--j)
            {
               UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ;
               UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h);
            }
            UCI.Divide(g);
         }
         else UCI.Zero();
         UCI.First() += 1.0;
         if (i==0) break;
         UCI.UpDiag();
      }
   }

   eps *= x;
   for (int k=n-1; k>=0; k--)
   {
      Real z = -FloatingPointPrecision::Maximum(); // to keep Gnu happy
      Real y; int limit = 50; int l = 0;
      while (limit--)
      {
         Real c, s; int i; int l1=k; bool tfc=false;
         for (l=k; l>=0; l--)
         {
//          if (fabs(E.element(l))<=eps) goto test_f_convergence;
            if (fabs(E.element(l))<=eps) { REPORT tfc=true; break; }
            if (fabs(Q.element(l-1))<=eps) { REPORT l1=l; break; }
            REPORT
         }
         if (!tfc)
         {
            REPORT
            l=l1; l1=l-1; s = -1.0; c = 0.0;
            for (i=l; i<=k; i++)
            {
               f = - s * E.element(i); E.element(i) *= c;
//             if (fabs(f)<=eps) goto test_f_convergence;
               if (fabs(f)<=eps) { REPORT break; }
               g = Q.element(i); h = pythag(g,f,c,s); Q.element(i) = h;
               if (withU)
               {
                  REPORT
                  RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1);
                  ComplexScale(UCJ, UCI, c, s);
               }
            }
         }
//       test_f_convergence: z = Q.element(k); if (l==k) goto convergence;
         z = Q.element(k);  if (l==k) { REPORT break; }

         x = Q.element(l); y = Q.element(k-1);
         g = E.element(k-1); h = E.element(k);
         f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y);
         if (f>1)         { REPORT g = f * sqrt(1 + square(1/f)); }
         else if (f<-1)   { REPORT g = -f * sqrt(1 + square(1/f)); }
         else             { REPORT g = sqrt(f*f + 1); }
            { REPORT f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x; }

         c = 1.0; s = 1.0;
         for (i=l+1; i<=k; i++)
         {
            g = E.element(i); y = Q.element(i); h = s*g; g *= c;
            z = pythag(f,h,c,s); E.element(i-1) = z;
            f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c;
            if (withV)
            {
               REPORT
               RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1);
               ComplexScale(VCI, VCJ, c, s);
            }
            z = pythag(f,h,c,s); Q.element(i-1) = z;
            f = c*g + s*y; x = -s*g + c*y;
            if (withU)
            {
               REPORT
               RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1);
               ComplexScale(UCI, UCJ, c, s);
            }
         }
         E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x;
      }
      if (l!=k) { Throw(ConvergenceException(A)); }
// convergence:
      if (z < 0.0)
      {
         REPORT
         Q.element(k) = -z;
         if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); }
      }
   }
   if (withU & withV) SortSV(Q, U, V);
   else if (withU) SortSV(Q, U);
   else if (withV) SortSV(Q, V);
   else SortDescending(Q);
}

SVD() [2/3]

void SVD	(	const Matrix &	,
		DiagonalMatrix &
	)

Definition at line 196 of file svd.cpp.

{ REPORT Matrix U; SVD(A, D, U, U, false, false); }

SVD() [3/3]

void SVD ( const Matrix &  A, DiagonalMatrix &  D, Matrix &  U, bool  withU = true  )

inline

Definition at line 66 of file newmatap.h.

                      { SVD(A, D, U, U, withU, false); }

UpdateCholesky()

void UpdateCholesky	(	UpperTriangularMatrix &	chol,
		RowVector	r1Modification
	)

Definition at line 96 of file cholesky.cpp.

References GivensRotation(), GeneralMatrix::Nrows(), and pythag().

{
   int ncols = chol.Nrows();
   ColumnVector cGivens(ncols); cGivens = 0.0;
   ColumnVector sGivens(ncols); sGivens = 0.0;
    
   for(int j = 1; j <= ncols; ++j) // process the jth column of chol
   {
      // apply the previous Givens rotations k = 1,...,j-1 to column j
      for(int k = 1; k < j; ++k)
         GivensRotation(cGivens(k), sGivens(k), chol(k,j), r1Modification(j));

      // determine the jth Given's rotation
      pythag(chol(j,j), r1Modification(j), cGivens(j), sGivens(j));

      // apply the jth Given's rotation
      {
         Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * r1Modification(j);
         chol(j,j) = tmp0; r1Modification(j) = 0.0;
      }

   }

}

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;
         }
      }
   }
}