tmt8.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 
00010 #include "include.h"
00011 
00012 #include "newmatap.h"
00013 
00014 #include "tmt.h"
00015 
00016 #ifdef use_namespace
00017 using namespace NEWMAT;
00018 #endif
00019 
00020 
00021 
00022 // **************************** test program ******************************
00023 
00024 
00025 void Transposer(const GenericMatrix& GM1, GenericMatrix&GM2)
00026    { GM2 = GM1.t(); }
00027 
00028 // this is a routine in "Numerical Recipes in C" format
00029 // if R is a row vector, C a column vector and D diagonal
00030 // make matrix DCR
00031 
00032 static void DCR(Real d[], Real c[], int m, Real r[], int n, Real **dcr)
00033 {
00034    int i, j;
00035    for (i = 1; i <= m; i++) for (j = 1; j <= n; j++)
00036    dcr[i][j] = d[i] * c[i] * r[j];
00037 }
00038 
00039 ReturnMatrix TestReturn(const GeneralMatrix& gm) { return gm; }
00040 
00041 void trymat8()
00042 {
00043 //   cout << "\nEighth test of Matrix package\n";
00044    Tracer et("Eighth test of Matrix package");
00045    Tracer::PrintTrace();
00046 
00047    int i;
00048 
00049 
00050    DiagonalMatrix D(6);
00051    for (i=1;i<=6;i++)  D(i,i)=i*i+i-10;
00052    DiagonalMatrix D2=D;
00053    Matrix MD=D;
00054 
00055    DiagonalMatrix D1(6); for (i=1;i<=6;i++) D1(i,i)=-100+i*i*i;
00056    Matrix MD1=D1;
00057    Print(Matrix(D*D1-MD*MD1));
00058    Print(Matrix((-D)*D1+MD*MD1));
00059    Print(Matrix(D*(-D1)+MD*MD1));
00060    DiagonalMatrix DX=D;
00061    {
00062       Tracer et1("Stage 1");
00063       DX=(DX+D1)*DX; Print(Matrix(DX-(MD+MD1)*MD));
00064       DX=D;
00065       DX=-DX*DX+(DX-(-D1))*((-D1)+DX);
00066       // Matrix MX = Matrix(MD1);
00067       // MD1=DX+(MX.t())*(MX.t()); Print(MD1);
00068       MD1=DX+(Matrix(MD1).t())*(Matrix(MD1).t()); Print(MD1);
00069       DX=D; DX=DX; DX=D2-DX; Print(DiagonalMatrix(DX));
00070       DX=D;
00071    }
00072    {
00073       Tracer et1("Stage 2");
00074       D.Release(2);
00075       D1=D; D2=D;
00076       Print(DiagonalMatrix(D1-DX));
00077       Print(DiagonalMatrix(D2-DX));
00078       MD1=1.0;
00079       Print(Matrix(MD1-1.0));
00080    }
00081    {
00082       Tracer et1("Stage 3");
00083       //GenericMatrix
00084       LowerTriangularMatrix LT(4);
00085       LT << 1 << 2 << 3 << 4 << 5 << 6  << 7 << 8 << 9 << 10;
00086       UpperTriangularMatrix UT = LT.t() * 2.0;
00087       GenericMatrix GM1 = LT;
00088       LowerTriangularMatrix LT1 = GM1-LT; Print(LT1);
00089       GenericMatrix GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
00090       GM2 = GM1; LT1 = GM2; LT1 = LT1-LT; Print(LT1);
00091       GM2 = GM1*2; LT1 = GM2; LT1 = LT1-LT*2; Print(LT1);
00092       GM1.Release();
00093       GM1=GM1; LT1=GM1-LT; Print(LT1); LT1=GM1-LT; Print(LT1);
00094       GM1.Release();
00095       GM1=GM1*4; LT1=GM1-LT*4; Print(LT1);
00096       LT1=GM1-LT*4; Print(LT1); GM1.CleanUp();
00097       GM1=LT; GM2=UT; GM1=GM1*GM2; Matrix M=GM1; M=M-LT*UT; Print(M);
00098       Transposer(LT,GM2); LT1 = LT - GM2.t(); Print(LT1);
00099       GM1=LT; Transposer(GM1,GM2); LT1 = LT - GM2.t(); Print(LT1);
00100       GM1 = LT; GM1 = GM1 + GM1; LT1 = LT*2-GM1; Print(LT1);
00101       DiagonalMatrix D; D << LT; GM1 = D; LT1 = GM1; LT1 -= D; Print(LT1);
00102       UpperTriangularMatrix UT1 = GM1; UT1 -= D; Print(UT1);
00103    }
00104    {
00105       Tracer et1("Stage 4");
00106       // Another test of SVD
00107       Matrix M(12,12); M = 0;
00108       M(1,1) = M(2,2) = M(4,4) = M(6,6) =
00109          M(7,7) = M(8,8) = M(10,10) = M(12,12) = -1;
00110       M(1,6) = M(1,12) = -5.601594;
00111       M(3,6) = M(3,12) = -0.000165;
00112       M(7,6) = M(7,12) = -0.008294;
00113       DiagonalMatrix D;
00114       SVD(M,D);
00115       SortDescending(D);
00116       // answer given by matlab
00117       DiagonalMatrix DX(12);
00118       DX(1) = 8.0461;
00119       DX(2) = DX(3) = DX(4) = DX(5) = DX(6) = DX(7) = 1;
00120       DX(8) = 0.1243;
00121       DX(9) = DX(10) = DX(11) = DX(12) = 0;
00122       D -= DX; Clean(D,0.0001); Print(D);
00123    }
00124 #ifndef DONT_DO_NRIC
00125    {
00126       Tracer et1("Stage 5");
00127       // test numerical recipes in C interface
00128       DiagonalMatrix D(10);
00129       D << 1 << 4 << 6 << 2 << 1 << 6 << 4 << 7 << 3 << 1;
00130       ColumnVector C(10);
00131       C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
00132       RowVector R(6);
00133       R << 2 << 3 << 5 << 7 << 11 << 13;
00134       nricMatrix M(10, 6);
00135       DCR( D.nric(), C.nric(), 10, R.nric(), 6, M.nric() );
00136       M -= D * C * R;  Print(M);
00137 
00138       D.ReSize(5);
00139       D << 1.25 << 4.75 << 9.5 << 1.25 << 3.75;
00140       C.ReSize(5);
00141       C << 1.5 << 7.5 << 4.25 << 0.0 << 7.25;
00142       R.ReSize(9);
00143       R << 2.5 << 3.25 << 5.5 << 7 << 11.25 << 13.5 << 0.0 << 1.5 << 3.5;
00144       Matrix MX = D * C * R;
00145       M.ReSize(MX);
00146       DCR( D.nric(), C.nric(), 5, R.nric(), 9, M.nric() );
00147       M -= MX;  Print(M);
00148       
00149       // test swap
00150       nricMatrix A(3,4); nricMatrix B(4,5);
00151       A.Row(1) << 2 << 7 << 3 << 6;
00152       A.Row(2) << 6 << 2 << 5 << 9;
00153       A.Row(3) << 1 << 0 << 1 << 6;
00154       B.Row(1) << 2 << 8 << 4 << 5 << 3;
00155       B.Row(2) << 1 << 7 << 5 << 3 << 9;
00156       B.Row(3) << 7 << 8 << 2 << 1 << 6;
00157       B.Row(4) << 5 << 2 << 9 << 0 << 9;
00158       nricMatrix A1(1,2); nricMatrix B1;
00159       nricMatrix X(3,5); Matrix X1 = A * B;
00160       swap(A, A1); swap(B1, B);
00161       for (int i = 1; i <= 3; ++i) for (int j = 1; j <= 5; ++j)
00162       {
00163          X.nric()[i][j] = 0.0;
00164          for (int k = 1; k <= 4; ++k)
00165             X.nric()[i][j] += A1.nric()[i][k] * B1.nric()[k][j];
00166       }
00167       X1 -= X; Print(X1); 
00168    }
00169 #endif
00170    {
00171       Tracer et1("Stage 6");
00172       // test dotproduct
00173       DiagonalMatrix test(5); test = 1;
00174       ColumnVector C(10);
00175       C << 3 << 7 << 5 << 1 << 4 << 2 << 3 << 9 << 1 << 3;
00176       RowVector R(10);
00177       R << 2 << 3 << 5 << 7 << 11 << 13 << -3 << -4 << 2 << 4;
00178       test(1) = (R * C).AsScalar() - DotProduct(C, R);
00179       test(2) = C.SumSquare() - DotProduct(C, C);
00180       test(3) = 6.0 * (C.t() * R.t()).AsScalar() - DotProduct(2.0 * C, 3.0 * R);
00181       Matrix MC = C.AsMatrix(2,5), MR = R.AsMatrix(5,2);
00182       test(4) = DotProduct(MC, MR) - (R * C).AsScalar();
00183       UpperTriangularMatrix UT(5);
00184       UT << 3 << 5 << 2 << 1 << 7
00185               << 1 << 1 << 8 << 2
00186                    << 7 << 0 << 1
00187                         << 3 << 5
00188                              << 6;
00189       LowerTriangularMatrix LT(5);
00190       LT << 5
00191          << 2 << 3
00192          << 1 << 0 << 7
00193          << 9 << 8 << 1 << 2
00194          << 0 << 2 << 1 << 9 << 2;
00195       test(5) = DotProduct(UT, LT) - Sum(SP(UT, LT));
00196       Print(test);
00197       // check row-wise load;
00198       LowerTriangularMatrix LT1(5);
00199       LT1.Row(1) << 5;
00200       LT1.Row(2) << 2   << 3;
00201       LT1.Row(3) << 1   << 0   << 7;
00202       LT1.Row(4) << 9   << 8   << 1   << 2;
00203       LT1.Row(5) << 0   << 2   << 1   << 9   << 2;
00204       Matrix M = LT1 - LT; Print(M);
00205       // check solution with identity matrix
00206       IdentityMatrix IM(5); IM *= 2;
00207       LinearEquationSolver LES1(IM);
00208       LowerTriangularMatrix LTX = LES1.i() * LT;
00209       M = LTX * 2 - LT; Print(M);
00210       DiagonalMatrix D = IM;
00211       LinearEquationSolver LES2(IM);
00212       LTX = LES2.i() * LT;
00213       M = LTX * 2 - LT; Print(M);
00214       UpperTriangularMatrix UTX = LES1.i() * UT;
00215       M = UTX * 2 - UT; Print(M);
00216       UTX = LES2.i() * UT;
00217       M = UTX * 2 - UT; Print(M);
00218    }
00219 
00220    {
00221       Tracer et1("Stage 7");
00222       // Some more GenericMatrix stuff with *= |= &=
00223       // but don't any additional checks
00224       BandMatrix BM1(6,2,3);
00225       BM1.Row(1) << 3 << 8 << 4 << 1;
00226       BM1.Row(2) << 5 << 1 << 9 << 7 << 2;
00227       BM1.Row(3) << 1 << 0 << 6 << 3 << 1 << 3;
00228       BM1.Row(4)      << 4 << 2 << 5 << 2 << 4;
00229       BM1.Row(5)           << 3 << 3 << 9 << 1;
00230       BM1.Row(6)                << 4 << 2 << 9;
00231       BandMatrix BM2(6,1,1);
00232       BM2.Row(1) << 2.5 << 7.5;
00233       BM2.Row(2) << 1.5 << 3.0 << 8.5;
00234       BM2.Row(3)        << 6.0 << 6.5 << 7.0;
00235       BM2.Row(4)               << 2.5 << 2.0 << 8.0;
00236       BM2.Row(5)                      << 0.5 << 4.5 << 3.5;
00237       BM2.Row(6)                             << 9.5 << 7.5;
00238       Matrix RM1 = BM1, RM2 = BM2;
00239       Matrix X;
00240       GenericMatrix GRM1 = RM1, GBM1 = BM1, GRM2 = RM2, GBM2 = BM2;
00241       Matrix Z(6,0); Z = 5; Print(Z);
00242       GRM1 |= Z; GBM1 |= Z; GRM2 &= Z.t(); GBM2 &= Z.t();
00243       X = GRM1 - BM1; Print(X); X = GBM1 - BM1; Print(X);
00244       X = GRM2 - BM2; Print(X); X = GBM2 - BM2; Print(X);
00245 
00246       GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
00247       GRM1 *= GRM2; GBM1 *= GBM2;
00248       X = GRM1 - BM1 * BM2; Print(X);
00249       X = RM1 * RM2 - GBM1; Print(X);
00250 
00251       GRM1 = RM1; GBM1 = BM1; GRM2 = RM2; GBM2 = BM2;
00252       GRM1 *= GBM2; GBM1 *= GRM2;          // Bs and Rs swapped on LHS
00253       X = GRM1 - BM1 * BM2; Print(X);
00254       X = RM1 * RM2 - GBM1; Print(X);
00255 
00256       X = BM1.t(); BandMatrix BM1X = BM1.t();
00257       GRM1 = RM1; X -= GRM1.t(); Print(X); X = BM1X - BM1.t(); Print(X);
00258 
00259       // check that linear equation solver works with Identity Matrix
00260       IdentityMatrix IM(6); IM *= 2;
00261       GBM1 = BM1; GBM1 *= 4; GRM1 = RM1; GRM1 *= 4;
00262       DiagonalMatrix D = IM;
00263       LinearEquationSolver LES1(D);
00264       BandMatrix BX;
00265       BX = LES1.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00266       LinearEquationSolver LES2(IM);
00267       BX = LES2.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00268       BX = D.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00269       BX = IM.i() * GBM1; BX -= BM1 * 2; X = BX; Print(X);
00270       BX = IM.i(); BX *= GBM1; BX -= BM1 * 2; X = BX; Print(X);
00271 
00272       // try symmetric band matrices
00273       SymmetricBandMatrix SBM; SBM << SP(BM1, BM1.t());
00274       SBM << IM.i() * SBM;
00275       X = 2 * SBM - SP(RM1, RM1.t()); Print(X);
00276 
00277       // Do this again with more general D
00278       D << 2.5 << 7.5 << 2 << 5 << 4.5 << 7.5;
00279       BX = D.i() * BM1; X = BX - D.i() * RM1;
00280       Clean(X,0.00000001); Print(X);
00281       BX = D.i(); BX *= BM1; X = BX - D.i() * RM1;
00282       Clean(X,0.00000001); Print(X);
00283       SBM << SP(BM1, BM1.t());
00284       BX = D.i() * SBM; X = BX - D.i() * SP(RM1, RM1.t());
00285       Clean(X,0.00000001); Print(X);
00286 
00287       // test return
00288       BX = TestReturn(BM1); X = BX - BM1;
00289       if (BX.BandWidth() != BM1.BandWidth()) X = 5;
00290       Print(X);
00291    }
00292 
00293 //   cout << "\nEnd of eighth test\n";
00294 }
00295 
00296 
00297 


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