tmtc.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 
00010 
00011 #include "include.h"
00012 #include "newmatap.h"
00013 
00014 #include "tmt.h"
00015 
00016 #ifdef use_namespace
00017 using namespace NEWMAT;
00018 #endif
00019 
00020 
00021 
00022 
00023 void trymatc()
00024 {
00025 //   cout << "\nTwelfth test of Matrix package\n";
00026    Tracer et("Twelfth test of Matrix package");
00027    Tracer::PrintTrace();
00028    DiagonalMatrix D(15); D=1.5;
00029    Matrix A(15,15);
00030    int i,j;
00031    for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
00032    { A = A + D; }
00033    ColumnVector B(15);
00034    for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
00035    {
00036       Tracer et1("Stage 1");
00037       ColumnVector B1=B;
00038       B=(A*2.0).i() * B1;
00039       Matrix X = A*B-B1/2.0;
00040       Clean(X, 0.000000001); Print(X);
00041       A.ReSize(3,5);
00042       for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
00043 
00044       B = A.AsColumn()+10000;
00045       RowVector R = (A+10000).AsColumn().t();
00046       Print( RowVector(R-B.t()) );
00047    }
00048 
00049    {
00050       Tracer et1("Stage 2");
00051       B = A.AsColumn()+10000;
00052       Matrix XR = (A+10000).AsMatrix(15,1).t();
00053       Print( RowVector(XR-B.t()) );
00054    }
00055 
00056    {
00057       Tracer et1("Stage 3");
00058       B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
00059       Matrix MR = (A+10000).AsColumn().t();
00060       Print( RowVector(MR-B.t()) );
00061 
00062       B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
00063       MR = A.AsColumn().t();
00064       Print( RowVector(MR-B.t()) );
00065    }
00066 
00067    {
00068       Tracer et1("Stage 4");
00069       B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
00070       RowVector R = A.AsColumn().t();
00071       Print( RowVector(R-B.t()) );
00072    }
00073 
00074    {
00075       Tracer et1("Stage 5");
00076       RowVector R = (A.AsColumn()-5000).t();
00077       B = ((R.t()+10000) - A.AsColumn())-5000;
00078       Print( RowVector(B.t()) );
00079    }
00080 
00081    {
00082       Tracer et1("Stage 6");
00083       B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
00084       Print(ColumnVector(B1-B));
00085    }
00086 
00087    {
00088       Tracer et1("Stage 7");
00089       Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
00090       for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
00091       Print(B);
00092    }
00093 
00094    {
00095       Tracer et1("Stage 8");
00096       A.ReSize(7,7); D.ReSize(7);
00097       for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
00098       for (i=1; i<=7; i++) D(i,i) = i;
00099       UpperTriangularMatrix U; U << A;
00100       Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
00101       A.Inject(D); Print(Matrix(X-A));
00102       X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
00103       Print(Matrix(X-A));
00104    }
00105 
00106    {
00107       Tracer et1("Stage 9");
00108       A.ReSize(7,5);
00109       for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
00110       Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
00111       Matrix X = A;
00112       Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
00113       Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
00114    }
00115 
00116    {
00117       Tracer et1("Stage 10");
00118       // some tests on submatrices
00119       UpperTriangularMatrix U(20);
00120       for (i=1; i<=20; i++) for (j=i; j<=20; j++) U(i,j)=100 * i + j;
00121       UpperTriangularMatrix V = U.SymSubMatrix(1,5);
00122       UpperTriangularMatrix U1 = U;
00123       U1.SubMatrix(4,8,5,9) /= 2;
00124       U1.SubMatrix(4,8,5,9) += 388 * V;
00125       U1.SubMatrix(4,8,5,9) *= 2;
00126       U1.SubMatrix(4,8,5,9) += V;
00127       U1 -= U; UpperTriangularMatrix U2 = U1;
00128       U1 << U1.SubMatrix(4,8,5,9);
00129       U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
00130       U1 -= (777*V); Print(U1);
00131 
00132       U1 = U; U1.SubMatrix(4,8,5,9) -= U.SymSubMatrix(1,5);
00133       U1 -= U;  U2 = U1; U1 << U1.SubMatrix(4,8,5,9);
00134       U2.SubMatrix(4,8,5,9) -= U1; Print(U2);
00135       U1 += V; Print(U1);
00136 
00137       U1 = U;
00138       U1.SubMatrix(3,10,15,19) += 29;
00139       U1 -= U;
00140       Matrix X = U1.SubMatrix(3,10,15,19); X -= 29; Print(X);
00141       U1.SubMatrix(3,10,15,19) *= 0; Print(U1);
00142 
00143       LowerTriangularMatrix L = U.t();
00144       LowerTriangularMatrix M = L.SymSubMatrix(1,5);
00145       LowerTriangularMatrix L1 = L;
00146       L1.SubMatrix(5,9,4,8) /= 2;
00147       L1.SubMatrix(5,9,4,8) += 388 * M;
00148       L1.SubMatrix(5,9,4,8) *= 2;
00149       L1.SubMatrix(5,9,4,8) += M;
00150       L1 -= L; LowerTriangularMatrix L2 = L1;
00151       L1 << L1.SubMatrix(5,9,4,8);
00152       L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
00153       L1 -= (777*M); Print(L1);
00154 
00155       L1 = L; L1.SubMatrix(5,9,4,8) -= L.SymSubMatrix(1,5);
00156       L1 -= L; L2 =L1; L1 << L1.SubMatrix(5,9,4,8);
00157       L2.SubMatrix(5,9,4,8) -= L1; Print(L2);
00158       L1 += M; Print(L1);
00159 
00160       L1 = L;
00161       L1.SubMatrix(15,19,3,10) -= 29;
00162       L1 -= L;
00163       X = L1.SubMatrix(15,19,3,10); X += 29; Print(X);
00164       L1.SubMatrix(15,19,3,10) *= 0; Print(L1);
00165    }
00166 
00167    {
00168       Tracer et1("Stage 11");
00169       // more tests on submatrices
00170       Matrix M(20,30);
00171       for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
00172       Matrix M1 = M;
00173 
00174       for (j=1; j<=30; j++)
00175          { ColumnVector CV = 3 * M1.Column(j); M.Column(j) += CV; }
00176       for (i=1; i<=20; i++)
00177          { RowVector RV = 5 * M1.Row(i); M.Row(i) -= RV; }
00178 
00179       M += M1; Print(M);
00180  
00181    }
00182 
00183    {
00184       Tracer et1("Stage 12");
00185       // more tests on Release
00186       Matrix M(20,30);
00187       for (i=1; i<=20; i++) for (j=1; j<=30; j++) M(i,j)=100 * i + j;
00188       Matrix M1 = M;
00189       M.Release();
00190       Matrix M2 = M;
00191       Matrix X = M;   Print(X);
00192       X = M1 - M2;    Print(X);
00193 
00194 #ifndef DONT_DO_NRIC
00195       nricMatrix N = M1;
00196       nricMatrix N1 = N;
00197       N.Release();
00198       nricMatrix N2 = N;
00199       nricMatrix Y = N;   Print(Y);
00200       Y = N1 - N2;        Print(Y);
00201       
00202       N = M1 / 2; N1 = N * 2; N.Release(); N2 = N * 2; Y = N; Print(N);
00203       Y = (N1 - M1) | (N2 - M1); Print(Y);
00204 #endif
00205 
00206    }
00207    
00208    {
00209       Tracer et("Stage 13");
00210       // test sum of squares of rows or columns
00211       MultWithCarry mwc;
00212       DiagonalMatrix DM; Matrix X;
00213       // rectangular matrix
00214       Matrix A(20, 15);
00215       FillWithValues(mwc, A);
00216       // sum of squares of rows
00217       DM << A * A.t();
00218       ColumnVector CV = A.sum_square_rows();
00219       X = CV - DM.AsColumn(); Clean(X, 0.000000001); Print(X);
00220       DM << A.t() * A;
00221       RowVector RV = A.sum_square_columns();
00222       X = RV - DM.AsRow(); Clean(X, 0.000000001); Print(X);
00223       X = RV - A.t().sum_square_rows().t(); Clean(X, 0.000000001); Print(X);
00224       X = CV - A.t().sum_square_columns().t(); Clean(X, 0.000000001); Print(X);
00225       // UpperTriangularMatrix
00226       A.ReSize(17,17); FillWithValues(mwc, A);
00227       UpperTriangularMatrix UT; UT << A;
00228       Matrix A1 = UT;
00229       X = UT.sum_square_rows() - A1.sum_square_rows(); Print(X);
00230       X = UT.sum_square_columns() - A1.sum_square_columns(); Print(X);
00231       // LowerTriangularMatrix
00232       LowerTriangularMatrix LT; LT << A;
00233       A1 = LT;
00234       X = LT.sum_square_rows() - A1.sum_square_rows(); Print(X);
00235       X = LT.sum_square_columns() - A1.sum_square_columns(); Print(X);
00236       // SymmetricMatrix
00237       SymmetricMatrix SM; SM << A;
00238       A1 = SM;
00239       X = SM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00240       X = SM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00241       // DiagonalMatrix
00242       DM << A;
00243       A1 = DM;
00244       X = DM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00245       X = DM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00246       // BandMatrix
00247       BandMatrix BM(17, 3, 5); BM.Inject(A);
00248       A1 = BM;
00249       X = BM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00250       X = BM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00251       // SymmetricBandMatrix
00252       SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
00253       A1 = SBM;
00254       X = SBM.sum_square_rows() - A1.sum_square_rows(); Print(X);
00255       X = SBM.sum_square_columns() - A1.sum_square_columns(); Print(X);
00256       // IdentityMatrix
00257       IdentityMatrix IM(29);
00258       X = IM.sum_square_rows() - 1; Print(X);
00259       X = IM.sum_square_columns() - 1; Print(X);
00260       // Matrix with zero rows
00261       A1.ReSize(0,10);
00262       X.ReSize(1,10); X = 0; X -= A1.sum_square_columns(); Print(X);
00263       X.ReSize(0,1); X -= A1.sum_square_rows(); Print(X);
00264       // Matrix with zero columns
00265       A1.ReSize(10,0);
00266       X.ReSize(10,1); X = 0; X -= A1.sum_square_rows(); Print(X);
00267       X.ReSize(1,0); X -= A1.sum_square_columns(); Print(X);
00268       
00269    }
00270    
00271    {
00272       Tracer et("Stage 14");
00273       // test extend orthonormal
00274       MultWithCarry mwc;
00275       Matrix A(20,5); FillWithValues(mwc, A);
00276       // Orthonormalise
00277       UpperTriangularMatrix R;
00278       Matrix A_old = A;
00279       QRZ(A,R);
00280       // Check decomposition
00281       Matrix X = A * R - A_old; Clean(X, 0.000000001); Print(X);
00282       // Check orthogonality
00283       X = A.t() * A - IdentityMatrix(5);
00284       Clean(X, 0.000000001); Print(X);
00285       // Try orthonality extend 
00286       SquareMatrix A1(20);
00287       A1.Columns(1,5) = A;
00288       extend_orthonormal(A1,5);
00289       // check columns unchanged
00290       X = A - A1.Columns(1,5); Print(X);
00291       // Check orthogonality
00292       X = A1.t() * A1 - IdentityMatrix(20);
00293       Clean(X, 0.000000001); Print(X); 
00294       X = A1 * A1.t() - IdentityMatrix(20);
00295       Clean(X, 0.000000001); Print(X);
00296       // Test with smaller number of columns 
00297       Matrix A2(20,15);
00298       A2.Columns(1,5) = A;
00299       extend_orthonormal(A2,5);
00300       // check columns unchanged
00301       X = A - A2.Columns(1,5); Print(X);
00302       // Check orthogonality
00303       X = A2.t() * A2 - IdentityMatrix(15);
00304       Clean(X, 0.000000001); Print(X);
00305       // check it works with no columns to start with
00306       A2.ReSize(100,100);
00307       extend_orthonormal(A2,0);
00308       // Check orthogonality
00309       X = A2.t() * A2 - IdentityMatrix(100);
00310       Clean(X, 0.000000001); Print(X);
00311       X = A2 * A2.t() - IdentityMatrix(100);
00312       Clean(X, 0.000000001); Print(X);
00313    }
00314 
00315    {
00316       Tracer et("Stage 15");
00317       // test extend orthonormal for SVD
00318       MultWithCarry mwc;
00319       Matrix A(15,40); FillWithValues(mwc, A);
00320       Matrix U, V; DiagonalMatrix D;
00321       SVD(A.t(),D,V,U);
00322       Matrix VE(40,40); VE.columns(1,15) = V;
00323       extend_orthonormal(VE,15);
00324       Matrix DE(15,40); DE = 0; DE.columns(1,15) = D;
00325       Matrix B = U * DE * VE.t();
00326       B -= A;
00327       Clean(B, 0.000000001); Print(B);
00328    }
00329       
00330    {
00331       Tracer et("Stage 16");
00332       // test sum of rows or columns
00333       MultWithCarry mwc;
00334       ColumnVector CVX; Matrix X;
00335       // rectangular matrix
00336       Matrix A(20, 15);
00337       FillWithValues(mwc, A);
00338       // sum of rows
00339       ColumnVector Ones(15); Ones = 1;
00340       CVX = A * Ones;
00341       ColumnVector CV = A.sum_rows();
00342       X = CV - CVX; Clean(X, 0.000000001); Print(X);
00343       Ones.resize(20); Ones = 1;
00344       CVX << A.t() * Ones;
00345       RowVector RV = A.sum_columns();
00346       X = RV - CVX.AsRow(); Clean(X, 0.000000001); Print(X);
00347       X = RV - A.t().sum_rows().t(); Clean(X, 0.000000001); Print(X);
00348       X = CV - A.t().sum_columns().t(); Clean(X, 0.000000001); Print(X);
00349       // UpperTriangularMatrix
00350       A.ReSize(17,17); FillWithValues(mwc, A);
00351       UpperTriangularMatrix UT; UT << A;
00352       Matrix A1 = UT;
00353       X = UT.sum_rows() - A1.sum_rows(); Print(X);
00354       X = UT.sum_columns() - A1.sum_columns(); Print(X);
00355       // LowerTriangularMatrix
00356       LowerTriangularMatrix LT; LT << A;
00357       A1 = LT;
00358       X = LT.sum_rows() - A1.sum_rows(); Print(X);
00359       X = LT.sum_columns() - A1.sum_columns(); Print(X);
00360       // SymmetricMatrix
00361       SymmetricMatrix SM; SM << A;
00362       A1 = SM;
00363       X = SM.sum_rows() - A1.sum_rows(); Print(X);
00364       X = SM.sum_columns() - A1.sum_columns(); Print(X);
00365       // DiagonalMatrix
00366       DiagonalMatrix DM; DM << A;
00367       A1 = DM;
00368       X = DM.sum_rows() - A1.sum_rows(); Print(X);
00369       X = DM.sum_columns() - A1.sum_columns(); Print(X);
00370       // BandMatrix
00371       BandMatrix BM(17, 3, 5); BM.Inject(A);
00372       A1 = BM;
00373       X = BM.sum_rows() - A1.sum_rows(); Print(X);
00374       X = BM.sum_columns() - A1.sum_columns(); Print(X);
00375       // SymmetricBandMatrix
00376       SymmetricBandMatrix SBM(17, 4); SBM.Inject(A);
00377       A1 = SBM;
00378       X = SBM.sum_rows() - A1.sum_rows(); Print(X);
00379       X = SBM.sum_columns() - A1.sum_columns(); Print(X);
00380       // IdentityMatrix
00381       IdentityMatrix IM(29);
00382       X = IM.sum_rows() - 1; Print(X);
00383       X = IM.sum_columns() - 1; Print(X);
00384       // Matrix with zero rows
00385       A1.ReSize(0,10);
00386       X.ReSize(1,10); X = 0; X -= A1.sum_columns(); Print(X);
00387       X.ReSize(0,1); X -= A1.sum_rows(); Print(X);
00388       // Matrix with zero columns
00389       A1.ReSize(10,0);
00390       X.ReSize(10,1); X = 0; X -= A1.sum_rows(); Print(X);
00391       X.ReSize(1,0); X -= A1.sum_columns(); Print(X);
00392       
00393    }
00394    
00395  
00396 
00397 //   cout << "\nEnd of twelfth test\n";
00398 }
00399 
00400 


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