tmt1.cc
Go to the documentation of this file.
00001 
00002 #define WANT_STREAM
00003 
00004 
00005 
00006 #include "include.h"
00007 
00008 #include "newmat.h"
00009 
00010 #include "tmt.h"
00011 
00012 #ifdef use_namespace
00013 using namespace NEWMAT;
00014 #endif
00015 
00016 
00017 /**************************** test program ******************************/
00018 
00019 
00020 void trymat1()
00021 {
00022 //   cout << "\nFirst test of Matrix package\n\n";
00023    Tracer et("First test of Matrix package");
00024    Tracer::PrintTrace();
00025    {
00026       Tracer et1("Stage 1");
00027       int i,j;
00028 
00029       LowerTriangularMatrix L(10);
00030       for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
00031       SymmetricMatrix S(10);
00032       for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
00033       SymmetricMatrix S1 = S / 2.0;
00034       S = S1 * 2.0;
00035       UpperTriangularMatrix U=L.t()*2.0;
00036       Print(LowerTriangularMatrix(L-U.t()*0.5));
00037       DiagonalMatrix D(10);
00038       for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
00039       Matrix M=(S+U-D+L)*(L+U-D+S);
00040       DiagonalMatrix DD=D*D;
00041       LowerTriangularMatrix LD=L*D;
00042       // expressions split for Turbo C
00043       Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
00044       { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
00045       { M1 = M1 - U*D + DD - LD + S*S; }
00046       { M1 = M1 + U*S - D*S + L*S - 10.0; }
00047       M=M1-M;
00048       Print(M);
00049    }
00050    {
00051       Tracer et1("Stage 2");
00052       int i,j;
00053 
00054       LowerTriangularMatrix L(9);
00055       for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
00056       UpperTriangularMatrix U1(9);
00057       for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
00058       LowerTriangularMatrix LX(9);
00059       for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
00060       UpperTriangularMatrix UX(9);
00061       for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
00062       {
00063          L=L+LX/0.5;   L=L-LX*3.0;   L=LX*2.0+L;
00064          U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
00065       }
00066 
00067 
00068       SymmetricMatrix S(9);
00069       for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
00070       {
00071          SymmetricMatrix S1 = S;
00072          S=S1+5.0;
00073          S=S-3.0;
00074       }
00075 
00076       DiagonalMatrix D(9);
00077       for (i=1;i<=9;i++) D(i,i)=S(i,i);
00078       UpperTriangularMatrix U=L.t()*2.0;
00079       {
00080          U1=U1*2.0 - U;  Print(U1);
00081          L=L*2.0-D; U=U-D;
00082       }
00083       Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
00084    }
00085    {
00086       Tracer et1("Stage 3");
00087       int i,j;
00088       Matrix M(10,3), N(10,3);
00089       for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
00090          {  M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00091       Matrix MN = M + N, M1;
00092 
00093       M1 = M; M1 += N; M1 -= MN; Print(M1);
00094       M1 = M; M1 += M1; M1 = M1 - M * 2; Print(M1);
00095       M1 = M; M1 += N * 2; M1 -= (MN + N); Print(M1);
00096       M1 = M; M1 -= M1; Print(M1);
00097       M1 = M; M1 -= MN + M1; M1 += N + M; Print(M1);
00098       M1 = M; M1 -= 5; M1 -= M; M1 *= 0.2; M1 = M1 + 1; Print(M1);
00099       Matrix NT = N.t();
00100       M1 = M; M1 *= NT; M1 -= M * N.t(); Print(M1);
00101       M = M * M.t();
00102       DiagonalMatrix D(10); D = 2;
00103       M1 = M; M1 += D; M1 -= M; M1 = M1 - D; Print(M1);
00104       M1 = M; M1 -= D; M1 -= M; M1 = M1 + D; Print(M1);
00105       M1 = M; M1 *= D; M1 /= 2; M1 -= M; Print(M1);
00106       SymmetricMatrix SM; SM << M;
00107       // UpperTriangularMatrix SM; SM << M;
00108       SM += 10; M1 = SM - M; M1 /=10; M1 = M1 - 1; Print(M1);
00109    }
00110    {
00111       Tracer et1("Stage 4");
00112       int i,j;
00113       Matrix M(10,3), N(10,5);
00114       for (i = 1; i<=10; i++) for (j = 1; j<=3; j++) M(i,j) = 2*i-j;
00115       for (i = 1; i<=10; i++) for (j = 1; j<=5; j++) N(i,j) = i*j + 20;
00116       Matrix M1;
00117       M1 = M; M1 |= N; M1 &= N | M;
00118       M1 -= (M | N) & (N | M); Print(M1);
00119       M1 = M; M1 |= M1; M1 &= M1;
00120       M1 -= (M | M) & (M | M); Print(M1);
00121 
00122    }
00123    {
00124       Tracer et1("Stage 5");
00125       int i,j;
00126       BandMatrix BM1(10,2,3), BM2(10,4,1); Matrix M1(10,10), M2(10,10);
00127       for (i=1;i<=10;i++) for (j=1;j<=10;j++)
00128         { M1(i,j) = 0.5*i+j*j-50; M2(i,j) = (i*101 + j*103) % 13; }
00129       BM1.Inject(M1); BM2.Inject(M2);
00130       BandMatrix BM = BM1; BM += BM2;
00131       Matrix M1X = BM1; Matrix M2X = BM2; Matrix MX = BM;
00132       MX -= M1X + M2X; Print(MX);
00133       MX = BM1; MX += BM2; MX -= M1X; MX -= M2X; Print(MX);
00134       SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t(); 
00135       SymmetricBandMatrix SM2; SM2 << BM2 * BM2.t();
00136       SM1 *= 5.5;
00137       M1X *= M1X.t(); M1X *= 5.5; M2X *= M2X.t();
00138       SM1 -= SM2; M1 = SM1 - M1X + M2X; Print(M1);
00139       M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1); 
00140       M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1); 
00141       M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1); 
00142       
00143    }
00144    {
00145       Tracer et1("Stage 6");
00146       int i,j;
00147       Matrix M(10,10), N(10,10);
00148       for (i = 1; i<=10; i++) for (j = 1; j<=10; j++)
00149          {  M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00150       GenericMatrix GM = M;
00151       GM += N; Matrix M1 = GM - N - M; Print(M1);
00152       DiagonalMatrix D(10); D = 3;
00153       GM = D; GM += N; GM += M; GM += D;
00154       M1 = D*2 - GM + M + N; Print(M1);
00155       GM = D; GM *= 4; GM += 16; GM /= 8; GM -= 2;
00156       GM -= D / 2; M1 = GM; Print(M1);
00157       GM = D; GM *= M; GM *= N; GM /= 3; M1 = M*N - GM; Print(M1);
00158       GM = D; GM |= M; GM &= N | D; M1 = GM - ((D | M) & (N | D));
00159       Print(M1);
00160       GM = M; M1 = M; GM += 5; GM *= 3; M *= 3; M += 15; M1 = GM - M;
00161       Print(M1);
00162       D.ReSize(10); for (i = 1; i<=10; i++) D(i) = i;
00163       M1 = D + 10; GM = D; GM += 10; M1 -= GM; Print(M1);
00164       GM = M; GM -= D; M1 = GM; GM = D; GM -= M; M1 += GM; Print(M1);
00165       GM = M; GM *= D; M1 = GM; GM = D; GM *= M.t();
00166       M1 -= GM.t(); Print(M1);
00167       GM = M; GM += 2 * GM; GM -= 3 * M; M1 = GM; Print(M1);
00168       GM = M; GM |= GM; GM -= (M | M); M1 = GM; Print(M1);
00169       GM = M; GM &= GM; GM -= (M & M); M1 = GM; Print(M1);
00170       M1 = M; M1 = (M1.t() & M.t()) - (M | M).t(); Print(M1);
00171       M1 = M; M1 = (M1.t() | M.t()) - (M & M).t(); Print(M1);
00172 
00173    }
00174 
00175    {
00176       Tracer et1("Stage 7");
00177       // test for bug in MS VC5
00178       int n = 3;
00179       int k; int j;
00180       Matrix A(n,n), B(n,n);
00181 
00182       //first version - MS VC++ 5 mis-compiles if optimisation is on
00183       for (k=1; k<=n; k++)
00184       {
00185          for (j = 1; j <= n; j++) A(k,j) = ((k-1) * (2*j-1));
00186       }
00187 
00188       //second version
00189       for (k=1; k<=n; k++)
00190       {
00191          const int k1 = k-1;          // otherwise Visual C++ 5 fails
00192          for (j = 1; j <= n; j++) B(k,j) = (k1 * (2*j-1));
00193       }
00194 
00195       if (A != B)
00196       {
00197          cout << "\nVisual C++ version 5 compiler error?";
00198          cout << "\nTurn off optimisation";
00199       }
00200 
00201       A -= B; Print(A);
00202 
00203    }
00204 
00205 //   cout << "\nEnd of first test\n";
00206 }
00207 


rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:13