tmtd.cpp

非常好用的用C编写的矩阵类,可在不同编译器下编译使用.

/// \ingroup newmat
///@{

/// \file tmtd.cpp
/// Part of matrix library test program.


//#define WANT_STREAM

#include "include.h"
#include "newmatap.h"

#include "tmt.h"

#ifdef use_namespace
using namespace NEWMAT;
#endif

ReturnMatrix Inverter1(const CroutMatrix& X)
{
   Matrix Y = X.i();
   Y.Release();
   return Y.ForReturn();
}

// this version forces a copy
ReturnMatrix Inverter2(CroutMatrix X)
{
   Matrix Y = X.i();
   Y.Release();
   return Y.ForReturn();
}

ReturnMatrix Inverter1(const BandLUMatrix& X)
{
   Matrix Y = X.i();
   Y.Release();
   return Y.ForReturn();
}

// this version forces a copy
ReturnMatrix Inverter2(BandLUMatrix X)
{
   Matrix Y = X.i();
   Y.Release();
   return Y.ForReturn();
}

ReturnMatrix LU1(const Matrix& A)
{
   Tracer et1("LU1 - Crout");
   CroutMatrix X = A;
   return X.for_return();
}

ReturnMatrix LU2(const Matrix& A)
{
   Tracer et1("LU2 - Crout");
   CroutMatrix X = A; X.release();
   return X.for_return();
}

ReturnMatrix LU3(const Matrix& A)
{
   Tracer et1("LU3 - Crout");
   CroutMatrix* X = new CroutMatrix(A); MatrixErrorNoSpace(X);
   X->release_and_delete();
   return X->for_return();
}

ReturnMatrix LU1(const BandMatrix& A)
{
   Tracer et1("LU1 - BandLU");
   BandLUMatrix X = A;
   return X.for_return();
}

ReturnMatrix LU2(const BandMatrix& A)
{
   Tracer et1("LU2 - BandLU");
   BandLUMatrix X = A; X.release();
   return X.for_return();
}

ReturnMatrix LU3(const BandMatrix& A)
{
   Tracer et1("LU3 - BandLU");
   BandLUMatrix* X = new BandLUMatrix(A); MatrixErrorNoSpace(X);
   X->release_and_delete();
   return X->for_return();
}



void CircularShift(const Matrix& X1, int first, int last)
{
      Matrix X; UpperTriangularMatrix U1, U2;
      int n = X1.Ncols();

      // Try right circular shift of columns
      X = X1; QRZ(X, U1);
      RightCircularUpdateCholesky(U1, first, last);
      X = X1.Columns(1,first-1) | X1.Column(last)
         | X1.Columns(first,last-1) | X1.Columns(last+1,n);
      QRZ(X, U2);
      X = U1 - U2; Clean(X, 0.000000001); Print(X);

      // Try left circular shift of columns
      X = X1; QRZ(X, U1);
      LeftCircularUpdateCholesky(U1, first, last);
      X = X1.Columns(1,first-1) | X1.Columns(first+1,last)
         | X1.Column(first) | X1.Columns(last+1,n);
      QRZ(X, U2);
      X = U1 - U2; Clean(X, 0.000000001); Print(X);
}

class TestUpdateQRZ
{
   int m,n1,n2,n3;
   Matrix X1, X2, X3;
   MultWithCarry mwc;   // Uniform random number generator
public:
   void Reset();
   TestUpdateQRZ(int mx, int n1x, int n2x=0, int n3x=0)
      : m(mx), n1(n1x), n2(n2x), n3(n3x) { Reset(); }
   void DoTest();
   void ClearRow(int i)       { X1.Row(i) = 0.0; }
   void SetRow(int i, int j)  { X1.Row(i) = X1.Row(j); }
};



