tmt1.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 #define WANT_STREAM
00009 
00010 
00011 
00012 #include "include.h"
00013 
00014 #include "newmat.h"
00015 
00016 #include "tmt.h"
00017 
00018 #ifdef use_namespace
00019 using namespace NEWMAT;
00020 #endif
00021 
00022 
00023 /**************************** test program ******************************/
00024 
00025 
00026 void trymat1()
00027 {
00028 //   cout << "\nFirst test of Matrix package\n\n";
00029    Tracer et("First test of Matrix package");
00030    Tracer::PrintTrace();
00031    {
00032       Tracer et1("Stage 1");
00033       int i,j;
00034 
00035       LowerTriangularMatrix L(10);
00036       for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
00037       SymmetricMatrix S(10);
00038       for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
00039       SymmetricMatrix S1 = S / 2.0;
00040       S = S1 * 2.0;
00041       UpperTriangularMatrix U=L.t()*2.0;
00042       Print(LowerTriangularMatrix(L-U.t()*0.5));
00043       DiagonalMatrix D(10);
00044       for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
00045       Matrix M=(S+U-D+L)*(L+U-D+S);
00046       DiagonalMatrix DD=D*D;
00047       LowerTriangularMatrix LD=L*D;
00048       // expressions split for Turbo C
00049       Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
00050       { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
00051       { M1 = M1 - U*D + DD - LD + S*S; }
00052       { M1 = M1 + U*S - D*S + L*S - 10.0; }
00053       M=M1-M;
00054       Print(M);
00055    }
00056    {
00057       Tracer et1("Stage 2");
00058       int i,j;
00059 
00060       LowerTriangularMatrix L(9);
00061       for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
00062       UpperTriangularMatrix U1(9);
00063       for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
00064       LowerTriangularMatrix LX(9);
00065       for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
00066       UpperTriangularMatrix UX(9);
00067       for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
00068       {
00069          L=L+LX/0.5;   L=L-LX*3.0;   L=LX*2.0+L;
00070          U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
00071       }
00072 
00073 
00074       SymmetricMatrix S(9);
00075       for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
00076       {
00077          SymmetricMatrix S1 = S;
00078          S=S1+5.0;
00079          S=S-3.0;
00080       }
00081 
00082       DiagonalMatrix D(9);
00083       for (i=1;i<=9;i++) D(i,i)=S(i,i);
00084       UpperTriangularMatrix U=L.t()*2.0;
00085       {
00086          U1=U1*2.0 - U;  Print(U1);
00087          L=L*2.0-D; U=U-D;
00088       }
00089       Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
00090    }
00091    {
00092       Tracer et1("Stage 3");
00093       int i,j;
00094       Matrix M(10,3), N(10,3);
00095       for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
00096          {  M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00097       Matrix MN = M + N, M1;
00098 
00099       M1 = M; M1 += N; M1 -= MN; Print(M1);
00100       M1 = M; M1 += M1; M1 = M1 - M * 2; Print(M1);
00101       M1 = M; M1 += N * 2; M1 -= (MN + N); Print(M1);
00102       M1 = M; M1 -= M1; Print(M1);
00103       M1 = M; M1 -= MN + M1; M1 += N + M; Print(M1);
00104       M1 = M; M1 -= 5; M1 -= M; M1 *= 0.2; M1 = M1 + 1; Print(M1);
00105       Matrix NT = N.t();
00106       M1 = M; M1 *= NT; M1 -= M * N.t(); Print(M1);
00107       M = M * M.t();
00108       DiagonalMatrix D(10); D = 2;
00109       M1 = M; M1 += D; M1 -= M; M1 = M1 - D; Print(M1);
00110       M1 = M; M1 -= D; M1 -= M; M1 = M1 + D; Print(M1);
00111       M1 = M; M1 *= D; M1 /= 2; M1 -= M; Print(M1);
00112       SymmetricMatrix SM; SM << M;
00113       // UpperTriangularMatrix SM; SM << M;
00114       SM += 10; M1 = SM - M; M1 /=10; M1 = M1 - 1; Print(M1);
00115    }
00116    {
00117       Tracer et1("Stage 4");
00118       int i,j;
00119       Matrix M(10,3), N(10,5);
00120       for (i = 1; i<=10; i++) for (j = 1; j<=3; j++) M(i,j) = 2*i-j;
00121       for (i = 1; i<=10; i++) for (j = 1; j<=5; j++) N(i,j) = i*j + 20;
00122       Matrix M1;
00123       M1 = M; M1 |= N; M1 &= N | M;
00124       M1 -= (M | N) & (N | M); Print(M1);
00125       M1 = M; M1 |= M1; M1 &= M1;
00126       M1 -= (M | M) & (M | M); Print(M1);
00127 
00128    }
00129    {
00130       Tracer et1("Stage 5");
00131       int i,j;
00132       BandMatrix BM1(10,2,3), BM2(10,4,1); Matrix M1(10,10), M2(10,10);
00133       for (i=1;i<=10;i++) for (j=1;j<=10;j++)
00134         { M1(i,j) = 0.5*i+j*j-50; M2(i,j) = (i*101 + j*103) % 13; }
00135       BM1.Inject(M1); BM2.Inject(M2);
00136       BandMatrix BM = BM1; BM += BM2;
00137       Matrix M1X = BM1; Matrix M2X = BM2; Matrix MX = BM;
00138       MX -= M1X + M2X; Print(MX);
00139       MX = BM1; MX += BM2; MX -= M1X; MX -= M2X; Print(MX);
00140       SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t(); 
00141       SymmetricBandMatrix SM2; SM2 << BM2 * BM2.t();
00142       SM1 *= 5.5;
00143       M1X *= M1X.t(); M1X *= 5.5; M2X *= M2X.t();
00144       SM1 -= SM2; M1 = SM1 - M1X + M2X; Print(M1);
00145       M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1); 
00146       M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1); 
00147       M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1); 
00148       
00149    }
00150    {
00151       Tracer et1("Stage 6");
00152       int i,j;
00153       Matrix M(10,10), N(10,10);
00154       for (i = 1; i<=10; i++) for (j = 1; j<=10; j++)
00155          {  M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00156       GenericMatrix GM = M;
00157       GM += N; Matrix M1 = GM - N - M; Print(M1);
00158       DiagonalMatrix D(10); D = 3;
00159       GM = D; GM += N; GM += M; GM += D;
00160       M1 = D*2 - GM + M + N; Print(M1);
00161       GM = D; GM *= 4; GM += 16; GM /= 8; GM -= 2;
00162       GM -= D / 2; M1 = GM; Print(M1);
00163       GM = D; GM *= M; GM *= N; GM /= 3; M1 = M*N - GM; Print(M1);
00164       GM = D; GM |= M; GM &= N | D; M1 = GM - ((D | M) & (N | D));
00165       Print(M1);
00166       GM = M; M1 = M; GM += 5; GM *= 3; M *= 3; M += 15; M1 = GM - M;
00167       Print(M1);
00168       D.ReSize(10); for (i = 1; i<=10; i++) D(i) = i;
00169       M1 = D + 10; GM = D; GM += 10; M1 -= GM; Print(M1);
00170       GM = M; GM -= D; M1 = GM; GM = D; GM -= M; M1 += GM; Print(M1);
00171       GM = M; GM *= D; M1 = GM; GM = D; GM *= M.t();
00172       M1 -= GM.t(); Print(M1);
00173       GM = M; GM += 2 * GM; GM -= 3 * M; M1 = GM; Print(M1);
00174       GM = M; GM |= GM; GM -= (M | M); M1 = GM; Print(M1);
00175       GM = M; GM &= GM; GM -= (M & M); M1 = GM; Print(M1);
00176       M1 = M; M1 = (M1.t() & M.t()) - (M | M).t(); Print(M1);
00177       M1 = M; M1 = (M1.t() | M.t()) - (M & M).t(); Print(M1);
00178 
00179    }
00180 
00181    {
00182       Tracer et1("Stage 7");
00183       // test for bug in MS VC5
00184       int n = 3;
00185       int k; int j;
00186       Matrix A(n,n), B(n,n);
00187 
00188       //first version - MS VC++ 5 mis-compiles if optimisation is on
00189       for (k=1; k<=n; k++)
00190       {
00191          for (j = 1; j <= n; j++) A(k,j) = ((k-1) * (2*j-1));
00192       }
00193 
00194       //second version
00195       for (k=1; k<=n; k++)
00196       {
00197          const int k1 = k-1;          // otherwise Visual C++ 5 fails
00198          for (j = 1; j <= n; j++) B(k,j) = (k1 * (2*j-1));
00199       }
00200 
00201       if (A != B)
00202       {
00203          cout << "\nVisual C++ version 5 compiler error?";
00204          cout << "\nTurn off optimisation";
00205       }
00206 
00207       A -= B; Print(A);
00208 
00209    }
00210 
00211    {
00212       Tracer et1("Stage 8");
00213       // Cross product
00214       ColumnVector i(3); i << 1 << 0 << 0;
00215       ColumnVector j(3); j << 0 << 1 << 0;
00216       ColumnVector k(3); k << 0 << 0 << 1;
00217       ColumnVector X;
00218       X = CrossProduct(i,j) - k; Print(X);
00219       X = CrossProduct(j,k) - i; Print(X);
00220       X = CrossProduct(k,i) - j; Print(X);
00221       X = CrossProduct(j,i) + k; Print(X);
00222       X = CrossProduct(k,j) + i; Print(X);
00223       X = CrossProduct(i,k) + j; Print(X);
00224       X = CrossProduct(i,i); Print(X);
00225       X = CrossProduct(j,j); Print(X);
00226       X = CrossProduct(k,k); Print(X);
00227 
00228       ColumnVector A(3); A << 2.25 << 1.75 << -7.5;
00229       ColumnVector B(3); B << -0.5 << 4.75 << 3.25;
00230       ColumnVector C(3); C << 9.25 << 3.5  << 1.25;
00231 
00232       Real d0 = (A | B | C).Determinant();    // Vector triple product
00233       Real d1 = DotProduct(CrossProduct(A, B), C);
00234       Real d2 = DotProduct(CrossProduct(B, C), A);
00235       Real d3 = DotProduct(CrossProduct(C, A), B);
00236       X << (d1 - d0) << (d2 - d0) << (d3 - d0);
00237       Clean(X, 0.000000001); Print(X);
00238 
00239       X = CrossProduct(A, CrossProduct(B, C))
00240         + CrossProduct(B, CrossProduct(C, A))
00241         + CrossProduct(C, CrossProduct(A, B));
00242       Print(X);
00243 
00244       RowVector XT = CrossProduct(A.AsRow(), B.AsRow());
00245       XT -= CrossProduct(A, B).AsRow();
00246       Print(XT);
00247    }
00248 
00249    {
00250       Tracer et1("Stage 9");
00251       // More cross product
00252       int i, j;
00253       Matrix M(10,3), N(10,3);
00254       for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
00255          {  M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00256 
00257       Matrix CP1 = CrossProductRows(M, N);
00258       Matrix CP2(10,3);
00259       for (i = 1; i<=10; i++)
00260          CP2.Row(i) = CrossProduct(M.Row(i), N.Row(i));
00261       CP2 -= CP1; Print(CP2);
00262 
00263       CP2 = CrossProductColumns(M.t(), N.t());
00264       CP2 -= CP1.t(); Print(CP2);
00265    }
00266 
00267    {
00268       Tracer et1("Stage 10");
00269       // Make sure RNG works
00270       MultWithCarry mwc;
00271       ColumnVector cv(10);
00272       for (int i = 1; i <= 10; ++i) cv(i) = mwc.Next();
00273       cv *= 100.0;
00274       cv(1) -= 6.27874; 
00275       cv(2) -= 42.1718; 
00276       cv(3) -= 80.2854; 
00277       cv(4) -= 12.961;  
00278       cv(5) -= 17.7499; 
00279       cv(6) -= 13.2657; 
00280       cv(7) -= 50.4923; 
00281       cv(8) -= 26.095;  
00282       cv(9) -= 57.9147; 
00283       cv(10) -= 30.1778;        
00284       Clean(cv, 0.0001); Print(cv);
00285    }
00286 
00287 
00288 //   cout << "\nEnd of first test\n";
00289 }
00290 
00291 


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