tmt6.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 #define WANT_MATH
00010 
00011 
00012 #include "include.h"
00013 
00014 #include "newmatap.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 // slow sort program
00027 
00028 static void SimpleSortDescending(Real* first, const int length)
00029 {
00030    int i = length;
00031    while (--i)
00032    {
00033       Real x = *first; Real* f = first; Real* g = f;
00034       int j = i;
00035       while (j--) if (x < *(++f)) { g = f; x = *g; }
00036       *g = *first; *first++ = x;
00037    }
00038 }
00039 
00040 static void TestSort(int n)
00041 {
00042    // make some data
00043    RowVector X(n);
00044    int i;
00045    for (i = 1; i <= n; i++)
00046       X(i) = sin((Real)i) + 0.3 * cos(i/5.0) - 0.6 * sin(i/7.0) + 0.2 * sin(2.0 * i);
00047    RowVector X_Sorted = X; SimpleSortDescending(X_Sorted.Store(), n);
00048    RowVector A = X + X.Reverse(); SimpleSortDescending(A.Store(), n);
00049 
00050    // test descending sort
00051 
00052    RowVector Y = X; SortDescending(Y); Y -= X_Sorted; Print(Y);
00053    Y = X_Sorted; SortDescending(Y); Y -= X_Sorted; Print(Y);
00054    Y = X_Sorted.Reverse(); SortDescending(Y); Y -= X_Sorted; Print(Y);
00055    Y = X + X.Reverse(); SortDescending(Y); Y -= A; Print(Y);
00056 
00057    // test ascending sort
00058 
00059    Y = X; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
00060    Y = X_Sorted; SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
00061    Y = X_Sorted.Reverse(); SortAscending(Y); Y -= X_Sorted.Reverse(); Print(Y);
00062    Y = X + X.Reverse(); SortAscending(Y); Y -= A.Reverse(); Print(Y);
00063 }
00064 
00065 
00066 void trymat6()
00067 {
00068    Tracer et("Sixth test of Matrix package");
00069    Tracer::PrintTrace();
00070 
00071    int i,j;
00072 
00073 
00074    DiagonalMatrix D(6);
00075    UpperTriangularMatrix U(6);
00076    for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; }
00077    LowerTriangularMatrix L=(U*3.0).t();
00078    SymmetricMatrix S(6);
00079    for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
00080    Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S;
00081    Matrix M(6,6);
00082    for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
00083    {
00084       Tracer et1("Stage 1");
00085       Print(Matrix(MS+(-MS)));
00086       Print(Matrix((S+M)-(MS+M)));
00087       Print(Matrix((M+U)-(M+MU)));
00088       Print(Matrix((M+L)-(M+ML)));
00089    }
00090    {
00091       Tracer et1("Stage 2");
00092       Print(Matrix((M+D)-(M+MD)));
00093       Print(Matrix((U+D)-(MU+MD)));
00094       Print(Matrix((D+L)-(ML+MD)));
00095       Print(Matrix((-U+D)+MU-MD));
00096       Print(Matrix((-L+D)+ML-MD));
00097    }
00098    {
00099       Tracer et1("Stage 3 - concatenate");
00100       RowVector A(5);
00101       A << 1 << 2 << 3 << 4 << 5;
00102       RowVector B(5);
00103       B << 3 << 1 << 4 << 1 << 5;
00104       Matrix C(3,5);
00105       C <<  2 <<  3 <<  5 <<  7 << 11
00106         << 13 << 17 << 19 << 23 << 29
00107         << 31 << 37 << 41 << 43 << 47;
00108       Matrix X1 = A & B & C;
00109       Matrix X2 = (A.t() | B.t() | C.t()).t();
00110       Matrix X3(5,5);
00111       X3.Row(1)=A; X3.Row(2)=B; X3.Rows(3,5)=C;
00112       Print(Matrix(X1-X2));
00113       Print(Matrix(X1-X3));
00114       LowerTriangularMatrix LT1; LT1 << (A & B & C);
00115       UpperTriangularMatrix UT1; UT1 << (A.t() | B.t() | C.t());
00116       Print(LowerTriangularMatrix(LT1-UT1.t()));
00117       DiagonalMatrix D1; D1 << (A.t() | B.t() | C.t());
00118       ColumnVector At = A.t();
00119       ColumnVector Bt = B.t();
00120       Matrix Ct = C.t();
00121       LowerTriangularMatrix LT2; LT2 << (At | Bt | Ct);
00122       UpperTriangularMatrix UT2; UT2 << (At.t() & Bt.t() & Ct.t());
00123       Matrix ABt = At | Bt;
00124       DiagonalMatrix D2; D2 << (ABt | Ct);
00125       Print(LowerTriangularMatrix(LT2-UT2.t()));
00126       Print(DiagonalMatrix(D1-D2));
00127       Print(Matrix(LT1+UT2-D2-X1));
00128       Matrix M1 = LT1 | UT2; Matrix M2 = UT1 & LT2;
00129       Print(Matrix(M1-M2.t()));
00130       M1 = UT2 | LT1; M2 = LT2 & UT1;
00131       Print(Matrix(M1-M2.t()));
00132       M1 = (LT1 | UT2) & (UT2 | LT1);
00133       M2 = (UT1 & LT2) | (LT2 & UT1);
00134       Print(Matrix(M1-M2.t()));
00135       SymmetricMatrix SM1; SM1 << (M1 + M1.t());
00136       SymmetricMatrix SM2; SM2 << ((SM1 | M1) & (M1.t() | SM1));
00137       Matrix M3(20,20);
00138       M3.SubMatrix(1,10,1,10) = SM1;
00139       M3.SubMatrix(1,10,11,20) = M1;
00140       M3.SubMatrix(11,20,1,10) = M2;
00141       M3.SubMatrix(11,20,11,20) = SM1;
00142       Print(Matrix(M3-SM2));
00143 
00144       SymmetricMatrix SM(15); SM = 0; SM.SymSubMatrix(1,10) = SM1;
00145       M3.ReSize(15,15); M3 = 0; M3.SubMatrix(1,10,1,10) = SM1;
00146       M3 -= SM; Print(M3);
00147       SM = 0; SM.SymSubMatrix(6,15) = SM1;
00148       M3.ReSize(15,15); M3 = 0; M3.SubMatrix(6,15,6,15) = SM1;
00149       M3 = M3.t() - SM; Print(M3);
00150    }
00151    {
00152       Tracer et1("Stage 4 - sort");
00153       TestSort(1); TestSort(2); TestSort(3); TestSort(4);
00154       TestSort(15); TestSort(16); TestSort(17); TestSort(18);
00155       TestSort(99); TestSort(100); TestSort(101);
00156    }
00157 
00158 
00159 //   cout << "\nEnd of sixth test\n";
00160 }
00161 
00162 


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