kni: tmt9.cpp Source File

Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 //#define WANT_STREAM
00008 
00009 #include "include.h"
00010 #include "newmatap.h"
00011 
00012 #include "tmt.h"
00013 
00014 #ifdef use_namespace
00015 using namespace NEWMAT;
00016 #endif
00017 
00018 
00019 /**************************** test program ******************************/
00020 
00021 
00022 void trymat9()
00023 {
00024    Tracer et("Ninth test of Matrix package");
00025    Tracer::PrintTrace();
00026 
00027 
00028    int i; int j;
00029    Matrix A(7,7); Matrix X(7,3);
00030    for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
00031    for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
00032    Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
00033    {
00034       Tracer et1("Stage 1");
00035       Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
00036       Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
00037       Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
00038    }
00039    for (i=1;i<=7;i++) D(i,i)=i*i+1;
00040    DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
00041    {
00042       Tracer et1("Stage 2");
00043       Matrix DXE = D.i() * X * E;
00044       DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE); 
00045       E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
00046    }
00047    DiagonalMatrix F=D;
00048    {
00049       Tracer et1("Stage 3");
00050       F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
00051       F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
00052    }
00053    {
00054       Tracer et1("Stage 4");
00055       F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
00056    }
00057    {
00058       Tracer et1("Stage 5");
00059       // testing equal
00060       ColumnVector A(18), B(18);
00061       Matrix X(3,3);
00062       X << 3 << 5  << 7 << 5 << 8 << 2 << 7 << 2 << 9;
00063       SymmetricMatrix S; S << X;
00064       B(1) = S == X;         A(1) = true;
00065       B(2) = S == (X+1);     A(2) = false;
00066       B(3) = (S+2) == (X+2); A(3) = true;
00067       Matrix Y = X;
00068       B(4) = X == Y;         A(4) = true;
00069       B(5) = (X*2) == (Y*2); A(5) = true;
00070       Y(3,3) = 10;
00071       B(6) = X == Y;         A(6) = false;
00072       B(7) = (X*2) == (Y*2); A(7) = false;
00073       B(8) = S == Y;         A(8) = false;
00074       B(9) = S == S;         A(9) = true;
00075       Matrix Z = X.SubMatrix(1,2,2,3);
00076       B(10) = X == Z;        A(10) = false;
00077       GenericMatrix GS = S;
00078       GenericMatrix GX = X;
00079       GenericMatrix GY = Y;
00080       B(11) = GS == GX;      A(11) = true;
00081       B(12) = GS == GY;      A(12) = false;
00082       CroutMatrix CS = S;
00083       CroutMatrix CX = X;
00084       CroutMatrix CY = Y;
00085       B(13) = CS == CX;      A(13) = true;
00086       B(14) = CS == CY;      A(14) = false;
00087       B(15) = X == CX;       A(15) = false;
00088       B(16) = X == A;        A(16) = false;
00089       B(17) = X == (X | X);  A(17) = false;
00090       B(18) = CX == X;       A(18) = false;
00091       A = A - B; Print(A);
00092    }
00093    {
00094       Tracer et1("Stage 6");
00095       // testing equal
00096       ColumnVector A(22), B(22);
00097       BandMatrix X(6,2,1);
00098       X(1,1)=23; X(1,2)=21;
00099       X(2,1)=12; X(2,2)=17; X(2,3)=45;
00100       X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
00101                  X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
00102                             X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
00103                                        X(6,4)=49; X(6,5)=41; X(6,6)=17;
00104       SymmetricBandMatrix S1(6,2); S1.Inject(X);
00105       BandMatrix U(6,2,3); U = 0.0; U.Inject(X);
00106       B(1) = U == X;         A(1) = true;
00107       B(2) = U == (X*3);     A(2) = false;
00108       B(3) = (U*5) == (X*5); A(3) = true;
00109       Matrix Y = X;
00110       B(4) = X == Y;         A(4) = true;
00111       B(5) = (X*2) == (Y*2); A(5) = true;
00112       Y(6,6) = 10;
00113       B(6) = X == Y;         A(6) = false;
00114       B(7) = (X*2) == (Y*2); A(7) = false;
00115       B(8) = U == Y;         A(8) = false;
00116       B(9) = U == U;         A(9) = true;
00117       Matrix Z = X.SubMatrix(1,2,2,3);
00118       B(10) = X == Z;        A(10) = false;
00119       GenericMatrix GU = U;
00120       GenericMatrix GX = X;
00121       GenericMatrix GY = Y;
00122       B(11) = GU == GX;      A(11) = true;
00123       B(12) = GU == GY;      A(12) = false;
00124       X = X + X.t(); U = U + U.t();
00125       SymmetricBandMatrix S(6,2); S.Inject(X);
00126       Matrix D = S-X; Print(D);
00127       BandLUMatrix BS = S;
00128       BandLUMatrix BX = X;
00129       BandLUMatrix BU = U;
00130       CroutMatrix CX = X;
00131       B(13) = BS == BX;      A(13) = true;
00132       B(14) = BX == BU;      A(14) = false;
00133       B(15) = X == BX;       A(15) = false;
00134       B(16) = X != BX;       A(16) = true;
00135       B(17) = BX != BS;      A(17) = false;
00136       B(18) = (2*X) != (X*2);A(18) = false;
00137       B(19) = (X*2) != (X+2);A(19) = true;
00138       B(20) = BX == CX;      A(20) = false;
00139       B(21) = CX == BX;      A(21) = false;
00140       B(22) = BX == X;       A(22) = false;
00141       A = A - B; Print(A);
00142       DiagonalMatrix I(6); I=1.0;
00143       D = BS.i() * X - I;  Clean(D,0.00000001); Print(D);
00144       D = BX.i() * X - I;  Clean(D,0.00000001); Print(D);
00145       D = BU.i() * X - I;  Clean(D,0.00000001); Print(D);
00146 
00147       // test row wise load
00148       SymmetricBandMatrix X1(6,2);
00149       X1.Row(1) << 23;
00150       X1.Row(2) << 12 << 17;
00151       X1.Row(3) << 35 << 19 << 24;
00152       X1.Row(4)       << 17 << 11 << 19;
00153       X1.Row(5)             << 10 << 44 << 23;
00154       X1.Row(6)                   << 49 << 41 << 17;
00155       Matrix M = X1 - S1; Print(M);
00156 
00157       // check out submatrix
00158       SymmetricBandMatrix X2(20,3); X2 = 0.0;
00159       X2.SubMatrix(2,7,2,7) = X1; X2.SymSubMatrix(11,16) = 2 * X1;
00160       Matrix MX1 = X1;
00161       Matrix MX2(20,20); MX2 = 0;
00162       MX2.SymSubMatrix(2,7) = MX1; MX2.SubMatrix(11,16,11,16) = MX1 * 2;
00163       MX2 -= X2; Print(MX2);
00164 
00165       BandMatrix X4(20,3,3); X4 = 0.0;
00166       X4.SubMatrix(2,7,3,8) = X1; X4.SubMatrix(11,16,10,15) = 2 * X1;
00167       MX1 = X1;
00168       Matrix MX4(20,20); MX4 = 0;
00169       MX4.SubMatrix(2,7,3,8) = MX1; MX4.SubMatrix(11,16,10,15) = MX1 * 2;
00170       MX4 -= X4; Print(MX4);
00171 
00172       MX1 = X1.i() * X1 - IdentityMatrix(6);
00173       Clean(MX1,0.00000001); Print(MX1);
00174 
00175    }
00176 
00177    {
00178       Tracer et1("Stage 7");
00179       // testing equal
00180       ColumnVector A(12), B(12);
00181       BandMatrix X(6,2,1);
00182       X(1,1)=23; X(1,2)=21;
00183       X(2,1)=12; X(2,2)=17; X(2,3)=45;
00184       X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
00185                  X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
00186                             X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
00187                                        X(6,4)=49; X(6,5)=41; X(6,6)=17;
00188       Matrix Y = X;
00189       LinearEquationSolver LX = X;
00190       LinearEquationSolver LY = Y;
00191       CroutMatrix CX = X;
00192       CroutMatrix CY = Y;
00193       BandLUMatrix BX = X;
00194       B(1) = LX == CX;       A(1) = false;
00195       B(2) = LY == CY;       A(2) = true;
00196       B(3) = X == Y;         A(3) = true;
00197       B(4) = BX == LX;       A(4) = true;
00198       B(5) = CX == CY;       A(5) = true;
00199       B(6) = LX == LY;       A(6) = false;
00200       B(7) = BX == BX;       A(7) = true;
00201       B(8) = CX == CX;       A(8) = true;
00202       B(9) = LX == LX;       A(9) = true;
00203       B(10) = LY == LY;      A(10) = true;
00204       CroutMatrix CX1 = X.SubMatrix(1,4,1,4);
00205       B(11) = CX == CX1;     A(11) = false;
00206       BandLUMatrix BX1 = X.SymSubMatrix(1,4);    // error with SubMatrix
00207       B(12) = BX == BX1;     A(12) = false;
00208       A = A - B; Print(A);
00209       DiagonalMatrix I(6); I=1.0; Matrix D;
00210       D = LX.i() * X - I;  Clean(D,0.00000001); Print(D);
00211       D = LY.i() * X - I;  Clean(D,0.00000001); Print(D);
00212       I.ReSize(4); I = 1;
00213       D = CX1.i() * X.SymSubMatrix(1,4) - I;  Clean(D,0.00000001); Print(D);
00214       D = BX1.i() * X.SubMatrix(1,4,1,4) - I;  Clean(D,0.00000001); Print(D);
00215    }
00216 
00217    {
00218       Tracer et1("Stage 8");
00219       // test copying CroutMatrix and BandLUMatrix - see also tmtd.cpp 
00220       MultWithCarry MWC;
00221       SymmetricBandMatrix SBM(50, 10);
00222       for (int i = 1; i <= 50; ++i) for (int j = 1; j <= i; ++j)
00223          if (i - j <= 10) SBM(i, j) = MWC.Next();
00224       CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
00225       CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
00226       CM1.release(); BM1.release();
00227       CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
00228       Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
00229       Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
00230       X *= SBM; X -= IdentityMatrix(50); Clean(X,0.00000001); Print(X);
00231       LogAndSign x = log_determinant(SBM);
00232       LogAndSign y = log_determinant(CM2);
00233       LogAndSign z = log_determinant(BM2);
00234       RowVector D(4);
00235       D(1) = y.value() - x.value();
00236       D(2) = z.value() - x.value();
00237       D(3) = y.sign() - x.sign();
00238       D(4) = z.sign() - x.sign();
00239       Clean(D,0.00000001); Print(D);
00240    }
00241 
00242    {
00243       Tracer et1("Stage 9");
00244       // do it again odd matrix size 
00245       MultWithCarry MWC;
00246       SymmetricBandMatrix SBM(51, 10);
00247       for (int i = 1; i <= 51; ++i) for (int j = 1; j <= i; ++j)
00248          if (i - j <= 10) SBM(i, j) = MWC.Next();
00249       CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
00250       CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
00251       CM1.release(); BM1.release();
00252       CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
00253       Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
00254       Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
00255       X *= SBM; X -= IdentityMatrix(51); Clean(X,0.00000001); Print(X);
00256       LogAndSign x = log_determinant(SBM);
00257       LogAndSign y = log_determinant(CM2);
00258       LogAndSign z = log_determinant(BM2);
00259       RowVector D(4);
00260       D(1) = y.value() - x.value();
00261       D(2) = z.value() - x.value();
00262       D(3) = y.sign() - x.sign();
00263       D(4) = z.sign() - x.sign();
00264       Clean(D,0.00000001); Print(D);
00265    }
00266    
00267    
00268 //   cout << "\nEnd of ninth test\n";
00269 }
00270 
00271