void trymatd()
{
   Tracer et("Thirteenth test of Matrix package");
   Tracer::PrintTrace();
   Matrix X(5,20);
   int i,j;
   for (j=1;j<=20;j++) X(1,j) = j+1;
   for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
   SymmetricMatrix S; S << X * X.t();
   Matrix SM = X * X.t() - S;
   Print(SM);
   LowerTriangularMatrix L = Cholesky(S);
   Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
   Print(Diff);
   {
      Tracer et1("Stage 1");
      LowerTriangularMatrix L1(5);
      Matrix Xt = X.t(); Matrix Xt2 = Xt;
      QRZT(X,L1);
      Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
      UpperTriangularMatrix Ut(5);
      QRZ(Xt,Ut);
      Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
      Matrix Y(3,20);
      for (j=1;j<=20;j++) Y(1,j) = 22-j;
      for (i=2;i<=3;i++) for (j=1;j<=20; j++)
         Y(i,j) = (long)Y(i-1,j) * j % 101;
      Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
      QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
      Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Y2 * Xt2 * S.i() - M * L.i();
      Clean(Diff,0.000000001); Print(Diff);
   }

   ColumnVector C1(5);
   {
      Tracer et1("Stage 2");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A = X;
      ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }

   {
      Tracer et1("Stage 3");
      X.ReSize(7,7);
      for (j=1;j<=7;j++) X(1,j) = j+1;
      for (i=2;i<=7;i++) for (j=1;j<=7; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(7);
      for (i=1;i<=7;i++) C1(i) = i*i;
      RowVector R1 = C1.t();
      Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
      Print(Diff);
   }

   {
      Tracer et1("Stage 4");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(5);
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A1 = X*X;
      ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }


   {
      Tracer et1("Stage 5");
      int n = 40;
      SymmetricBandMatrix B(n,2); B = 0.0;
      for (i=1; i<=n; i++)
      {
         B(i,i) = 6;
         if (i<=n-1) B(i,i+1) = -4;
         if (i<=n-2) B(i,i+2) = 1;
      }
      B(1,1) = 5; B(n,n) = 5;
      SymmetricMatrix A = B;
      ColumnVector X(n);
      X(1) = 429;
      for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
      X = X / 100000L;
      // the matrix B is rather ill-conditioned so the difficulty is getting
      // good agreement (we have chosen X very small) may not be surprising;
      // maximum element size in B.i() is around 1400
      ColumnVector Y1 = A.i() * X;
      LowerTriangularMatrix C1 = Cholesky(A);
      ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
      Clean(Y2, 0.000000001); Print(Y2);
      UpperTriangularMatrix CU = C1.t().i();
      LowerTriangularMatrix CL = C1.i();
      Y2 = CU * (CL * X) - Y1;
      Clean(Y2, 0.000000001); Print(Y2);
      Y2 = B.i() * X - Y1; Clean(Y2, 0.000000001); Print(Y2);

      LowerBandMatrix C2 = Cholesky(B);
      Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
      ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
      Clean(Y3, 0.000000001); Print(Y3);
      CU = C1.t().i();
      CL = C1.i();
      Y3 = CU * (CL * X) - Y1;
      Clean(Y3, 0.000000001); Print(Y3);

      Y3 = B.i() * X - Y1; Clean(Y3, 0.000000001); Print(Y3);

      SymmetricMatrix AI = A.i();
      Y2 = AI*X - Y1; Clean(Y2, 0.000000001); Print(Y2);
      SymmetricMatrix BI = B.i();
      BandMatrix C = B; Matrix CI = C.i();
      M = A.i() - CI; Clean(M, 0.000000001); Print(M);
      M = B.i() - CI; Clean(M, 0.000000001); Print(M);
      M = AI-BI; Clean(M, 0.000000001); Print(M);
      M = AI-CI; Clean(M, 0.000000001); Print(M);

      M = A; AI << M; M = AI-A; Clean(M, 0.000000001); Print(M);
      C = B; BI << C; M = BI-B; Clean(M, 0.000000001); Print(M);
   }

   {
      Tracer et1("Stage 5");
      SymmetricMatrix A(4), B(4);
      A << 5
        << 1 << 4
        << 2 << 1 << 6
        << 1 << 0 << 1 << 7;
      B << 8
        << 1 << 5
        << 1 << 0 << 9
        << 2 << 1 << 0 << 6;