MOOCHO (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <iomanip>
00043 #include <ostream>
00044 #include <vector>
00045 #include <typeinfo>
00046 
00047 #include "DenseLinAlgPack_TestDenseLinAlgPack.hpp"
00048 #include "DenseLinAlgPack_DMatrixClass.hpp"
00049 #include "DenseLinAlgPack_DVectorOut.hpp"
00050 #include "DenseLinAlgPack_DMatrixOut.hpp"
00051 #include "DenseLinAlgPack_MatVecCompare.hpp"
00052 
00053 namespace {
00054 
00055 using DenseLinAlgPack::sqrt_eps;
00056 
00057 // Check consistency of row(), col(), diag() and operator()().
00058 template<class M_t>
00059 void check_access( M_t& M, typename M_t::size_type row_offset, typename M_t::size_type col_offset
00060   , std::ostream* out, bool* success )
00061 {
00062   if(out)
00063     *out  << "Checking M(i,j) == M.row(i)(j) == M.col(j)(i) == "
00064         << "M.diag(...)(...) == "
00065         << "(i + "<<row_offset<<") + 0.1*(j+"<<col_offset<<") : ";
00066 
00067   bool result = true;
00068 
00069   for( typename M_t::size_type i = 1; i <= M.rows(); ++i ) {
00070     for( typename M_t::size_type j = 1; j <= M.rows(); ++j ) {
00071       const typename M_t::value_type
00072         Mij = M(i,j);
00073       typename M_t::value_type
00074         val = (i+row_offset)+0.1*(j+col_offset);
00075       if( ::fabs(Mij-val) > sqrt_eps ) {
00076         result = false;
00077         if(out) *out << "(M("<<i<<","<<j<<") -> "<<Mij<<") != "<<val<<std::endl;
00078       }
00079       if( Mij != (val = M.row(i)(j)) ) {
00080         result = false;
00081         if(out) *out << "M("<<i<<","<<j<<") != (M.row("<<i<<")("<<j<<") -> "<<val<<")\n";
00082       }
00083       if( Mij != (val = M.col(j)(i)) ) {
00084         result = false;
00085         if(out) *out << "M("<<i<<","<<j<<") != (M.col("<<j<<")("<<i<<") -> "<<val<<")\n";
00086       }
00087       const int k = ( i > j ? -i + j : j - i );
00088       const typename M_t::size_type k_i = ( i > j ? j : i );
00089       if( Mij != (val = M.diag(k)(k_i) ) ) {
00090         result = false;
00091         if(out) *out << "M("<<i<<","<<j<<") != (M.diag("<<k<<")("<<k_i<<") -> "<<val<<")\n";
00092       }
00093     }
00094   }
00095   if(out) *out << result << std::endl;
00096   if(!result) *success = false;
00097 }
00098 
00099 // Print out a string for overlap
00100 const char* overlap_str( DenseLinAlgPack::EOverLap overlap ) {
00101   switch(overlap) {
00102     case DenseLinAlgPack::NO_OVERLAP:
00103       return "NO_OVERLAP";
00104     case DenseLinAlgPack::SOME_OVERLAP:
00105       return "SOME_OVERLAP";
00106     case DenseLinAlgPack::SAME_MEM:
00107       return "SAME_MEM";
00108   }
00109   return "Invalid value for EOverLap";
00110 }
00111 
00112 } // end namespace
00113 
00114 bool DenseLinAlgPack::TestingPack::TestGenMatrixClass(std::ostream* out)
00115 {
00116 
00117   using DenseLinAlgPack::comp;
00118   using DenseLinAlgPack::sqrt_eps;
00119 
00120   bool success = true;
00121   bool result;
00122 
00123   if(out)
00124     *out  << "\n****************************************************"
00125         << "\n*** Testing DMatrix and DMatrixSlice classes ***"
00126         << "\n****************************************************\n"
00127         << std::boolalpha;
00128 
00129   try {
00130 
00131   const size_type
00132     m = 6,
00133     n = 8;
00134 
00135   const value_type
00136     ptr[m*n] =
00137        {  1.1,  2.1,  3.1,  4.1,  5.1,  6.1,
00138         1.2,  2.2,  3.2,  4.2,  5.2,  6.2,
00139         1.3,  2.3,  3.3,  4.3,  5.3,  6.3,
00140         1.4,  2.4,  3.4,  4.4,  5.4,  6.4,
00141         1.5,  2.5,  3.5,  4.5,  5.5,  6.5,
00142         1.6,  2.6,  3.6,  4.6,  5.6,  6.6,
00143         1.7,  2.7,  3.7,  4.7,  5.7,  6.7,
00144         1.8,  2.8,  3.8,  4.8,  5.8,  6.8 };
00145 
00146   // /////////////////////////////
00147   // Test Constructors
00148 
00149   if(out)
00150     *out  << "\n***\n*** Testing constructors\n***\n";
00151 
00152   // DMatrixSlice
00153   if(out) *out << "\nGenMatrixSlice gms1;\n";
00154   DMatrixSlice gms1;
00155   if(out) *out << "gms1 =\n" << gms1;
00156   update_success( result = (gms1.rows() == 0 && gms1.cols() == 0 ), &success );
00157   if(out)
00158     *out  << "((gms1.rows() -> "<<gms1.rows()
00159         << ") == 0 && (gms1.cols() -> "<<gms1.cols()<<") == 0 ) : "
00160         << result << std::endl;
00161 
00162 
00163   if(out) *out << "\nGenMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );\n";
00164   const DMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );
00165   if(out) *out << "gms2 =\n" << gms2;
00166 
00167   if(out) *out << "\nGenMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );\n";
00168   const DMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );
00169   if(out) *out << "gms3 =\n" << gms3;
00170 
00171   // DMatrix
00172 
00173   if(out) *out << "\nGenMatrix gm1;\n";
00174   DMatrix gm1;  
00175   if(out) *out << "gm1 =\n" << gm1();
00176   update_success( result = (gm1.rows() == 0 && gm1.cols() == 0 ), &success );
00177   if(out)
00178     *out  << "((gm1.rows() -> "<<gm1.rows()
00179         << ") == 0 && (gm1.cols() -> "<<gm1.cols()<<") == 0 ) : "
00180         << result << std::endl;
00181   
00182   if(out) *out << "\nGenMatrix gm2(m,n);\n";
00183   DMatrix gm2(m,n);
00184   if(out) *out << "gm2 =\n" << gm2();
00185 
00186   if(out) *out << "\nGenMatrix gm3(1.0,m,n);\n";
00187   DMatrix gm3(1.0,m,n);
00188   if(out) *out << "gm3 =\n" << gm3();
00189   update_success( result = comp( gm3(), 1.0 ), &success );
00190   if(out) *out << "gm3 == 1.0 : " << result << std::endl;
00191 
00192   if(out) *out << "\nGenMatrix gm4(ptr,m,n);\n";
00193   DMatrix gm4(ptr,m,n);
00194   if(out) *out << "gm4 =\n" << gm4();
00195 
00196   if(out) *out << "\nGenMatrix gm5(gms2);\n";
00197   DMatrix gm5(gms2);
00198   if(out) *out << "gm5 =\n" << gm5();
00199 
00200   // ////////////////////////////
00201   // Test DMatrixSlice binding
00202 
00203   if(out)
00204     *out  << "\n***\n*** Testing DMatrixSlice binding\n***\n";
00205 
00206   if(out) *out << "\ngms1.bind(gm4());\n";
00207   gms1.bind(gm4());
00208   if(out) *out << "gms1 =\n" << gms1();
00209 
00210   // ////////////////////////////
00211   // Test DMatrix resizing
00212 
00213   if(out)
00214     *out  << "\n***\n*** Testing DMatrix resizing\n***\n";
00215 
00216   if(out) *out << "\ngm1.resize(m,n,1.0);\n";
00217   gm1.resize(m,n,1.0);
00218   if(out) *out << "gm1 =\n" << gm1();
00219   update_success( result = comp( gm1(), 1.0 ), &success );
00220   if(out) *out << "gm1 == 1.0 : " << result << std::endl;
00221 
00222   // ///////////////////////////////////////////////
00223   // Test row, col, diag access and element access
00224 
00225   // DMatrixSlice
00226 
00227   if(out)
00228     *out  << "\n***\n*** Testing row, col, diag access and element access\n***\n";
00229 
00230   if(out) *out << "\nLet M = gms1\n";
00231   check_access( gms1, 0, 0, out, &success );
00232 
00233   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gms1)\n";
00234   check_access( const_cast<const DMatrixSlice&>(gms1), 0, 0, out, &success );
00235 
00236   // DMatrix
00237 
00238   if(out) *out << "\nLet M = gm4\n";
00239   check_access( gm4, 0, 0, out, &success );
00240 
00241   if(out) *out << "\nLet M = const_cast<const DMatrix&>(gm4)\n";
00242   check_access( const_cast<const DMatrix&>(gm4), 0, 0, out, &success );
00243 
00244   // ////////////////////////////
00245   // Test submatrix access
00246 
00247   if(out)
00248     *out  << "\n***\n*** Testing submatrix access\n***\n";
00249 
00250   if(out) *out << "\nRange1D r_rng(2,m-1), c_rng(2,n-1);\n";
00251   Range1D r_rng(2,m-1), c_rng(2,n-1);
00252 
00253   // DMatrixSlice
00254 
00255   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng)\n";
00256   gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
00257   if(out) *out << "M =\n" << gms1;
00258   check_access( gms1, 1, 1, out, &success );
00259 
00260   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(r_rng,c_rng))\n";
00261   gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
00262   if(out) *out << "M =\n" << gms1;
00263   check_access( gms1, 1, 1, out, &success );
00264 
00265   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1)\n";
00266   gms1.bind(const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
00267   if(out) *out << "M =\n" << gms1;
00268   check_access( gms1, 1, 1, out, &success );
00269 
00270   if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(2,m-1,2,n-1))\n";
00271   gms1.bind( const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
00272   if(out) *out << "M =\n" << gms1;
00273   check_access( gms1, 1, 1, out, &success );
00274 
00275   // DMatrix
00276 
00277   if(out) *out << "\nLet M = gm4(r_rng,c_rng)\n";
00278   gms1.bind( gm4(r_rng,c_rng) );
00279   if(out) *out << "M =\n" << gms1;
00280   check_access( gms1, 1, 1, out, &success );
00281 
00282   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(r_rng,c_rng)\n";
00283   gms1.bind( const_cast<const DMatrix&>(gm4)(r_rng,c_rng) );
00284   if(out) *out << "M =\n" << gms1;
00285   check_access( gms1, 1, 1, out, &success );
00286 
00287   if(out) *out << "\nLet M = gm4(2,m-1,2,n-1)\n";
00288   gms1.bind( gm4(2,m-1,2,n-1) );
00289   if(out) *out << "M =\n" << gms1;
00290   check_access( gms1, 1, 1, out, &success );
00291 
00292   if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(2,m-1,2,n-1)\n";
00293   gms1.bind( const_cast<const DMatrix&>(gm4)(2,m-1,2,n-1) );
00294   if(out) *out << "M =\n" << gms1;
00295   check_access( gms1, 1, 1, out, &success );
00296 
00297   // ////////////////////
00298   // Test matrix overlap
00299 
00300   if(out)
00301     *out  << "\n***\n*** matrix overlap\n***\n";
00302 
00303   EOverLap ovlap;
00304 
00305   // DMatrixSlice
00306 
00307   if(out) *out << "(gms2.overlap(gms2) -> ";
00308   ovlap = gms2.overlap(gms2);
00309   result = update_success( ovlap == SAME_MEM, &success );
00310   if(out) *out  << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
00311 
00312   if(out) *out << "(gms2.overlap(gms2(r_rng,c_rng)) -> ";
00313   ovlap = gms2.overlap(gms2(r_rng,c_rng));
00314   result = update_success( ovlap == SOME_OVERLAP, &success );
00315   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00316 
00317   if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n)) -> ";
00318   ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n));
00319   result = update_success( ovlap == SOME_OVERLAP, &success );
00320   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00321 
00322   if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n)) -> ";
00323   ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n));
00324   result = update_success( ovlap == NO_OVERLAP, &success );
00325   if(out) *out  << overlap_str(ovlap) << ") == NO_OVERLAP : " << result << std::endl;
00326 
00327   // DMatrix
00328 
00329   if(out) *out << "(gm4.overlap(gm4) -> ";
00330   ovlap = gm4.overlap(gm4());
00331   result = update_success( ovlap == SAME_MEM, &success );
00332   if(out) *out  << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
00333 
00334   if(out) *out << "(gm4.overlap(gm4(r_rng,c_rng)) -> ";
00335   ovlap = gm4.overlap(gm4(r_rng,c_rng));
00336   result = update_success( ovlap == SOME_OVERLAP, &success );
00337   if(out) *out  << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
00338 
00339   // //////////////////////////////////////////////////////
00340   // Test vector overlap (continuation of vector testing)
00341   //
00342   // ToDo: Finish this someday once you get it figured out.
00343 
00344   // ///////////////////////////
00345   // Test assignment operators
00346 
00347   if(out)
00348     *out  << "\n***\n*** assignment operators\n***\n";
00349 
00350   // DMatrixSlice
00351 
00352   if(out) *out << "\ngms1.bind(gm1());\n";
00353   gms1.bind(gm1());
00354 
00355   if(out) *out << "\ngms1 = 2.0;\n";
00356   gms1 = 2.0;
00357   if(out) *out << "gms1 =\n" << gms1;
00358   update_success( result = comp(gms1,2.0), &success );
00359   if(out) *out << "gms1 == 2.0 : " << result << std::endl; 
00360 
00361   if(out) *out << "\ngms1 = gms2;\n";
00362   gms1 = gms2;
00363   if(out) *out << "gms1 =\n" << gms1;
00364   update_success( result = comp(gms1,gms2), &success );
00365   if(out) *out << "gms1 == gms2 : " << result << std::endl; 
00366 
00367   // DMatrix
00368 
00369   if(out) *out << "\ngm1 = 3.0;\n";
00370   gm1 = 3.0;
00371   if(out) *out << "gm1 =\n" << gm1();
00372   update_success( result = comp(gm1(),3.0), &success );
00373   if(out) *out << "gm1 == 3.0 : " << result << std::endl; 
00374 
00375   if(out) *out << "\ngm1.resize(0,0); gm1 = gms2;\n";
00376   gm1.resize(0,0);
00377   gm1 = gms2;
00378   if(out) *out << "gm1 =\n" << gm1();
00379   update_success( result = comp(gm1(),gms2()), &success );
00380   if(out) *out << "gm1 == gms2 : " << result << std::endl; 
00381 
00382   if(out) *out << "\ngm1.resize(0,0); gm1 = gm4;\n";
00383   gm1.resize(0,0);
00384   gm1 = gm4;
00385   if(out) *out << "gm1 =\n" << gm1();
00386   update_success( result = comp(gm1(),gm4()), &success );
00387   if(out) *out << "gm1 == gm4 : " << result << std::endl; 
00388 
00389   } // end try
00390   catch( const std::exception& excpt ) {
00391     success = false;
00392     if(out)
00393       (*out)  << "\nError, a standard exception was thrown: "
00394           << typeName(excpt) << ": "
00395           << excpt.what() << std::endl; 
00396   }
00397   catch(...) {
00398     success = false;
00399     if(out)
00400       (*out)  << "\nError, an unknown exception was thrown\n";
00401   }
00402 
00403   if(out) {
00404     if(success)
00405       (*out)
00406         << "\n*** Congradulations, DMatrix and DMatrixSlice seem to check out. ***\n";
00407     else
00408       (*out)
00409         << "\n*** Oops, all of the tests for DMatrix and DMatrixSlice "
00410           "where not successful. ***\n";
00411   }
00412 
00413   return success;
00414 }
00415 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines