DenseLinAlgPack_TestGenMatrixClass.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <iomanip>
00030 #include <ostream>
00031 #include <vector>
00032 #include <typeinfo>
00033 
00034 #include "DenseLinAlgPack_TestDenseLinAlgPack.hpp"
00035 #include "DenseLinAlgPack_DMatrixClass.hpp"
00036 #include "DenseLinAlgPack_DVectorOut.hpp"
00037 #include "DenseLinAlgPack_DMatrixOut.hpp"
00038 #include "DenseLinAlgPack_MatVecCompare.hpp"
00039 
00040 namespace {
00041 
00042 using DenseLinAlgPack::sqrt_eps;
00043 
00044 // Check consistency of row(), col(), diag() and operator()().
00045 template<class M_t>
00046 void check_access( M_t& M, typename M_t::size_type row_offset, typename M_t::size_type col_offset
00047   , std::ostream* out, bool* success )
00048 {
00049   if(out)
00050     *out  << "Checking M(i,j) == M.row(i)(j) == M.col(j)(i) == "
00051         << "M.diag(...)(...) == "
00052         << "(i + "<<row_offset<<") + 0.1*(j+"<<col_offset<<") : ";
00053 
00054   bool result = true;
00055 
00056   for( typename M_t::size_type i = 1; i <= M.rows(); ++i ) {
00057     for( typename M_t::size_type j = 1; j <= M.rows(); ++j ) {
00058       const typename M_t::value_type
00059         Mij = M(i,j);
00060       typename M_t::value_type
00061         val = (i+row_offset)+0.1*(j+col_offset);
00062       if( ::fabs(Mij-val) > sqrt_eps ) {
00063         result = false;
00064         if(out) *out << "(M("<<i<<","<<j<<") -> "<<Mij<<") != "<<val<<std::endl;
00065       }
00066       if( Mij != (val = M.row(i)(j)) ) {
00067         result = false;
00068         if(out) *out << "M("<<i<<","<<j<<") != (M.row("<<i<<")("<<j<<") -> "<<val<<")\n";
00069       }
00070       if( Mij != (val = M.col(j)(i)) ) {
00071         result = false;
00072         if(out) *out << "M("<<i<<","<<j<<") != (M.col("<<j<<")("<<i<<") -> "<<val<<")\n";
00073       }
00074       const int k = ( i > j ? -i + j : j - i );
00075       const typename M_t::size_type k_i = ( i > j ? j : i );
00076       if( Mij != (val = M.diag(k)(k_i) ) ) {
00077         result = false;
00078         if(out) *out << "M("<<i<<","<<j<<") != (M.diag("<<k<<")("<<k_i<<") -> "<<val<<")\n";
00079       }
00080     }
00081   }
00082   if(out) *out << result << std::endl;
00083   if(!result) *success = false;
00084 }
00085 
00086 // Print out a string for overlap
00087 const char* overlap_str( DenseLinAlgPack::EOverLap overlap ) {
00088   switch(overlap) {
00089     case DenseLinAlgPack::NO_OVERLAP:
00090       return "NO_OVERLAP";
00091     case DenseLinAlgPack::SOME_OVERLAP:
00092       return "SOME_OVERLAP";
00093     case DenseLinAlgPack::SAME_MEM:
00094       return "SAME_MEM";
00095   }
00096   return "Invalid value for EOverLap";
00097 }
00098 
00099 } // end namespace
00100 
00101 bool DenseLinAlgPack::TestingPack::TestGenMatrixClass(std::ostream* out)
00102 {
00103 
00104   using DenseLinAlgPack::comp;
00105   using DenseLinAlgPack::sqrt_eps;
00106 
00107   bool success = true;
00108   bool result;
00109 
00110   if(out)
00111     *out  << "\n****************************************************"
00112         << "\n*** Testing DMatrix and DMatrixSlice classes ***"
00113         << "\n****************************************************\n"
00114         << std::boolalpha;
00115 
00116   try {
00117 
00118   const size_type
00119     m = 6,
00120     n = 8;
00121 
00122   const value_type
00123     ptr[m*n] =
00124        {  1.1,  2.1,  3.1,  4.1,  5.1,  6.1,
00125         1.2,  2.2,  3.2,  4.2,  5.2,  6.2,
00126         1.3,  2.3,  3.3,  4.3,  5.3,  6.3,
00127         1.4,  2.4,  3.4,  4.4,  5.4,  6.4,
00128         1.5,  2.5,  3.5,  4.5,  5.5,  6.5,
00129         1.6,  2.6,  3.6,  4.6,  5.6,  6.6,
00130         1.7,  2.7,  3.7,  4.7,  5.7,  6.7,
00131         1.8,  2.8,  3.8,  4.8,  5.8,  6.8 };
00132 
00133   // /////////////////////////////
00134   // Test Constructors
00135 
00136   if(out)
00137     *out  << "\n***\n*** Testing constructors\n***\n";
00138 
00139   // DMatrixSlice
00140   if(out) *out << "\nGenMatrixSlice gms1;\n";
00141   DMatrixSlice gms1;
00142   if(out) *out << "gms1 =\n" << gms1;
00143   update_success( result = (gms1.rows() == 0 && gms1.cols() == 0 ), &success );
00144   if(out)
00145     *out  << "((gms1.rows() -> "<<gms1.rows()
00146         << ") == 0 && (gms1.cols() -> "<<gms1.cols()<<") == 0 ) : "
00147         << result << std::endl;
00148 
00149 
00150   if(out) *out << "\nGenMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );\n";
00151   const DMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );
00152   if(out) *out << "gms2 =\n" << gms2;
00153 
00154   if(out) *out << "\nGenMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );\n";
00155   const DMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );
00156   if(out) *out << "gms3 =\n" << gms3;
00157 
00158   // DMatrix
00159 
00160   if(out) *out << "\nGenMatrix gm1;\n";
00161   DMatrix gm1;  
00162   if(out) *out << "gm1 =\n" << gm1();
00163   update_success( result = (gm1.rows() == 0 && gm1.cols() == 0 ), &success );
00164   if(out)
00165     *out  << "((gm1.rows() -> "<<gm1.rows()
00166         << ") == 0 && (gm1.cols() -> "<<gm1.cols()<<") == 0 ) : "
00167         << result << std::endl;
00168   
00169   if(out) *out << "\nGenMatrix gm2(m,n);\n";
00170   DMatrix gm2(m,n);
00171   if(out) *out << "gm2 =\n" << gm2();
00172 
00173   if(out) *out << "\nGenMatrix gm3(1.0,m,n);\n";
00174   DMatrix gm3(1.0,m,n);
00175   if(out) *out << "gm3 =\n" << gm3();
00176   update_success( result = comp( gm3(), 1.0 ), &success );
00177   if(out) *out << "gm3 == 1.0 : " << result << std::endl;
00178 
00179   if(out) *out << "\nGenMatrix gm4(ptr,m,n);\n";
00180   DMatrix gm4(ptr,m,n);
00181   if(out) *out << "gm4 =\n" << gm4();
00182 
00183   if(out) *out << "\nGenMatrix gm5(gms2);\n";
00184   DMatrix gm5(gms2);
00185   if(out) *out << "gm5 =\n" << gm5();
00186 
00187   // ////////////////////////////
00188   // Test DMatrixSlice binding
00189 
00190   if(out)
00191     *out  << "\n***\n*** Testing DMatrixSlice binding\n***\n";
00192 
00193   if(out) *out << "\ngms1.bind(gm4());\n";
00194   gms1.bind(gm4());
00195   if(out) *out << "gms1 =\n" << gms1();
00196 
00197   // ////////////////////////////
00198   // Test DMatrix resizing
00199 
00200   if(out)
00201     *out  << "\n***\n*** Testing DMatrix resizing\n***\n";
00202 
00203   if(out) *out << "\ngm1.resize(m,n,1.0);\n";
00204   gm1.resize(m,n,1.0);
00205   if(out) *out << "gm1 =\n" << gm1();
00206   update_success( result = comp( gm1(), 1.0 ), &success );
00207   if(out) *out << "gm1 == 1.0 : " << result << std::endl;
00208 
00209   // ///////////////////////////////////////////////
00210   // Test row, col, diag access and element access
00211 
00212   // DMatrixSlice
00213 
00214   if(out)
00215     *out  << "\n***\n*** Testing row, col, diag access and element access\n***\n";
00216 
00217   if(out) *out << "\nLet M = gms1\n";
00218   check_access( gms1, 0, 0, out, &success );
00219 
00220   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gms1)\n";
00221   check_access( const_cast<const DMatrixSlice&>(gms1), 0, 0, out, &success );
00222 
00223   // DMatrix
00224 
00225   if(out) *out << "\nLet M = gm4\n";
00226   check_access( gm4, 0, 0, out, &success );
00227 
00228   if(out) *out << "\nLet M = const_cast<const DMatrix&>(gm4)\n";
00229   check_access( const_cast<const DMatrix&>(gm4), 0, 0, out, &success );
00230 
00231   // ////////////////////////////
00232   // Test submatrix access
00233 
00234   if(out)
00235     *out  << "\n***\n*** Testing submatrix access\n***\n";
00236 
00237   if(out) *out << "\nRange1D r_rng(2,m-1), c_rng(2,n-1);\n";
00238   Range1D r_rng(2,m-1), c_rng(2,n-1);
00239 
00240   // DMatrixSlice
00241 
00242   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng)\n";
00243   gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
00244   if(out) *out << "M =\n" << gms1;
00245   check_access( gms1, 1, 1, out, &success );
00246 
00247   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(r_rng,c_rng))\n";
00248   gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
00249   if(out) *out << "M =\n" << gms1;
00250   check_access( gms1, 1, 1, out, &success );
00251 
00252   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1)\n";
00253   gms1.bind(const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
00254   if(out) *out << "M =\n" << gms1;
00255   check_access( gms1, 1, 1, out, &success );
00256 
00257   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(2,m-1,2,n-1))\n";
00258   gms1.bind( const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
00259   if(out) *out << "M =\n" << gms1;
00260   check_access( gms1, 1, 1, out, &success );
00261 
00262   // DMatrix
00263 
00264   if(out) *out << "\nLet M = gm4(r_rng,c_rng)\n";
00265   gms1.bind( gm4(r_rng,c_rng) );
00266   if(out) *out << "M =\n" << gms1;
00267   check_access( gms1, 1, 1, out, &success );
00268 
00269   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(r_rng,c_rng)\n";
00270   gms1.bind( const_cast<const DMatrix&>(gm4)(r_rng,c_rng) );
00271   if(out) *out << "M =\n" << gms1;
00272   check_access( gms1, 1, 1, out, &success );
00273 
00274   if(out) *out << "\nLet M = gm4(2,m-1,2,n-1)\n";
00275   gms1.bind( gm4(2,m-1,2,n-1) );
00276   if(out) *out << "M =\n" << gms1;
00277   check_access( gms1, 1, 1, out, &success );
00278 
00279   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(2,m-1,2,n-1)\n";
00280   gms1.bind( const_cast<const DMatrix&>(gm4)(2,m-1,2,n-1) );
00281   if(out) *out << "M =\n" << gms1;
00282   check_access( gms1, 1, 1, out, &success );
00283 
00284   // ////////////////////
00285   // Test matrix overlap
00286 
00287   if(out)
00288     *out  << "\n***\n*** matrix overlap\n***\n";
00289 
00290   EOverLap ovlap;
00291 
00292   // DMatrixSlice
00293 
00294   if(out) *out << "(gms2.overlap(gms2) -> ";
00295   ovlap = gms2.overlap(gms2);
00296   result = update_success( ovlap == SAME_MEM, &success );
00297   if(out) *out  << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
00298 
00299   if(out) *out << "(gms2.overlap(gms2(r_rng,c_rng)) -> ";
00300   ovlap = gms2.overlap(gms2(r_rng,c_rng));
00301   result = update_success( ovlap == SOME_OVERLAP, &success );
00302   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00303 
00304   if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n)) -> ";
00305   ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n));
00306   result = update_success( ovlap == SOME_OVERLAP, &success );
00307   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00308 
00309   if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n)) -> ";
00310   ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n));
00311   result = update_success( ovlap == NO_OVERLAP, &success );
00312   if(out) *out  << overlap_str(ovlap) << ") == NO_OVERLAP : " << result << std::endl;
00313 
00314   // DMatrix
00315 
00316   if(out) *out << "(gm4.overlap(gm4) -> ";
00317   ovlap = gm4.overlap(gm4());
00318   result = update_success( ovlap == SAME_MEM, &success );
00319   if(out) *out  << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
00320 
00321   if(out) *out << "(gm4.overlap(gm4(r_rng,c_rng)) -> ";
00322   ovlap = gm4.overlap(gm4(r_rng,c_rng));
00323   result = update_success( ovlap == SOME_OVERLAP, &success );
00324   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00325 
00326   // //////////////////////////////////////////////////////
00327   // Test vector overlap (continuation of vector testing)
00328   //
00329   // ToDo: Finish this someday once you get it figured out.
00330 
00331   // ///////////////////////////
00332   // Test assignment operators
00333 
00334   if(out)
00335     *out  << "\n***\n*** assignment operators\n***\n";
00336 
00337   // DMatrixSlice
00338 
00339   if(out) *out << "\ngms1.bind(gm1());\n";
00340   gms1.bind(gm1());
00341 
00342   if(out) *out << "\ngms1 = 2.0;\n";
00343   gms1 = 2.0;
00344   if(out) *out << "gms1 =\n" << gms1;
00345   update_success( result = comp(gms1,2.0), &success );
00346   if(out) *out << "gms1 == 2.0 : " << result << std::endl; 
00347 
00348   if(out) *out << "\ngms1 = gms2;\n";
00349   gms1 = gms2;
00350   if(out) *out << "gms1 =\n" << gms1;
00351   update_success( result = comp(gms1,gms2), &success );
00352   if(out) *out << "gms1 == gms2 : " << result << std::endl; 
00353 
00354   // DMatrix
00355 
00356   if(out) *out << "\ngm1 = 3.0;\n";
00357   gm1 = 3.0;
00358   if(out) *out << "gm1 =\n" << gm1();
00359   update_success( result = comp(gm1(),3.0), &success );
00360   if(out) *out << "gm1 == 3.0 : " << result << std::endl; 
00361 
00362   if(out) *out << "\ngm1.resize(0,0); gm1 = gms2;\n";
00363   gm1.resize(0,0);
00364   gm1 = gms2;
00365   if(out) *out << "gm1 =\n" << gm1();
00366   update_success( result = comp(gm1(),gms2()), &success );
00367   if(out) *out << "gm1 == gms2 : " << result << std::endl; 
00368 
00369   if(out) *out << "\ngm1.resize(0,0); gm1 = gm4;\n";
00370   gm1.resize(0,0);
00371   gm1 = gm4;
00372   if(out) *out << "gm1 =\n" << gm1();
00373   update_success( result = comp(gm1(),gm4()), &success );
00374   if(out) *out << "gm1 == gm4 : " << result << std::endl; 
00375 
00376   } // end try
00377   catch( const std::exception& excpt ) {
00378     success = false;
00379     if(out)
00380       (*out)  << "\nError, a standard exception was thrown: "
00381           << typeid(excpt).name() << ": "
00382           << excpt.what() << std::endl; 
00383   }
00384   catch(...) {
00385     success = false;
00386     if(out)
00387       (*out)  << "\nError, an unknown exception was thrown\n";
00388   }
00389 
00390   if(out) {
00391     if(success)
00392       (*out)
00393         << "\n*** Congradulations, DMatrix and DMatrixSlice seem to check out. ***\n";
00394     else
00395       (*out)
00396         << "\n*** Oops, all of the tests for DMatrix and DMatrixSlice "
00397           "where not successful. ***\n";
00398   }
00399 
00400   return success;
00401 }
00402 

Generated on Thu Sep 18 12:35:16 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1