tmt3.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 "newmat.h"
00013 
00014 #include "tmt.h"
00015 
00016 #ifdef use_namespace
00017 using namespace NEWMAT;
00018 #endif
00019 
00020 
00021 /**************************** test program ******************************/
00022 
00023 void trymat3()
00024 {
00025    Tracer et("Third test of Matrix package");
00026    Tracer::PrintTrace();
00027 
00028    {
00029       Tracer et1("Stage 1");
00030       int i,j;
00031       SymmetricMatrix S(7);
00032       for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
00033                 S=-S+2.0;
00034 
00035       DiagonalMatrix D(7);
00036       for (i=1;i<=7;i++) D(i,i)=S(i,i);
00037 
00038       Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0;  M4=M4-4.0; Print(M4); }
00039       SymmetricMatrix S2=D; Matrix M2=S2;  { M2=-D+M2; Print(M2); }
00040       UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
00041       LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
00042       M2=D; M2=M2-D; Print(M2);
00043       for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
00044       U2=L2.t(); D=D.t(); S=S.t();
00045       M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
00046       M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
00047    }
00048 
00049    {
00050       Tracer et1("Stage 2");
00051       int i,j;
00052       DiagonalMatrix D(6);
00053       for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
00054       UpperTriangularMatrix U2(7); LowerTriangularMatrix L2(7);
00055       for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
00056                 { U2=L2.t(); }
00057       DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
00058       Matrix M2(6,7);
00059       for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
00060       Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
00061       Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
00062       MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
00063       {
00064          D.ReSize(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
00065          LowerTriangularMatrix LD(1); LD=D;
00066          UpperTriangularMatrix UD(1); UD=D;
00067          M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
00068          M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
00069          M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
00070          M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
00071       }
00072    }
00073 
00074    {
00075       Tracer et1("Stage 3");
00076       // test inverse * scalar
00077       DiagonalMatrix D(6);
00078       for (int i=1;i<=6;i++) D(i)=i*i;
00079       DiagonalMatrix E = D.i() * 4.0;
00080       DiagonalMatrix I(6); I = 1.0;
00081       E=D*E-I*4.0; Print(E);
00082       E = D.i() / 0.25; E=D*E-I*4.0; Print(E);
00083    }
00084    {
00085       Tracer et1("Stage 4");
00086       Matrix sigma(3,3); Matrix sigmaI(3,3);
00087       sigma = 0; sigma(1,1) = 1.0; sigma(2,2) = 1.0; sigma(3,3) = 1.0;
00088       sigmaI = sigma.i();
00089       sigmaI -= sigma;  Clean(sigmaI, 0.000000001); Print(sigmaI);
00090    }
00091    {
00092       Tracer et1("Stage 5");
00093       Matrix X(5,5); DiagonalMatrix DM(5);
00094       for (int i=1; i<=5; i++) for (int j=1; j<=5; j++)
00095          X(i,j) = (23*i+59*j) % 43;
00096       DM << 1 << 8 << -7 << 2 << 3;
00097       Matrix Y = X.i() * DM; Y = X * Y - DM;
00098       Clean(Y, 0.000000001); Print(Y);
00099    }
00100    {
00101       Tracer et1("Stage 6");          // test reverse function
00102       ColumnVector CV(10), RCV(10);
00103       CV  <<  2 << 7 << 1  << 6 << -3 <<  1 << 8 << -4 << 0 << 17;
00104       RCV << 17 << 0 << -4 << 8 << 1  << -3 << 6 <<  1 << 7 << 2;
00105       ColumnVector X = CV - RCV.Reverse(); Print(X);
00106       RowVector Y = CV.t() - RCV.t().Reverse(); Print(Y);
00107       DiagonalMatrix D = CV.AsDiagonal() - RCV.AsDiagonal().Reverse();
00108       Print(D);
00109       X = CV & CV.Rows(1,9).Reverse();
00110       ColumnVector Z(19);
00111       Z.Rows(1,10) = RCV.Reverse(); Z.Rows(11,19) = RCV.Rows(2,10);
00112       X -= Z; Print(X); Z -= Z.Reverse(); Print(Z);
00113       Matrix A(3,3); A << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9;
00114       Matrix B(3,3); B << 9 << 8 << 7 << 6 << 5 << 4 << 3 << 2 << 1;
00115       Matrix Diff = A - B.Reverse(); Print(Diff);
00116       Diff = (-A).Reverse() + B; Print(Diff);
00117       UpperTriangularMatrix U;
00118       U << A.Reverse(); Diff = U; U << B; Diff -= U; Print(Diff);
00119       U << (-A).Reverse(); Diff = U; U << B; Diff += U; Print(Diff);
00120    }
00121    {
00122       Tracer et1("Stage 7");           // test IsSingular function
00123       ColumnVector XX(4);
00124       Matrix A(3,3);
00125       A = 0;
00126       CroutMatrix B1 = A;
00127       XX(1) = B1.IsSingular() ? 0 : 1;
00128       A << 1 << 3 << 6
00129         << 7 << 11 << 13
00130         << 2 << 4  << 1;
00131       CroutMatrix B2(A);
00132       XX(2) = B2.IsSingular() ? 1 : 0;
00133       BandMatrix C(3,1,1); C.Inject(A);
00134       BandLUMatrix B3(C);
00135       XX(3) = B3.IsSingular() ? 1 : 0;
00136       C = 0;
00137       BandLUMatrix B4(C);
00138       XX(4) = B4.IsSingular() ? 0 : 1;
00139       Print(XX);
00140    }
00141    {
00142       Tracer et1("Stage 8");           // inverse with vector of 0s
00143       Matrix A(3,3); Matrix Z(3,3); ColumnVector X(6);
00144       A <<  1 <<  3 <<  6
00145         <<  7 << 11 << 13
00146         <<  2 <<  4 <<  1;
00147       Z = 0;
00148       Matrix B = (A | Z) & (Z | A);   // 6 * 6 matrix
00149       X = 0.0;
00150       X = B.i() * X;
00151       Print(X);
00152       // also check inverse with non-zero Y
00153       Matrix Y(3,3);
00154       Y << 0.0 << 1.0 << 1.0
00155         << 5.0 << 0.0 << 5.0
00156         << 3.0 << 3.0 << 0.0;
00157       Matrix YY = Y & Y;        // stack Y matrices
00158       YY = B.i() * YY;
00159       Matrix Y1 = A.i() * Y;
00160       YY -= Y1 & Y1; Clean(YY, 0.000000001); Print(YY);
00161       Y1 = A * Y1 - Y; Clean(Y1, 0.000000001); Print(Y1);
00162    }
00163 
00164 
00165 }
00166 
00167 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13