tmt.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 #define WANT_STREAM
00008 #define WANT_TIME
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 namespace NEWMAT {
00019 #endif
00020 
00021 
00022 /**************************** test program ******************************/
00023 
00024 
00025 class PrintCounter
00026 {
00027    int count;
00028    const char* s;
00029 public:
00030    ~PrintCounter();
00031    PrintCounter(const char * sx) : count(0), s(sx) {}
00032    void operator++() { count++; }
00033 };
00034 
00035 PrintCounter PCZ("Number of non-zero matrices (should be 1) = ");
00036 PrintCounter PCN("Number of matrices tested                 = ");
00037 
00038 PrintCounter::~PrintCounter()
00039 { cout << s << count << "\n"; }
00040 
00041 
00042 void Print(const Matrix& X)
00043 {
00044    ++PCN;
00045    cout << "\nMatrix type: " << X.Type().Value() << " (";
00046    cout << X.Nrows() << ", ";
00047    cout << X.Ncols() << ")\n\n";
00048    if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
00049    int nr=X.Nrows(); int nc=X.Ncols();
00050    for (int i=1; i<=nr; i++)
00051    {
00052       for (int j=1; j<=nc; j++)  cout << X(i,j) << "\t";
00053       cout << "\n";
00054    }
00055    cout << flush; ++PCZ;
00056 }
00057 
00058 void Print(const UpperTriangularMatrix& X)
00059 {
00060    ++PCN;
00061    cout << "\nMatrix type: " << X.Type().Value() << " (";
00062    cout << X.Nrows() << ", ";
00063    cout << X.Ncols() << ")\n\n";
00064    if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
00065    int nr=X.Nrows(); int nc=X.Ncols();
00066    for (int i=1; i<=nr; i++)
00067    {
00068       int j;
00069       for (j=1; j<i; j++) cout << "\t";
00070       for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
00071       cout << "\n";
00072    }
00073    cout << flush; ++PCZ;
00074 }
00075 
00076 void Print(const DiagonalMatrix& X)
00077 {
00078    ++PCN;
00079    cout << "\nMatrix type: " << X.Type().Value() << " (";
00080    cout << X.Nrows() << ", ";
00081    cout << X.Ncols() << ")\n\n";
00082    if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
00083    int nr=X.Nrows(); int nc=X.Ncols();
00084    for (int i=1; i<=nr; i++)
00085    {
00086       for (int j=1; j<i; j++) cout << "\t";
00087       if (i<=nc) cout << X(i,i) << "\t";
00088       cout << "\n";
00089    }
00090    cout << flush; ++PCZ;
00091 }
00092 
00093 void Print(const SymmetricMatrix& X)
00094 {
00095    ++PCN;
00096    cout << "\nMatrix type: " << X.Type().Value() << " (";
00097    cout << X.Nrows() << ", ";
00098    cout << X.Ncols() << ")\n\n";
00099    if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
00100    int nr=X.Nrows(); int nc=X.Ncols();
00101    for (int i=1; i<=nr; i++)
00102    {
00103       int j;
00104       for (j=1; j<i; j++) cout << X(j,i) << "\t";
00105       for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
00106       cout << "\n";
00107    }
00108    cout << flush; ++PCZ;
00109 }
00110 
00111 void Print(const LowerTriangularMatrix& X)
00112 {
00113    ++PCN;
00114    cout << "\nMatrix type: " << X.Type().Value() << " (";
00115    cout << X.Nrows() << ", ";
00116    cout << X.Ncols() << ")\n\n";
00117    if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
00118    int nr=X.Nrows();
00119    for (int i=1; i<=nr; i++)
00120    {
00121       for (int j=1; j<=i; j++) cout << X(i,j) << "\t";
00122       cout << "\n";
00123    }
00124    cout << flush; ++PCZ;
00125 }
00126 
00127 
00128 void Clean(Matrix& A, Real c)
00129 {
00130    int nr = A.Nrows(); int nc = A.Ncols();
00131    for (int i=1; i<=nr; i++)
00132    {
00133       for (int j=1; j<=nc; j++)
00134       { Real a = A(i,j); if ((a < c) && (a > -c)) A(i,j) = 0.0; }
00135    }
00136 }
00137 
00138 void Clean(DiagonalMatrix& A, Real c)
00139 {
00140    int nr = A.Nrows();
00141    for (int i=1; i<=nr; i++)
00142    { Real a = A(i,i); if ((a < c) && (a > -c)) A(i,i) = 0.0; }
00143 }
00144 
00145 void PentiumCheck(Real N, Real D)
00146 {
00147    Real R = N / D;
00148    R = R * D - N;
00149    if ( R > 1 || R < -1)
00150       cout << "Pentium error detected: % error = " << 100 * R / N << "\n";
00151 }
00152 
00153 // random number generator class
00154 // See newran03 documentation for details
00155 
00156 MultWithCarry::MultWithCarry(double s)
00157 {
00158    if (s>=1.0 || s<=0.0)
00159       Throw(Logic_error("MultWithCarry: seed out of range"));
00160    x = (unsigned long)(s * 4294967296.0);
00161    crry = 1234567;
00162 }
00163 
00164 
00165 // carry out 32bit * 32bit multiply in software
00166 
00167 void MultWithCarry::NextValue()
00168 {
00169    unsigned long  mult = 2083801278;
00170    unsigned long  m_hi = mult >> 16;
00171    unsigned long  m_lo = mult & 0xFFFF;
00172 
00173    unsigned long  x_hi = x >> 16;
00174    unsigned long  x_lo = x & 0xFFFF;
00175 
00176    unsigned long  c_hi = crry >> 16;
00177    unsigned long  c_lo = crry & 0xFFFF;
00178 
00179    x = x_lo * m_lo + c_lo;
00180    unsigned long axc = x_lo * m_hi + x_hi * m_lo + c_hi + (x >> 16);
00181    crry = x_hi * m_hi + (axc >> 16);
00182 
00183    x = (x & 0xFFFF) + (axc << 16);
00184 
00185 }
00186 
00187 Real MultWithCarry::Next() { NextValue(); return ((double)x + 0.5) / 4294967296.0; }
00188 
00189 // fill a matrix with values from the MultWithCarry random number generator
00190 void FillWithValues(MultWithCarry&MWC, Matrix& M)
00191 {
00192    int nr = M.nrows();
00193    int nc = M.ncols();
00194    for (int i = 1; i <= nr; ++i) for (int j = 1; j <= nc; ++ j)
00195    M(i, j) = MWC.Next();
00196 }
00197    
00198 
00199 #ifdef use_namespace
00200 }
00201 using namespace NEWMAT;
00202 #endif
00203 
00204 
00205 //*************************** main program **********************************
00206 
00207 void TestTypeAdd();                            // test +
00208 void TestTypeMult();                           // test *
00209 void TestTypeConcat();                         // test |
00210 void TestTypeSP();                             // test SP
00211 void TestTypeKP();                             // test KP
00212 void TestTypeOrder();                          // test >=
00213 
00214 
00215 int main()
00216 {
00217    time_lapse tl;      // measure program run time
00218    Real* s1; Real* s2; Real* s3; Real* s4;
00219    cout << "\nBegin test\n";   // Forces cout to allocate memory at beginning
00220    cout << "Now print a real number: " << 3.14159265 << endl;
00221    // Throw exception to set up exception buffer
00222 #ifndef DisableExceptions
00223    Try { Throw(BaseException("Just a dummy\n")); }
00224    CatchAll {}
00225 #else
00226    cout << "Not doing exceptions\n";
00227 #endif
00228    { Matrix A1(40,200); s1 = A1.Store(); }
00229    { Matrix A1(1,1); s3 = A1.Store(); }
00230    {
00231       Tracer et("Matrix test program");
00232 
00233       Matrix A(25,150);
00234       {
00235          int i;
00236          RowVector A(8);
00237          for (i=1;i<=7;i++) A(i)=0.0; A(8)=1.0;
00238          Print(A);
00239       }
00240       cout << "\n";
00241 
00242       TestTypeAdd(); TestTypeMult(); TestTypeConcat();
00243       TestTypeSP(); TestTypeKP(); TestTypeOrder();
00244 
00245       Try { 
00246          trymat1();
00247          trymat2();
00248          trymat3();
00249          trymat4();
00250          trymat5();
00251          trymat6();
00252          trymat7();
00253          trymat8();
00254          trymat9();
00255          trymata();
00256          trymatb();
00257          trymatc();
00258          trymatd();
00259          trymate();
00260          trymatf();
00261          trymatg();
00262          trymath();
00263          trymati();
00264          trymatj();
00265          trymatk();
00266          trymatl();
00267          trymatm();
00268 
00269          cout << "\nEnd of tests\n";
00270       }
00271       CatchAll
00272       {
00273          cout << "\nTest program fails - exception generated\n\n";
00274          cout << BaseException::what();
00275       }
00276 
00277 
00278    }
00279 
00280    { Matrix A1(40,200); s2 = A1.Store(); }
00281    cout << "\n(The following memory checks are probably not valid with all\n";
00282    cout << "compilers - see documentation)\n";
00283    cout << "\nChecking for lost memory (large block): "
00284       << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
00285    if (s1 != s2) cout << " - see section 2.8\n\n"; else cout << " - ok\n\n";
00286    { Matrix A1(1,1); s4 = A1.Store(); }
00287    cout << "\nChecking for lost memory (small block): "
00288       << (unsigned long)s3 << " " << (unsigned long)s4 << " ";
00289    if (s3 != s4) cout << " - see section 2.8\n\n"; else cout << " - ok\n\n";
00290 
00291    // check for Pentium bug
00292    PentiumCheck(4195835L,3145727L);
00293    PentiumCheck(5244795L,3932159L);
00294 
00295 #ifdef DO_FREE_CHECK
00296    FreeCheck::Status();
00297 #endif
00298    return 0;
00299 }
00300 
00301 
00302 
00303 
00304 //************************ test type manipulation **************************/
00305 
00306 
00307 // These functions may cause problems for Glockenspiel 2.0c; they are used
00308 // only for testing so you can delete them
00309 
00310 
00311 void TestTypeAdd()
00312 {
00313    MatrixType list[13];
00314    list[0] = MatrixType::UT;
00315    list[1] = MatrixType::LT;
00316    list[2] = MatrixType::Rt;
00317    list[3] = MatrixType::Sq;
00318    list[4] = MatrixType::Sm;
00319    list[5] = MatrixType::Sk;
00320    list[6] = MatrixType::Dg;
00321    list[7] = MatrixType::Id;
00322    list[8] = MatrixType::BM;
00323    list[9] = MatrixType::UB;
00324    list[10] = MatrixType::LB;
00325    list[11] = MatrixType::SB;
00326    list[12] = MatrixType::KB;
00327 
00328    cout << "+     ";
00329    int i;
00330    for (i=0; i<MatrixType::nTypes(); i++) cout << list[i].Value() << " ";
00331    cout << "\n";
00332    for (i=0; i<MatrixType::nTypes(); i++)
00333         {
00334                 cout << list[i].Value() << " ";
00335       for (int j=0; j<MatrixType::nTypes(); j++)
00336          cout << (list[j]+list[i]).Value() << " ";
00337       cout << "\n";
00338    }
00339    cout << "\n";
00340 }
00341 
00342 void TestTypeMult()
00343 {
00344    MatrixType list[13];
00345    list[0] = MatrixType::UT;
00346    list[1] = MatrixType::LT;
00347    list[2] = MatrixType::Rt;
00348    list[3] = MatrixType::Sq;
00349    list[4] = MatrixType::Sm;
00350    list[5] = MatrixType::Sk;
00351    list[6] = MatrixType::Dg;
00352    list[7] = MatrixType::Id;
00353    list[8] = MatrixType::BM;
00354    list[9] = MatrixType::UB;
00355    list[10] = MatrixType::LB;
00356    list[11] = MatrixType::SB;
00357    list[12] = MatrixType::KB;
00358 
00359    cout << "*     ";
00360    int i;
00361    for (i=0; i<MatrixType::nTypes(); i++)
00362       cout << list[i].Value() << " ";
00363    cout << "\n";
00364    for (i=0; i<MatrixType::nTypes(); i++)
00365    {
00366                 cout << list[i].Value() << " ";
00367       for (int j=0; j<MatrixType::nTypes(); j++)
00368          cout << (list[j]*list[i]).Value() << " ";
00369       cout << "\n";
00370    }
00371    cout << "\n";
00372 }
00373 
00374 void TestTypeConcat()
00375 {
00376    MatrixType list[13];
00377    list[0] = MatrixType::UT;
00378    list[1] = MatrixType::LT;
00379    list[2] = MatrixType::Rt;
00380    list[3] = MatrixType::Sq;
00381    list[4] = MatrixType::Sm;
00382    list[5] = MatrixType::Sk;
00383    list[6] = MatrixType::Dg;
00384    list[7] = MatrixType::Id;
00385    list[8] = MatrixType::BM;
00386    list[9] = MatrixType::UB;
00387    list[10] = MatrixType::LB;
00388    list[11] = MatrixType::SB;
00389    list[12] = MatrixType::KB;
00390 
00391    cout << "|     ";
00392    int i;
00393    for (i=0; i<MatrixType::nTypes(); i++)
00394       cout << list[i].Value() << " ";
00395    cout << "\n";
00396    for (i=0; i<MatrixType::nTypes(); i++)
00397    {
00398                 cout << list[i].Value() << " ";
00399       for (int j=0; j<MatrixType::nTypes(); j++)
00400          cout << (list[j] | list[i]).Value() << " ";
00401       cout << "\n";
00402    }
00403    cout << "\n";
00404 }
00405 
00406 void TestTypeSP()
00407 {
00408    MatrixType list[13];
00409    list[0] = MatrixType::UT;
00410    list[1] = MatrixType::LT;
00411    list[2] = MatrixType::Rt;
00412    list[3] = MatrixType::Sq;
00413    list[4] = MatrixType::Sm;
00414    list[5] = MatrixType::Sk;
00415    list[6] = MatrixType::Dg;
00416    list[7] = MatrixType::Id;
00417    list[8] = MatrixType::BM;
00418    list[9] = MatrixType::UB;
00419    list[10] = MatrixType::LB;
00420    list[11] = MatrixType::SB;
00421    list[12] = MatrixType::KB;
00422 
00423    cout << "SP    ";
00424    int i;
00425    for (i=0; i<MatrixType::nTypes(); i++)
00426                 cout << list[i].Value() << " ";
00427    cout << "\n";
00428    for (i=0; i<MatrixType::nTypes(); i++)
00429    {
00430                 cout << list[i].Value() << " ";
00431       for (int j=0; j<MatrixType::nTypes(); j++)
00432          cout << (list[j].SP(list[i])).Value() << " ";
00433       cout << "\n";
00434    }
00435    cout << "\n";
00436 }
00437 
00438 void TestTypeKP()
00439 {
00440    MatrixType list[13];
00441    list[0] = MatrixType::UT;
00442    list[1] = MatrixType::LT;
00443    list[2] = MatrixType::Rt;
00444    list[3] = MatrixType::Sq;
00445    list[4] = MatrixType::Sm;
00446    list[5] = MatrixType::Sk;
00447    list[6] = MatrixType::Dg;
00448    list[7] = MatrixType::Id;
00449    list[8] = MatrixType::BM;
00450    list[9] = MatrixType::UB;
00451    list[10] = MatrixType::LB;
00452    list[11] = MatrixType::SB;
00453    list[12] = MatrixType::KB;
00454 
00455    cout << "KP    ";
00456    int i;
00457    for (i=0; i<MatrixType::nTypes(); i++)
00458                 cout << list[i].Value() << " ";
00459    cout << "\n";
00460    for (i=0; i<MatrixType::nTypes(); i++)
00461    {
00462                 cout << list[i].Value() << " ";
00463       for (int j=0; j<MatrixType::nTypes(); j++)
00464          cout << (list[j].KP(list[i])).Value() << " ";
00465       cout << "\n";
00466    }
00467    cout << "\n";
00468 }
00469 
00470 void TestTypeOrder()
00471 {
00472    MatrixType list[13];
00473    list[0] = MatrixType::UT;
00474    list[1] = MatrixType::LT;
00475    list[2] = MatrixType::Rt;
00476    list[3] = MatrixType::Sq;
00477    list[4] = MatrixType::Sm;
00478    list[5] = MatrixType::Sk;
00479    list[6] = MatrixType::Dg;
00480    list[7] = MatrixType::Id;
00481    list[8] = MatrixType::BM;
00482    list[9] = MatrixType::UB;
00483    list[10] = MatrixType::LB;
00484    list[11] = MatrixType::SB;
00485    list[12] = MatrixType::KB;
00486 
00487    cout << ">=    ";
00488    int i;
00489    for (i = 0; i<MatrixType::nTypes(); i++)
00490       cout << list[i].Value() << " ";
00491    cout << "\n";
00492    for (i=0; i<MatrixType::nTypes(); i++)
00493    {
00494       cout << list[i].Value() << " ";
00495       for (int j=0; j<MatrixType::nTypes(); j++)
00496          cout << ((list[j]>=list[i]) ? "Yes   " : "No    ");
00497       cout << "\n";
00498    }
00499    cout << "\n";
00500 }
00501 
00502 
00503 //************** elapsed time class ****************
00504 
00505 time_lapse::time_lapse()
00506 {
00507    start_time = ((double)clock())/(double)CLOCKS_PER_SEC;
00508 }
00509 
00510 time_lapse::~time_lapse()
00511 {
00512    double time = ((double)clock())/(double)CLOCKS_PER_SEC - start_time;
00513    cout << "Elapsed (processor) time = " << setprecision(2) << time << " seconds" << endl;
00514    cout << endl;
00515 }
00516 
00517 
00518 
00519 
00521 
00522 
00523 


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