tmt4.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 
00010 
00011 #include "include.h"
00012 
00013 #include "newmat.h"
00014 
00015 #include "tmt.h"
00016 
00017 #ifdef use_namespace
00018 using namespace NEWMAT;
00019 #endif
00020 
00021 
00022 /**************************** test program ******************************/
00023 
00024 
00025 void trymat4()
00026 {
00027 //   cout << "\nFourth test of Matrix package\n";
00028    Tracer et("Fourth test of Matrix package");
00029    Tracer::PrintTrace();
00030 
00031 
00032    {
00033       Tracer et1("Stage 1");
00034       int i, j;
00035       Matrix M(10,10);
00036       UpperTriangularMatrix U(10);
00037       for (i=1;i<=10;i++) for (j=1;j<=10;j++) M(i,j) = 100*i+j;
00038       U << -M;
00039       Matrix X1 = M.Rows(2,4);
00040       Matrix Y1 = U.t().Rows(2,4);
00041       Matrix X = U; { Print(Matrix(X.Columns(2,4).t()-Y1)); }
00042       RowVector RV = M.Row(5);
00043       {
00044          X.ReSize(3,10);
00045          X.Row(1) << M.Row(2); X.Row(2) << M.Row(3); X.Row(3) << M.Row(4);
00046          Print(Matrix(X-X1));
00047       }
00048       {
00049          UpperTriangularMatrix V = U.SymSubMatrix(3,5);
00050          Matrix MV = U.SubMatrix(3,5,3,5); { Print(Matrix(MV-V)); }
00051          Matrix X2 = M.t().Columns(2,4); { Print(Matrix(X2-X1.t())); }
00052          Matrix Y2 = U.Columns(2,4); { Print(Matrix(Y2-Y1.t())); }
00053          ColumnVector CV = M.t().Column(5); { Print(ColumnVector(CV-RV.t())); }
00054          X.ReSize(10,3); M = M.t();
00055          X.Column(1) << M.Column(2); X.Column(2) << M.Column(3);
00056          X.Column(3) << M.Column(4);
00057          Print(Matrix(X-X2));
00058       }
00059    }
00060 
00061    {
00062       Tracer et1("Stage 2");
00063       int i, j;
00064       Matrix M; Matrix X; M.ReSize(5,8);
00065       for (i=1;i<=5;i++) for (j=1;j<=8;j++) M(i,j) = 100*i+j;
00066       {
00067          X = M.Columns(5,8); M.Columns(5,8) << M.Columns(1,4);
00068              M.Columns(1,4) << X;
00069          X = M.Columns(3,4); M.Columns(3,4) << M.Columns(1,2);
00070              M.Columns(1,2) << X;
00071          X = M.Columns(7,8); M.Columns(7,8) << M.Columns(5,6);
00072              M.Columns(5,6) << X;
00073       }
00074       {
00075          X = M.Column(2); M.Column(2) = M.Column(1); M.Column(1) = X;
00076          X = M.Column(4); M.Column(4) = M.Column(3); M.Column(3) = X;
00077          X = M.Column(6); M.Column(6) = M.Column(5); M.Column(5) = X;
00078          X = M.Column(8); M.Column(8) = M.Column(7); M.Column(7) = X;
00079          X.ReSize(5,8);
00080       }
00081       for (i=1;i<=5;i++) for (j=1;j<=8;j++) X(i,9-j) = 100*i+j;
00082       Print(Matrix(X-M));
00083    }
00084    {
00085       Tracer et1("Stage 3");
00086       // try submatrices of zero dimension
00087       int i, j;
00088       Matrix A(4,5); Matrix B, C;
00089       for (i=1; i<=4; i++) for (j=1; j<=5; j++)
00090          A(i,j) = 100+i*10+j;
00091       B = A + 100;
00092       C = A | B.Columns(4,3); Print(Matrix(A - C));
00093       C = A | B.Columns(1,0); Print(Matrix(A - C));
00094       C = A | B.Columns(6,5); Print(Matrix(A - C));
00095       C = A & B.Rows(2,1); Print(Matrix(A - C));
00096    }
00097    {
00098       Tracer et1("Stage 4");
00099       BandMatrix BM(5,3,2);
00100       BM(1,1) = 1; BM(1,2) = 2; BM(1,3) = 3;
00101       BM(2,1) = 4; BM(2,2) = 5; BM(2,3) = 6; BM(2,4) = 7;
00102       BM(3,1) = 8; BM(3,2) = 9; BM(3,3) =10; BM(3,4) =11; BM(3,5) =12;
00103       BM(4,1) =13; BM(4,2) =14; BM(4,3) =15; BM(4,4) =16; BM(4,5) =17;
00104                    BM(5,2) =18; BM(5,3) =19; BM(5,4) =20; BM(5,5) =21;
00105       SymmetricBandMatrix SM(5,3);
00106       SM.Inject(BandMatrix(BM + BM.t()));
00107       Matrix A = BM + 1;
00108       Matrix M = A + A.t() - 2;
00109       Matrix C = A.i() * BM;
00110       C = A * C - BM; Clean(C, 0.000000001); Print(C);
00111       C = A.i() * SM;
00112       C = A * C - M; Clean(C, 0.000000001); Print(C);
00113 
00114       // check row-wise load
00115       BandMatrix BM1(5,3,2);
00116       BM1.Row(1) <<  1 <<  2 <<  3;
00117       BM1.Row(2) <<  4 <<  5 <<  6 <<  7;
00118       BM1.Row(3) <<  8 <<  9 << 10 << 11 << 12;
00119       BM1.Row(4) << 13 << 14 << 15 << 16 << 17;
00120       BM1.Row(5)       << 18 << 19 << 20 << 21;
00121       Matrix M1 = BM1 - BM; Print(M1);
00122    }
00123    {
00124       Tracer et1("Stage 5");
00125       Matrix X(4,4);
00126       X << 1 << 2 << 3 << 4
00127         << 5 << 6 << 7 << 8
00128         << 9 <<10 <<11 <<12
00129         <<13 <<14 <<15 <<16;
00130       Matrix Y(4,0);
00131       Y = X | Y;
00132       X -= Y; Print(X);
00133 
00134       DiagonalMatrix D(1);
00135       D << 23;                       // matrix input with just one value
00136       D(1) -= 23; Print(D);
00137 
00138    }
00139    {
00140       Tracer et1("Stage 6");
00141       Matrix h (2,2);
00142       h << 1.0 << 2.0 << 0.0 << 1.0 ;
00143       RowVector c(2);
00144       c << 0.0 << 1.0;
00145       h -= c & c;
00146       h -= c.t().Reverse() | c.Reverse().t();
00147       Print(h);
00148    }
00149    {
00150       Tracer et1("Stage 7");
00151       // Check row-wise input for diagonal matrix
00152       DiagonalMatrix D(4);
00153       D << 18 << 23 << 31 << 17;
00154       DiagonalMatrix D1(4);
00155       D1.Row(1) << 18; D1.Row(2) << 23; D1.Row(3) << 31; D1.Row(4) << 17;
00156       D1 -= D; Print(D1);
00157       D1(1) = 18; D1(2) = 23; D1(3) = 31; D1(4) = 17;
00158       D1 -= D; Print(D1);
00159    }
00160    
00161    {
00162       Tracer et1("Stage 8");
00163       // test swap functions
00164       MultWithCarry MWC;
00165       Matrix A(3,4); Matrix B(5,6);
00166       FillWithValues(MWC, A); FillWithValues(MWC, B);
00167       Matrix A1 = A; Matrix B1 = B; A.Release(); B.Release(2);
00168       swap(A, B);
00169       int a = A.size() - B1.size(), b = B.size() - A1.size();
00170       Matrix D = A - B1; Print(D);
00171       D = B - A1; Print(D);
00172       Print(B);   // should be zero because of release
00173       D = A - B1; Print(D);
00174       Print(A);   // now should be zero
00175       D.ReSize(1,2); D(1,1) = a; D(1,2) = b; Print(D);
00176       
00177       A.ReSize(20,20); FillWithValues(MWC, A);
00178      
00179       UpperTriangularMatrix UA; UA << A; UpperTriangularMatrix UA1 = UA;
00180       UpperTriangularMatrix UB;
00181       swap(UA, UB); Print(UA); UB -= UA1; Print(UB);
00182       
00183       LowerTriangularMatrix LA; LA << A; LowerTriangularMatrix LA1 = LA;
00184       LowerTriangularMatrix LB;
00185       swap(LB, LA); Print(LA); LB -= LA1; Print(LB);
00186 
00187       SymmetricMatrix SA; SA << A; SymmetricMatrix SA1 = SA;
00188       SymmetricMatrix SB;
00189       swap(SA, SB); Print(SA); SB -= SA1; Print(SB);
00190       
00191       DiagonalMatrix DA; DA << A; DiagonalMatrix DA1 = DA;
00192       DiagonalMatrix DB;
00193       swap(DB, DA); Print(DA); DB -= DA1; Print(DB);
00194       
00195       RowVector RVA = A.Row(1); RowVector RVA1 = RVA;
00196       RowVector RVB;
00197       swap(RVB, RVA); Print(RVA); RVB -= RVA1; Print(RVB);
00198       
00199       ColumnVector CVA = A.Column(1); ColumnVector CVA1 = CVA;
00200       ColumnVector CVB;
00201       swap(CVA, CVB); Print(CVA); CVB -= CVA1; Print(CVB);
00202       
00203       BandMatrix BA(20, 7, 4); BA.Inject(A); BandMatrix BA1 = BA;
00204       BandMatrix BB;
00205       swap(BA, BB); D = BA; Print(D); BB -= BA1; D = BB; Print(D);
00206       
00207       LowerBandMatrix LBA(20, 6); LBA.Inject(A); LowerBandMatrix LBA1 = LBA;
00208       LowerBandMatrix LBB;
00209       swap(LBB, LBA); D = LBA; Print(D); LBB -= LBA1; D = LBB; Print(D);
00210       
00211       UpperBandMatrix UBA(20, 9); UBA.Inject(A); UpperBandMatrix UBA1 = UBA;
00212       UpperBandMatrix UBB;
00213       swap(UBA, UBB); D = UBA; Print(D); UBB -= UBA1; D = UBB; Print(D);
00214       
00215       SymmetricBandMatrix SBA(20, 4); SBA.Inject(A);
00216       SymmetricBandMatrix SBA1 = SBA;
00217       SymmetricBandMatrix SBB;
00218       
00219       swap(SBB, SBA); D = SBA; Print(D);
00220       SBB -= SBA1; D = SBB; Print(D);
00221       
00222       B.ReSize(10,10); FillWithValues(MWC, B);
00223       
00224       CroutMatrix CA = A; IdentityMatrix IA(20);
00225       CroutMatrix CB = B; IdentityMatrix IB(10);
00226       swap(CA, CB); swap(IA, IB);
00227       D = CA.i() * B - IA; Clean(D,0.00000001); Print(D);
00228       D = CB.i() * A - IB; Clean(D,0.00000001); Print(D);
00229       
00230       BA.ReSize(20, 5, 7); BA.Inject(A); BandLUMatrix BLUA = BA;
00231       BB.ReSize(10, 3, 4); BB.Inject(B); BandLUMatrix BLUB = BB;
00232       swap(BLUA, BLUB);
00233       D = BLUA.i() * BB - IA; Clean(D,0.00000001); Print(D);
00234       D = BLUB.i() * BA - IB; Clean(D,0.00000001); Print(D);
00235 
00236       
00237       SBA.ReSize(20, 5); SBA.Inject(A); BandLUMatrix SBLUA = SBA;
00238       SBB.ReSize(10, 3); SBB.Inject(B); BandLUMatrix SBLUB = SBB;
00239       swap(SBLUA, SBLUB);
00240       D = SBLUA.i() * SBB - IA; Clean(D,0.00000001); Print(D);
00241       D = SBLUB.i() * SBA - IB; Clean(D,0.00000001); Print(D);
00242       
00243       UA << A;
00244       GenericMatrix GUA = UA; GenericMatrix GB = B; swap(GUA, GB);
00245       D = GB - UA; Print(D); D = B - GUA; Print(D);
00246       
00247    }
00248 
00249 //   cout << "\nEnd of fourth test\n";
00250 }
00251 
00252 


kni
Author(s): Martin Günther
autogenerated on Thu Aug 27 2015 13:40:07