tmta.cc
Go to the documentation of this file.
00001 
00002 //#define WANT_STREAM
00003 
00004 
00005 #include "include.h"
00006 
00007 #include "newmatap.h"
00008 
00009 #include "tmt.h"
00010 
00011 #ifdef use_namespace
00012 using namespace NEWMAT;
00013 #endif
00014 
00015 
00016 /**************************** test program ******************************/
00017 
00018 
00019 static void process(const GeneralMatrix& A,
00020    const ColumnVector& X1, const ColumnVector& X2)
00021 {
00022       Matrix B = A;
00023       LinearEquationSolver L(A);
00024       Matrix Y(4,2);
00025       Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
00026       Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
00027       Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
00028 }
00029 
00030 
00031 
00032 void trymata()
00033 {
00034 //   cout << "\nTenth test of Matrix package\n";
00035    Tracer et("Tenth test of Matrix package");
00036    Tracer::PrintTrace();
00037    int i; int j;
00038    UpperTriangularMatrix U(8);
00039    for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
00040    Matrix X(8,6);
00041    for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
00042    Matrix Y = U.i()*X; Matrix MU=U;
00043    Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
00044    Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
00045    UpperTriangularMatrix UX(8);
00046    for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
00047    UX(4,4)=0; UX(4,5)=0;
00048    UpperTriangularMatrix UY = U.i() * UX;
00049    { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
00050    LowerTriangularMatrix LY = U.t().i() * UX.t();
00051    { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
00052    DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
00053    { X=D.i()*MU; }
00054    { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
00055    { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
00056 //   X=MU.t();
00057 //   LY=D.i()*U.t(); Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
00058 //   LowerTriangularMatrix L=U.t();
00059 //   LY=D.i()*L; Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
00060    U.ReSize(8);
00061    for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
00062    MU = U;
00063    MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
00064    MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
00065 
00066    // test LINEQ
00067    {
00068       ColumnVector X1(4), X2(4);
00069       X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
00070       X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
00071 
00072 
00073       Matrix A(4,4);
00074       A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
00075       A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
00076       A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
00077       A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
00078       process(A,X1,X2);
00079 
00080       BandMatrix B(4,1,1);  B.Inject(A);
00081       process(B,X1,X2);
00082 
00083       UpperTriangularMatrix U(4);
00084       U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
00085                 U(2,2)=8; U(2,3)=7; U(2,4)=6;
00086                           U(3,3)=2; U(3,4)=7;
00087                                     U(4,4)=1;
00088       process(U,X1,X2);
00089 
00090       // check rowwise load
00091       UpperTriangularMatrix U1(4);
00092       U1.Row(1) << 1 << 2 << 3 << 4;
00093       U1.Row(2)      << 8 << 7 << 6;
00094       U1.Row(3)           << 2 << 7;
00095       U1.Row(4)                << 1;
00096 
00097       U1 -= U;
00098 
00099       Print(U1);
00100 
00101       LowerTriangularMatrix L = U.t();
00102       process(L,X1,X2);
00103    }
00104 
00105    // test inversion of poorly conditioned matrix
00106    // a user complained this didn't work under OS9
00107    {
00108       Matrix M(4,4);
00109 
00110       M <<  8.613057e+00 <<  8.693985e+00 << -2.322050e-01  << 0.000000e+00
00111         <<  8.693985e+00 <<  8.793868e+00 << -2.346310e-01  << 0.000000e+00
00112         << -2.322050e-01 << -2.346310e-01 <<  6.264000e-03  << 0.000000e+00
00113         <<  0.000000e+00 <<  0.000000e+00 <<  0.000000e+00  << 3.282806e+03 ;
00114       Matrix MI = M.i();
00115       DiagonalMatrix I(4); I = 1;
00116       Matrix Diff = MI *  M - I;
00117       Clean(Diff,0.00000001); Print(Diff);
00118       // Alternatively do Cholesky
00119       SymmetricMatrix SM; SM << M;
00120       LowerTriangularMatrix LT = Cholesky(SM).i();
00121       MI = LT.t() * LT; Diff = MI *  M - I;
00122       Clean(Diff,0.00000001); Print(Diff);
00123    }
00124 
00125 //   cout << "\nEnd of tenth test\n";
00126 }


rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:13