MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_TestGenMatrixOp.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 // 2/10/00: g++ 2.95.2 is giving me some trouble with update_success(...) in template
00030 // functions for some reason?
00031 //
00032 
00033 #include <iomanip>
00034 #include <ostream>
00035 #include <vector>
00036 #include <typeinfo>
00037 
00038 #include "DenseLinAlgPack_TestDenseLinAlgPack.hpp"
00039 #include "DenseLinAlgPack_DMatrixClass.hpp"
00040 #include "DenseLinAlgPack_DMatrixOp.hpp"
00041 #include "DenseLinAlgPack_DVectorOut.hpp"
00042 #include "DenseLinAlgPack_DMatrixOut.hpp"
00043 #include "DenseLinAlgPack_MatVecCompare.hpp"
00044 
00045 namespace {
00046 
00047 using DenseLinAlgPack::size_type;
00048 using DenseLinAlgPack::value_type;
00049 using DenseLinAlgPack::DVector;
00050 using DenseLinAlgPack::DVectorSlice;
00051 using DenseLinAlgPack::DMatrix;
00052 using DenseLinAlgPack::DMatrixSlice;
00053 
00054 // vs_lhs = alpha * op(M_rhs1) * vs_rhs2 + beta * vs_lhs
00055 template<class M_t>
00056 void test_MtV( M_t& M_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2
00057   , const DVectorSlice& expected_MtV, DVector* tmp1, DVector* tmp2
00058   , std::ostream* out, bool* success )
00059 {
00060   using BLAS_Cpp::rows;
00061   using DenseLinAlgPack::Vt_S;
00062   using DenseLinAlgPack::Vp_StV;
00063   using DenseLinAlgPack::Vp_StMtV;
00064   using DenseLinAlgPack::comp;
00065   using TestingHelperPack::update_success;
00066   
00067   // Check alpha = 1.0, 0.5.  Check beta = 0.0, 0.5, 1.0
00068   bool result = true;
00069   value_type
00070     alpha = 1.0,
00071     beta = 0.0;
00072   
00073   const size_type
00074     m = M_rhs1.rows(),
00075     n = M_rhs1.cols();
00076   
00077   tmp1->resize( rows(m,n,trans_rhs1) );
00078   tmp2->resize( tmp1->dim() );
00079 
00080   for( int i = 0; i < 2; ++i, alpha -= 0.5 ) {
00081     beta = 0.0;
00082     for( int j = 0; j < 3; ++j, beta += 0.5 ) {
00083       if(out)
00084         *out  << "vs_lhs = " << alpha << " * M_rhs1"
00085             << ( trans_rhs1 == BLAS_Cpp::trans ? '\'' : ' ' )
00086             << " * vs_rhs2 + " << beta << " * vs_lhs : ";
00087       // Initialize vs_lhs
00088       for( int k = 1; k <= tmp1->dim(); ++k )
00089         (*tmp1)(k) = (*tmp2)(k) = k;
00090       // Compute expected
00091       Vt_S( &(*tmp1)(), beta );
00092       Vp_StV( &(*tmp1)(), alpha, expected_MtV );
00093       // Computed
00094       Vp_StMtV( &(*tmp2)(), alpha, M_rhs1, trans_rhs1, vs_rhs2, beta );
00095       // Compare
00096       update_success( result = comp( *tmp1, *tmp2 ), success );
00097       if(out) *out << result << std::endl;
00098       if(out && !result) *out << "expected_vs_lhs =\n" << *tmp1
00099                   << "computed_vs_lhs =\n"<< *tmp2;
00100     }
00101   }
00102 }
00103 
00104 // Print a '\'' for trans and ' ' for no_trans
00105 inline char trans_char(BLAS_Cpp::Transp _trans) {
00106   return (_trans == BLAS_Cpp::trans ? '\'' : ' ' );
00107 }
00108 
00109 // Test DenseLinAlgPack compatable Level-3 BLAS (Mp_StMtM(...))
00110 template <class M_t>
00111 void test_MtM( M_t& B, BLAS_Cpp::Transp trans_B, const DMatrixSlice& I
00112   , BLAS_Cpp::Transp trans_I, std::ostream* out
00113   , DMatrix* C1, DMatrix* C2, bool* success )
00114 {
00115   using BLAS_Cpp::no_trans;
00116   using BLAS_Cpp::rows;
00117   using BLAS_Cpp::cols;
00118   using DenseLinAlgPack::assign;
00119   using DenseLinAlgPack::Mp_StM;
00120   using DenseLinAlgPack::Mp_StMtM;
00121   using DenseLinAlgPack::comp;
00122   using TestingHelperPack::update_success;
00123   
00124   value_type alpha = 1.0, beta = 0.0;
00125   bool result = true;
00126 
00127   for( int i = 1; i <= 3; ++i ) {
00128     // Set alpha + (beta)(0.5) = 1.0
00129     switch(i) {
00130       case 1:
00131         alpha = 1.0;
00132         beta = 0.0;
00133         break;
00134       case 2:
00135         alpha = 0.75;
00136         beta = 0.5;
00137         break;
00138       case 3:
00139         alpha = 0.5;
00140         beta = 1.0;
00141         break;
00142     }
00143     // Perform tests
00144     const char
00145       B_trans_chr = (trans_B==no_trans?' ':'\''),
00146       I_trans_chr = (trans_I==no_trans?' ':'\'');
00147     // expected result
00148     const size_type
00149       rows_B = rows( B.rows(), B.cols(), trans_B ),
00150       cols_B = cols( B.rows(), B.cols(), trans_B );
00151     C1->resize(rows_B,cols_B);
00152     (*C1) = 0.0;
00153     Mp_StM( &(*C1)(), 1.0, B, trans_B );
00154     // (left)
00155     if(out) *out << "C = " << alpha << " * B" << B_trans_chr
00156         << "* I" << I_trans_chr << " + " << beta << " * (0.5) * B" << B_trans_chr
00157         << " : ";
00158     C2->resize(rows_B,cols_B);
00159     (*C2) = 0.0;
00160     Mp_StM( &(*C2)(), 0.5, B, trans_B );
00161     size_type n = cols( B.rows(), B.cols(), trans_B );
00162     Mp_StMtM( &(*C2)(), alpha, B, trans_B, I(1,n,1,n), trans_I, beta );
00163     update_success( result = comp( (*C1)(), (*C2)() ), success );
00164     if(out) *out << result << std::endl;
00165     if( out && !result ) *out << "C1 =\n" << (*C1) << "C2 =\n" << (*C2);
00166     // (right)
00167     if(out) *out << "C = " << alpha << " * I" << I_trans_chr
00168         << "* B" << B_trans_chr << " + " << beta << " * (0.5) * B" << B_trans_chr
00169         << " : ";
00170     (*C2) = 0.0;
00171     Mp_StM( &(*C2)(), 0.5, B, trans_B );
00172     n = rows( B.rows(), B.cols(), trans_B );
00173     Mp_StMtM( &(*C2)(), alpha, I(1,n,1,n), trans_I, B, trans_B, beta );
00174     update_success( result = comp( (*C1)(), (*C2)() ), success );
00175     if(out) *out << result << std::endl;
00176     if( out && !result ) *out << "C1 =\n" << (*C1) << "C2 =\n" << (*C2);
00177 
00178   }
00179 }
00180 
00181 } // end namespace
00182 
00183 bool DenseLinAlgPack::TestingPack::TestGenMatrixOp(std::ostream* out)
00184 {
00185 
00186   using BLAS_Cpp::trans;
00187   using BLAS_Cpp::no_trans;
00188   using BLAS_Cpp::upper;
00189   using BLAS_Cpp::lower;
00190   using BLAS_Cpp::unit;
00191   using BLAS_Cpp::nonunit;
00192   using BLAS_Cpp::trans_to_bool;
00193   using DenseLinAlgPack::comp;
00194   using DenseLinAlgPack::sqrt_eps;
00195 
00196   bool success = true;
00197   bool result, result1, result2;
00198 
00199   if(out)
00200     *out  << "\n**************************************"
00201         << "\n*** Testing GenMatrixOp Operations ***"
00202         << "\n**************************************\n"
00203         << std::boolalpha << std::setw(10);
00204 
00205   // Arrays for looping through options.
00206 
00207   BLAS_Cpp::Uplo    a_uplo[2]     = { lower,    upper   };
00208   char        str_uplo[2][10]   = { "lower",  "upper"   };
00209 
00210   BLAS_Cpp::Diag    a_diag[2]     = { unit,   nonunit   };
00211   char        str_diag[2][10]   = { "unit",   "nonunit" };
00212 
00213   BLAS_Cpp::Transp  a_trans[2]      = { trans,    no_trans  };
00214 //  char        str_trans[2][10]  = { "trans",  "no_trans"  };
00215 
00216   try {
00217 
00218   const size_type
00219     m = 6,
00220     n = 8;
00221 
00222   const value_type
00223     ptr[m*n] =
00224        {  1.1,  2.1,  3.1,  4.1,  5.1,  6.1,
00225         1.2,  2.2,  3.2,  4.2,  5.2,  6.2,
00226         1.3,  2.3,  3.3,  4.3,  5.3,  6.3,
00227         1.4,  2.4,  3.4,  4.4,  5.4,  6.4,
00228         1.5,  2.5,  3.5,  4.5,  5.5,  6.5,
00229         1.6,  2.6,  3.6,  4.6,  5.6,  6.6,
00230         1.7,  2.7,  3.7,  4.7,  5.7,  6.7,
00231         1.8,  2.8,  3.8,  4.8,  5.8,  6.8 };
00232 
00233   DMatrix gm_lhs(m,n);
00234   const DMatrixSlice gms_rhs(const_cast<value_type*>(ptr),m*n,m,m,n);
00235 
00236   // ////////////////////////////////////////////////////////////////////////////////
00237   // Test Element-wise Assignment functions not already tested in TestGenMatrixClass
00238 
00239   if(out) *out << "\n***\n*** Test Element-wise Assignment functions not already "
00240           "tested in TestGenMatrixClass(...)\n***\n";
00241   
00242   // gms_lhs = op(gms_rhs)
00243   if(out) *out << "\ngms_lhs = op(gms_rhs)\n"
00244          << "\ngm_lhs.resize(n,m); assign( &gm_lhs(), gms_rhs, trans );\n";
00245   gm_lhs.resize(n,m);
00246   assign( &gm_lhs(), gms_rhs, trans );
00247   update_success( result = comp( gm_lhs(), no_trans, gms_rhs, trans ), &success );
00248   if(out && !result) *out << "gm_lhs =\n" << gm_lhs();
00249   if(out) *out << "gm_lhs == gms_rhs' : " << result << std::endl;
00250 
00251   // gms_lhs = op(gms_rhs)
00252   if(out) *out << "\ngm_lhs = op(gms_rhs)\n"
00253          << "\ngm_lhs.resize(1,1); assign( &gm_lhs, gms_rhs, trans );\n";
00254   gm_lhs.resize(1,1);
00255   assign( &gm_lhs, gms_rhs, trans );
00256   update_success( result = comp( gm_lhs(), no_trans, gms_rhs, trans ), &success );
00257   if(out && !result) *out << "gm_lhs =\n" << gm_lhs();
00258   if(out) *out << "gm_lhs == gms_rhs' : " << result << std::endl;
00259 
00260   // tri_ele_gms_lhs = alpha (elementwise)
00261 
00262   if(out) *out << "\ntri_ele_gms_lhs = alpha (elementwise)\n";
00263 
00264   if(out) *out << "\ngms_lhs.resize(m,m); gm_lhs = 1.0; assign(&nonconst_tri_ele(gm_lhs(),lower),2.0);\n";
00265   gm_lhs.resize(m,m);
00266   gm_lhs = 1.0;
00267   assign(&nonconst_tri_ele(gm_lhs(),lower),2.0);
00268   update_success( result1 = comp(tri_ele(gm_lhs(),lower),2.0), &success );
00269   update_success( result2 = comp(tri_ele(gm_lhs(1,m-1,2,m),upper),1.0), &success );
00270   if(out && (!result1 || !result2) )
00271     *out << "gm_lhs =\n" << gm_lhs();
00272   if(out) *out << "tri_ele(gm_lhs(),lower) == 2.0 : " << result1 << std::endl
00273          << "tri_ele(gm_lhs(1,m-1,2,m),upper) == 1.0 : " << result2 << std::endl;
00274 
00275   if(out) *out << "\nassign(&tri_ele(gm_lhs(1,m-1,2,m),upper),3.0);\n";
00276   assign(&nonconst_tri_ele(gm_lhs(1,m-1,2,m),upper),3.0);
00277   update_success( result1 = comp(tri_ele(gm_lhs(),lower),2.0), &success );
00278   update_success( result2 = comp(tri_ele(gm_lhs(1,m-1,2,m),upper),3.0), &success );
00279   if(out && (!result1 || !result2) )
00280     *out << "gm_lhs =\n" << gm_lhs();
00281   if(out) *out << "tri_ele(gm_lhs(),lower) == 2.0 : " << result1 << std::endl
00282          << "tri_ele(gm_lhs(1,m-1,2,m),upper) == 3.0 : " << result2 << std::endl;
00283 
00284   // tri_ele_gms_lhs = tri_ele_gms_rhs
00285 
00286   if(out) *out << "\ntri_ele_gms_lhs = tri_ele_gms_rhs\n"
00287          << "\nassign(&tri_ele(gm_lhs(2,m,1,m-1),lower),tri_ele(gms_rhs(1,m-1,2,m),upper));\n";
00288   assign(&nonconst_tri_ele(gm_lhs(2,m,1,m-1),lower),tri_ele(gms_rhs(1,m-1,2,m),upper));
00289   if(out) *out << "assign(&tri_ele(gm_lhs(1,m-1,2,m),upper),tri_ele(gms_rhs(2,m,1,m-1),lower));\n";
00290   assign(&nonconst_tri_ele(gm_lhs(1,m-1,2,m),upper),tri_ele(gms_rhs(2,m,1,m-1),lower));
00291   if(out) *out << "gm_lhs.diag() = gms_rhs.diag();\n";
00292   gm_lhs.diag() = gms_rhs.diag();
00293   update_success( result1 = comp( gm_lhs(1,m,1,m), no_trans, gms_rhs(1,m,1,m), trans ), &success );
00294   update_success( result2 = comp( tri_ele(gm_lhs(2,m,1,m-1),lower), tri_ele(gms_rhs(1,m-1,2,m),upper) ), &success );
00295   if(out && (!result1 || !result2) )
00296     *out << "gm_lhs =\n" << gm_lhs();
00297   if(out) *out << "gm_lhs(1,m,1,m) == gms_rhs(1,m,1,m)' : " << result1 << std::endl
00298          << "tri_ele(gm_lhs(2,m,1,m-1),lower) == tri_ele(gms_rhs(1,m-1,2,m),upper) : " << result2 << std::endl;
00299 
00300   // //////////////////////////////////////////
00301   // Test Element-wise Algebraic Operations
00302 
00303   if(out) *out << "\n***\n*** Test Element-wise Algebraic Operations\n***\n";
00304 
00305   // gms_lhs *= alpha
00306   if(out) *out << "\ngms_lhs *= alpha\n"
00307          << "\ngm_lhs = 1.0; Mt_S( &gm_lhs(), 2.0 );\n";
00308   gm_lhs = 1.0;
00309   Mt_S( &gm_lhs(), 2.0 );
00310   update_success( result = comp( gm_lhs(), 2.0 ), &success );
00311   if(out && !result)
00312     *out << "gm_lhs =\n" << gm_lhs();
00313   if(out) *out << "gm_lhs == 2.0: " << result << std::endl;
00314 
00315   // tri_lhs *= alpha
00316   {
00317     if(out) *out << "\ntri_lhs *= alpha\n"
00318            << "\ngm_lhs = 1.0;\nLet tm1 = tri_ele(gm_lhs(1,m,1,m),BLAS_Cpp::lower), "
00319             "tm2 = tri_ele(gm_lhs(1,m-1,2,m),BLAS_Cpp::upper)\n";
00320     gm_lhs = 1.0;
00321     DMatrixSliceTriEle
00322       tm1 = nonconst_tri_ele(gm_lhs(1,m,1,m),BLAS_Cpp::lower),
00323       tm2 = nonconst_tri_ele(gm_lhs(1,m-1,2,m),BLAS_Cpp::upper);
00324     if(out) *out << "Mt_S( &tm1, 2.0 );\n";
00325     Mt_S( &tm1, 2.0 );
00326     if(out) *out << "Mt_S( &tm2, 3.0 );\n";
00327     Mt_S( &tm2, 3.0 );
00328     update_success( result1 = comp( tm1, 2.0 ), &success );
00329     update_success( result2 = comp( tm2, 3.0 ), &success );
00330     if(out && (!result1 || !result2) )
00331       *out << "gm_lhs =\n" << gm_lhs();
00332     if(out) *out << "tm1 == 2.0 : " << result1 << std::endl
00333            << "tm2 == 3.0 : " << result2 << std::endl;
00334   }
00335 
00336   // tri_lhs += alpha * tri_rhs
00337   if(out) *out << "\ntri_lhs += alpha * tri_rhs\n"
00338             "gm_lhs = 1.0;\n";
00339   gm_lhs = 1.0;
00340   if(out) *out << "Mp_StM( &tri_ele(gm_lhs(2,m,1,m-1),lower), 2.0, tri_ele(gm_lhs(1,m-1,2,m),upper) );\n";
00341   Mp_StM( &nonconst_tri_ele(gm_lhs(2,m,1,m-1),lower), 2.0, tri_ele(gm_lhs(1,m-1,2,m),upper) );
00342   update_success( result1 = comp( tri_ele(gm_lhs(2,m,1,m-1),lower), 3.0 ), &success );
00343   update_success( result2 = comp( tri_ele(gm_lhs(1,m,1,m),upper), 1.0 ), &success );
00344   if(out && (!result1 || !result2) )
00345     *out << "gm_lhs =\n" << gm_lhs();
00346   if(out) *out << "tri_ele(gm_lhs(2,m,1,m-1),lower) == 3.0 : " << result1 << std::endl
00347          << "tri_ele(gm_lhs(1,m,1,m),upper) == 1.0 : " << result2 << std::endl;
00348 
00349   // LinAlgOpPack compatable (Mp_StM(...)).
00350   if(out) *out << "\n****** LinAlgOpPack compatable (Mp_StM(...))\n";
00351 
00352   // gms_lhs += alpha * op(gms_rhs)
00353   if(out) *out << "\ngms_lhs += alpha * op(gms_rhs)\n"
00354             "--------------------------\n"
00355           "gm_lhs.resize(m,n); gm_lhs = 0.0;"
00356           "Mp_StM( &gm_lhs(), 0.5, gm_rhs(), no_trans ); "
00357           "Mp_StM( &gm_lhs(), 0.5, gm_rhs(), no_trans );\n";
00358   gm_lhs.resize(m,n);
00359   gm_lhs = 0.0;
00360   Mp_StM( &gm_lhs(), 0.5, gms_rhs, no_trans );
00361   Mp_StM( &gm_lhs(), 0.5, gms_rhs, no_trans );
00362   update_success( result = comp( gm_lhs(), gms_rhs ), &success );
00363   if(out && !result)
00364     *out << "gm_lhs =\n" << gm_lhs();
00365   if(out) *out << "gm_lhs == gm_rhs : " << result << std::endl;
00366 
00367   // gms_lhs += alpha * op(sym_rhs)
00368   if(out) *out << "\ngms_lhs += alpha * op(sym_rhs)\n"
00369             "------------------------------\n";
00370   gm_lhs.resize(m,m);
00371   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00372     for(int i_trans = 0; i_trans < 2; ++i_trans) {
00373       const BLAS_Cpp::Uplo  _uplo = a_uplo[i_uplo];
00374       const BLAS_Cpp::Transp  _trans  = a_trans[i_trans];
00375       if(out)
00376         *out
00377           << "gms_lhs += alpha * sym(gms_rhs(1,m,1,m),"
00378           << str_uplo[i_uplo] << ")" << (_trans == trans ? '\'' : ' ' ) << " : ";
00379       gm_lhs = 0.0;
00380       const DMatrixSliceSym M = sym(gms_rhs(1,m,1,m),_uplo);
00381       Mp_StM( &gm_lhs(), 0.5, M, _trans );
00382       Mp_StM( &gm_lhs(), 0.5, M, _trans );
00383       update_success( result1 = comp( tri_ele(gm_lhs(),lower), tri_ele(M.gms(),_uplo) )
00384         , &success );
00385       update_success( result2 = comp( tri_ele(gm_lhs(),upper), tri_ele(M.gms(),_uplo) )
00386         , &success );
00387       if(out) *out << ( result1 && result2 ) << std::endl;
00388       if( out && ( !result1 || !result2 ) )
00389         *out << "gm_lhs =\n" << gm_lhs();
00390     }
00391   }}
00392 
00393   // gms_lhs += alpha * op(tri_rhs)
00394   if(out) *out << "\ngms_lhs += alpha * op(tri_rhs)\n"
00395             "------------------------------\n";
00396   gm_lhs.resize(m,m);
00397   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00398     for(int i_diag = 0; i_diag < 2; ++i_diag) {
00399       for(int i_trans = 0; i_trans < 2; ++i_trans) {
00400         const BLAS_Cpp::Uplo  _uplo = a_uplo[i_uplo];
00401         const BLAS_Cpp::Diag  _diag = a_diag[i_diag];
00402         const BLAS_Cpp::Transp  _trans  = a_trans[i_trans];
00403         if(out)
00404           *out
00405             << "gms_lhs += alpha * tri(gms_rhs(1,m,1,m),"
00406             << str_uplo[i_uplo] << "," << str_diag[i_diag] << ")"
00407             << (_trans == trans ? '\'' : ' ' ) << " : ";
00408         // compute
00409         gm_lhs = 0.0;
00410         const DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),_uplo,_diag);
00411         Mp_StM( &gm_lhs(), 0.5, M, _trans );
00412         Mp_StM( &gm_lhs(), 0.5, M, _trans );
00413         // test diagonal
00414         if( _diag == nonunit )
00415           result = comp( gm_lhs.diag(), gms_rhs.diag() );
00416         else
00417           result = comp( gm_lhs.diag(), 1.0 );
00418         update_success( result, &success ); 
00419         // test the rest
00420         const BLAS_Cpp::Uplo
00421           as_uplo = (   ( _uplo == lower && _trans == no_trans )
00422                 ||  ( _uplo == upper && _trans == trans )
00423                 ? lower : upper                     ),
00424           oth_as_uplo = ( as_uplo == lower ? upper : lower );
00425         const DMatrixSliceTriEle
00426           M_ele = tri_ele( ( _uplo==lower ? M.gms()(2,m,1,m-1) : M.gms()(1,m-1,2,m) )
00427                     , _uplo ),
00428           tri_reg_ele = tri_ele( ( as_uplo==lower ? gm_lhs(2,m,1,m-1) : gm_lhs(1,m-1,2,m) )
00429                     , as_uplo ),
00430           oth_tri_reg_ele = tri_ele( ( oth_as_uplo==lower ? gm_lhs(2,m,1,m-1) : gm_lhs(1,m-1,2,m) )
00431                     , oth_as_uplo );
00432         update_success( result1 = comp( tri_reg_ele, M_ele ), &success );
00433         update_success( result2 = comp( oth_tri_reg_ele, 0.0 ), &success );
00434         if(out) *out << ( result && result1 && result2 ) << std::endl;
00435         if( out && ( !result || !result1 || !result2 ) )
00436           *out << "gm_lhs =\n" << gm_lhs();
00437       }
00438     }
00439   }}
00440 
00441   // //////////////////////////////////////////
00442   // Level-2 BLAS
00443 
00444   if(out) *out << "\n***\n*** Level-2 BLAS \n***\n";
00445 
00446   DVector tmp1, tmp2, vs_rhs(n), expected_MtV;
00447   {for(int k = 1; k <= n; ++k) vs_rhs(k) = k; }
00448 
00449   // *** Triagnular matrices
00450   if(out) *out << "\n*** BLAS-2 Equivalent Triangular Matrix operations\n";
00451 
00452   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00453     for(int i_diag = 0; i_diag < 2; ++i_diag) {
00454       for(int i_trans = 0; i_trans < 2; ++i_trans) {
00455         if(out)
00456           *out << "\nLet M = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
00457              << "," << str_diag[i_diag] << ")"
00458              << (a_trans[i_trans] == trans ? '\'' : ' ' ) << std::endl;
00459         DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),a_uplo[i_uplo],a_diag[i_diag]);
00460         // v_lhs
00461         tmp1.resize(1);
00462         tmp2.resize(1);
00463         V_InvMtV( &tmp1, M, a_trans[i_trans], vs_rhs(1,m) );
00464         V_MtV( &tmp2, M, a_trans[i_trans], tmp1() );
00465         update_success( result = comp( tmp2, vs_rhs(1,m)), &success );
00466         if(out) *out << "(v_lhs)... M * (inv(M)*vs_rhs(1,m)) == vs_rhs(1,m) : "
00467                << result << std::endl;
00468         // vs_lhs
00469         V_InvMtV( &tmp1(), M, a_trans[i_trans], vs_rhs(1,m) );
00470         V_MtV( &tmp1(), M, a_trans[i_trans], tmp1() );
00471         update_success( result = comp( tmp1, vs_rhs(1,m)), &success );
00472         if(out) *out << "(vs_lhs)... M * (inv(M)*vs_rhs(1,m)) == vs_rhs(1,m) : "
00473                << result << std::endl;
00474       }
00475     }
00476   }}
00477 
00478   // *** LinAlgOpPack complient
00479   if(out) *out << "\n*** Test DenseLinAlgPack complient (Vp_StMtV(...))\n";
00480 
00481   // vs_lhs = alpha * op(gms_rhs1) * vs_rhs2 + beta * vs_lhs (xGEMV)
00482 
00483   if(out) *out << "\n****** General Matrix MtV\n";
00484 
00485   // no_trans
00486   if(out) *out << "\nvs_lhs = alpha * gms_rhs1 * vs_rhs2 + beta * vs_lhs (xGEMV)\n";
00487   expected_MtV.resize(m);
00488   expected_MtV = 0.0;
00489   {for(int i=1;i<=m;++i)
00490     for(int j=1;j<=n;++j)
00491       expected_MtV(i) += gms_rhs(i,j) * vs_rhs(j);
00492   }
00493   test_MtV( gms_rhs, no_trans, vs_rhs(1,n), expected_MtV, &tmp1, &tmp2
00494     , out, &success );
00495 
00496   // trans
00497   if(out) *out << "\nvs_lhs = alpha * gms_rhs1' * vs_rhs2 + beta * vs_lhs (xGEMV)\n";
00498   expected_MtV.resize(n);
00499   expected_MtV = 0.0;
00500   {for(int i=1;i<=n;++i)
00501     for(int j=1;j<=m;++j)
00502       expected_MtV(i) += gms_rhs(j,i) * vs_rhs(j);
00503   }
00504   test_MtV( gms_rhs, trans, vs_rhs(1,m), expected_MtV, &tmp1, &tmp2
00505     , out, &success );
00506 
00507   // vs_lhs = alpha * op(sym_rhs1) * vs_rhs2 + beta * vs_lhs (xSYMV)
00508 
00509   if(out) *out << "\n****** Symmetric Matrix MtV\n";
00510 
00511   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00512     for(int i_trans = 0; i_trans < 2; ++i_trans) {
00513       if(out)
00514         *out
00515           << "\nvs_lhs = alpha * sym(gms_rhs(1,m,1,m),"
00516           << str_uplo[i_uplo] << ")"
00517           << (a_trans[i_trans] == trans ? '\'' : ' ' )
00518           << " * vs_rhs2 + beta * vs_lhs (xSYMV)\n";
00519       // compute expected MtV
00520       expected_MtV.resize(m);
00521       expected_MtV = 0.0;
00522       switch(a_uplo[i_uplo]) {
00523         case lower: {
00524           for(int i=1;i<=m;++i)
00525             for(int j=1;j<=m;++j)
00526               expected_MtV(i)
00527                 += (j < i ? gms_rhs(i,j) : gms_rhs(j,i) ) * vs_rhs(j);
00528           break;
00529         }
00530         case upper: {
00531           for(int i=1;i<=m;++i)
00532             for(int j=1;j<=m;++j)
00533               expected_MtV(i)
00534                 += (j > i ? gms_rhs(i,j) : gms_rhs(j,i) ) * vs_rhs(j);
00535           break;
00536         }
00537       }
00538       test_MtV( sym(gms_rhs(1,m,1,m),a_uplo[i_uplo]), a_trans[i_trans]
00539         , vs_rhs(1,m), expected_MtV, &tmp1, &tmp2, out, &success );
00540     }
00541   }}
00542 
00543 
00544   // vs_lhs = alpha * op(tri_rhs1) * vs_rhs2 + beta * vs_lhs.
00545 
00546   if(out) *out << "\n****** Triangular Matrix MtV\n";
00547 
00548   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00549     for(int i_diag = 0; i_diag < 2; ++i_diag) {
00550       for(int i_trans = 0; i_trans < 2; ++i_trans) {
00551         const BLAS_Cpp::Uplo  _uplo = a_uplo[i_uplo];
00552         const BLAS_Cpp::Diag  _diag = a_diag[i_diag];
00553         const BLAS_Cpp::Transp  _trans  = a_trans[i_trans];
00554         if(out)
00555           *out
00556             << "\nvs_lhs = alpha * tri(gms_rhs(1,m,1,m),"
00557             << str_uplo[i_uplo] << "," << str_diag[i_diag] 
00558             << ")" << (_trans == trans ? '\'' : ' ' )
00559             << " * vs_rhs2 + beta * vs_lhs\n";
00560         DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),_uplo,_diag);
00561         //
00562         // Compute expected MtV
00563         //
00564         // Determine conceptually lower or upper triangular 
00565         const BLAS_Cpp::Uplo
00566           as_uplo = (   ( _uplo == lower && _trans == no_trans )
00567                 ||  ( _uplo == upper && _trans == trans )
00568                 ? lower : upper );
00569         // Compute expected
00570         expected_MtV.resize(m);
00571         expected_MtV = 0.0;
00572         switch(as_uplo) {
00573           case lower: {
00574             for(int i=1;i<=m;++i) {
00575               for(int j=1;j<i;++j) {
00576                 const int
00577                   _i = ( _trans == no_trans ? i : j ),
00578                   _j = ( _trans == no_trans ? j : i );
00579                 expected_MtV(i) += gms_rhs(_i,_j) * vs_rhs(j);
00580               }
00581               expected_MtV(i) += (_diag==unit?1.0:gms_rhs(i,i))*vs_rhs(i);
00582             }
00583             break;
00584           }
00585           case upper: {
00586             for(int i=1;i<=m;++i) {
00587               expected_MtV(i) += (_diag==unit?1.0:gms_rhs(i,i))*vs_rhs(i);
00588               for(int j=i+1;j<=m;++j) {
00589                 const int
00590                   _i = ( _trans == no_trans ? i : j ),
00591                   _j = ( _trans == no_trans ? j : i );
00592                 expected_MtV(i) += gms_rhs(_i,_j) * vs_rhs(j);
00593               }
00594             }
00595             break;
00596           }
00597         }
00598         test_MtV( tri(gms_rhs(1,m,1,m),_uplo,_diag), _trans, vs_rhs(1,m)
00599           , expected_MtV, &tmp1, &tmp2, out, &success );
00600       }
00601     }
00602   }}
00603 
00604   // *** Symmetric Matrix Updates.
00605 
00606   if(out) *out << "\n*** Symmetric Matrix Updates\n";
00607 
00608   // sym_lhs = alpha * vs_rhs * vs_rhs' + sym_lhs (BLAS xSYR).
00609   if(out) *out << "\n****** sym_lhs = alpha * vs_rhs * vs_rhs' + sym_lhs (BLAS xSYR)\n";
00610 
00611   gm_lhs.resize(m,m);
00612   gm_lhs = 0.0;
00613   DMatrix ex_gm_lhs(0.0,m,m);
00614   value_type alpha = 2.0;
00615     
00616   {for( int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00617     const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
00618     if(out)
00619       *out << "\nFor sym_lhs = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
00620          << ")\n";
00621     assign( &nonconst_tri_ele(gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
00622     assign( &nonconst_tri_ele(ex_gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
00623     // Compute expected
00624     for(int i = 1; i<=m; ++i) {
00625       for(int j = i; j<=m; ++j) { // upper triangular
00626         const int
00627           _i = (_uplo == upper ? i : j),  // adjust for lower triangular
00628           _j = (_uplo == upper ? j : i);
00629         ex_gm_lhs(_i,_j) += alpha * vs_rhs(i) * vs_rhs(j);
00630       }
00631     }
00632     // Compute update
00633     syr( alpha, vs_rhs(1,m), &nonconst_sym(gm_lhs(),_uplo) );
00634     // Compare
00635     update_success( result = comp( tri_ele(gm_lhs(),_uplo), tri_ele(ex_gm_lhs(),_uplo) )
00636       , &success );
00637     if( out && !result )
00638       *out << "gm_lhs =\n" << gm_lhs << "ex_gm_lhs =\n" << ex_gm_lhs();
00639     if(out) *out << "    gm_lhs == expected_gm_lhs : " << result << std::endl;
00640   }}
00641 
00642   // sym_lhs = alpha * vs_rhs1 * vs_rhs2' + alpha * vs_rhs2 * vs_rhs1' + sym_lhs (BLAS xSYR2).
00643   if(out) *out << "\n****** sym_lhs = alpha * vs_rhs1 * vs_rhs2' + alpha * vs_rhs2 * vs_rhs1' + sym_lhs (BLAS xSYR2)\n";
00644 
00645   gm_lhs = 0.0;
00646   ex_gm_lhs = 0.0;
00647   alpha = 2.0;
00648   DVector vs_rhs2 = vs_rhs.rev();
00649     
00650   {for( int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00651     const BLAS_Cpp::Uplo _uplo = a_uplo[i_uplo];
00652     if(out)
00653       *out << "\nFor sym_lhs = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
00654          << ")\n";
00655     assign( &nonconst_tri_ele(gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
00656     assign( &nonconst_tri_ele(ex_gm_lhs(),_uplo), tri_ele(gms_rhs(1,m,1,m),_uplo) );
00657     // Compute expected
00658     for(int i = 1; i<=m; ++i) {
00659       for(int j = i; j<=m; ++j) { // upper triangular
00660         const int
00661           _i = (_uplo == upper ? i : j),  // adjust for lower triangular
00662           _j = (_uplo == upper ? j : i);
00663         ex_gm_lhs(_i,_j) += alpha * (vs_rhs(i) * vs_rhs2(j) + vs_rhs2(i) * vs_rhs(j));
00664       }
00665     }
00666     // Compute update
00667     syr2( alpha, vs_rhs(1,m), vs_rhs2(1,m), &nonconst_sym(gm_lhs(),_uplo) );
00668     // Compare
00669     update_success( result = comp( tri_ele(gm_lhs(),_uplo), tri_ele(ex_gm_lhs(),_uplo) )
00670       , &success );
00671     if( out && !result )
00672       *out << "gm_lhs =\n" << gm_lhs() << "ex_gm_lhs =\n" << ex_gm_lhs();
00673     if(out) *out << "    gm_lhs == expected_gm_lhs : " << result << std::endl;
00674   }}
00675 
00676   // /////////////////////////////////////////
00677   // Level-3 BLAS
00678 
00679   if(out) *out << "\n***\n*** Level-3 BLAS\n***\n";
00680 
00681   // Triangular Matrices
00682   if(out) *out << "\n*** BLAS-2 Equivalent Triangular Matrix operations\n";
00683 
00684   DMatrix Tmp1, Tmp2;
00685 
00686   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00687     for(int i_diag = 0; i_diag < 2; ++i_diag) {
00688       for(int i_trans1 = 0; i_trans1 < 2; ++i_trans1) {
00689         for(int i_trans2 = 0; i_trans2 < 2; ++i_trans2) {
00690           const BLAS_Cpp::Transp
00691             _trans1 = a_trans[i_trans1],
00692             _trans2 = a_trans[i_trans2]; 
00693           if(out)
00694             *out << "\nLet M = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
00695                << "," << str_diag[i_diag] << ")"
00696                << trans_char(_trans1) << std::endl;
00697           DMatrixSliceTri M = tri(gms_rhs(1,m,1,m),a_uplo[i_uplo],a_diag[i_diag]);
00698           // gm_lhs (left)
00699           // gm_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) (left) (BLAS xTRMM).
00700           // gm_lhs = alpha * inv(op(tri_rhs1)) * op(gms_rhs2) (left) (BLAS xTRSM).
00701           Tmp1.resize(1,1);
00702           Tmp2.resize(1,1);
00703           M_StInvMtM( &Tmp1, 2.0, M, _trans1, gms_rhs(1,m,1,m), _trans2 );
00704           M_StMtM( &Tmp2, 0.5, M, _trans1, Tmp1(), no_trans );
00705           update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans2 )
00706             , &success );
00707           if(out) *out << "(gm_lhs,left)...0.5*M*(2*inv(M)*gms_rhs(1,m,1,m)"
00708                  << trans_char(_trans2) << ")"
00709                   " == gms_rhs(1,m,1,m)" << trans_char(_trans2) << " : "
00710                   << result << std::endl;
00711           // gms_lhs (left)
00712           // gms_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) (left) (BLAS xTRMM).
00713           // gms_lhs = alpha * inv(op(tri_rhs1)) * op(gms_rhs2) (left) (BLAS xTRSM).
00714           assign( &Tmp2, gms_rhs(1,m,1,m), _trans2 );
00715           M_StInvMtM( &Tmp1(), 2.0, M, _trans1, Tmp2(), no_trans );
00716           M_StMtM( &Tmp2, 0.5, M, _trans1, Tmp1(), no_trans );
00717           update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans2 )
00718             , &success );
00719           if(out) {
00720              *out
00721                << "(gms_lhs,left)...0.5*M*(2*inv(M)*gms_rhs(1,m,1,m)"
00722               << trans_char(_trans2) << ")"
00723                 " == gms_rhs(1,m,1,m)" << trans_char(_trans2) << " : "
00724               << result << std::endl;
00725             if(!result) {
00726               *out
00727                 << "\ngms_lhs =\n" << Tmp2
00728                 << "\ngms_rhs(1,m,1,m) =\n" << gms_rhs(1,m,1,m);
00729             }
00730           }
00731           // gm_lhs (right)
00732           // gm_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) (right) (BLAS xTRMM).
00733           // gm_lhs = alpha * op(gms_rhs1) * inv(op(tri_rhs2)) (right) (BLAS xTRSM).
00734           Tmp1.resize(1,1);
00735           Tmp2.resize(1,1);
00736           M_StMtInvM( &Tmp1, 2.0, gms_rhs(1,m,1,m), _trans1, M, _trans2 );
00737           M_StMtM( &Tmp2, 0.5, Tmp1(), no_trans, M, _trans2 );
00738           update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans1 )
00739             , &success );
00740           if(out) *out << "(gm_lhs,right)...5.0*(2.0*gms_rhs(1,m,1,m)"
00741                  << trans_char(_trans1) << "*inv(M))*M"
00742                   " == gms_rhs(1,m,1,m)" << trans_char(_trans1) << " : "
00743                   << result << std::endl;
00744           // gms_lhs (right)
00745           // gms_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) (right) (BLAS xTRMM).
00746           // gms_lhs = alpha * op(gms_rhs1) * inv(op(tri_rhs2)) (right) (BLAS xTRSM).
00747           assign( &Tmp2(), gms_rhs(1,m,1,m), _trans1 );
00748           M_StMtInvM( &Tmp2(), 2.0, Tmp2(), no_trans, M, _trans2 );
00749           M_StMtM( &Tmp2(), 0.5, Tmp2(), no_trans, M, _trans2 );
00750           update_success( result = comp( Tmp2(), no_trans, gms_rhs(1,m,1,m), _trans1 )
00751             , &success );
00752           if(out) *out << "(gms_lhs,right)...5.0*(2.0*gms_rhs(1,m,1,m)"
00753                  << trans_char(_trans1) << "*inv(M))*M"
00754                   " == gms_rhs(1,m,1,m)" << trans_char(_trans1) << " : "
00755                   << result << std::endl;
00756         }
00757       }
00758     }
00759   }}
00760 
00761   // *** LinAlgOpPack complient
00762   if(out) *out << "\n*** Test DenseLinAlgPack complient (Mp_StMtM(...))\n";
00763 
00764   DMatrix I(0.0,n,n);
00765   I.diag() = 1.0;
00766 
00767   // ****** Rectangular Matrices 
00768   if(out) *out << "\n****** Rectangular MtM (BLAS xGEMV)\n";
00769 
00770   // gms_lhs = alpha * op(gms_rhs1) * op(gms_rhs2) + beta * gms_lhs (BLAS xGEMV).
00771   {for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
00772     for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
00773       const BLAS_Cpp::Transp
00774         _trans_B = a_trans[i_trans_B],
00775         _trans_I = a_trans[i_trans_I];
00776       if(out) *out << "\nLet B = gms_rhs\n";
00777       test_MtM( gms_rhs, _trans_B, I, _trans_I, out, &Tmp1, &Tmp2, &success );
00778     }
00779   }}
00780 
00781   // ****** Symmetric Matrices
00782   if(out) *out << "\n****** Symmetric MtM\n";
00783 
00784   // gms_lhs = alpha * op(sym_rhs1) * op(gms_rhs2) + beta * gms_lhs (left) (BLAS xSYMM).
00785   // gms_lhs = alpha * op(gms_rhs1) * op(sym_rhs2) + beta * gms_lhs (right) (BLAS xSYMM).
00786   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00787     for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
00788       for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
00789         const BLAS_Cpp::Uplo
00790           _uplo = a_uplo[i_uplo];
00791         const BLAS_Cpp::Transp
00792           _trans_B = a_trans[i_trans_B],
00793           _trans_I = a_trans[i_trans_I]; 
00794         if(out) *out << "\nLet B = sym(gms_rhs(1,m,1,m)," << str_uplo[i_uplo] << ")\n";
00795         test_MtM( sym(gms_rhs(1,m,1,m),_uplo), _trans_B, I, _trans_I, out
00796           , &Tmp1, &Tmp2, &success );
00797       }
00798     }
00799   }}
00800 
00801   // ****** Triangular Matrices
00802   if(out) *out << "\n****** Triangular MtM\n";
00803 
00804   // gms_lhs = alpha * op(tri_rhs1) * op(gms_rhs2) + beta * gms_lhs (left).
00805   // gms_lhs = alpha * op(gms_rhs1) * op(tri_rhs2) + beta * gms_lhs (right).
00806 
00807   {for(int i_uplo = 0; i_uplo < 2; ++i_uplo) {
00808     for(int i_diag = 0; i_diag < 2; ++i_diag) {
00809       for(int i_trans_B = 0; i_trans_B < 2; ++i_trans_B) {
00810         for(int i_trans_I = 0; i_trans_I < 2; ++i_trans_I) {
00811           const BLAS_Cpp::Uplo
00812             _uplo = a_uplo[i_uplo];
00813           const BLAS_Cpp::Diag
00814             _diag = a_diag[i_diag];
00815           const BLAS_Cpp::Transp
00816             _trans_B = a_trans[i_trans_B],
00817             _trans_I = a_trans[i_trans_I]; 
00818           if(out) *out << "\nLet B = tri(gms_rhs(1,m,1,m)," << str_uplo[i_uplo]
00819                  << "," << str_diag[i_diag] << ")\n";
00820           test_MtM( tri(gms_rhs(1,m,1,m),_uplo,_diag), _trans_B, I, _trans_I, out
00821             , &Tmp1, &Tmp2, &success );
00822         }
00823       }
00824     }
00825   }}
00826 
00827   // *** Symmetric Matrix updating
00828   if(out) *out << "\n*** Symmetric Matrix updating\n";
00829 
00830   if(out) *out << "\nWarning! Not Tested!\n";
00831 
00832   // sym_lhs = alpha * op(gms_rhs) * op(gms_rhs')  + beta * sym_lhs (BLAS xSYRK).
00833   // sym_lhs = alpha * op(gms_rhs1) * op(gms_rhs2') + alpha * op(gms_rhs2) * op(gms_rhs1') + beta * sym_lhs (BLAS xSYR2K)
00834 
00835   } // end try
00836   catch( const std::exception& excpt ) {
00837     success = false;
00838     if(out)
00839       (*out)  << "\nError, a standard exception was thrown: "
00840           << typeName(excpt) << ": "
00841           << excpt.what() << std::endl; 
00842   }
00843   catch(...) {
00844     success = false;
00845     if(out)
00846       (*out)  << "\nError, an unknown exception was thrown\n";
00847   }
00848 
00849   if(out) {
00850     if(success)
00851       (*out)
00852         << "\n*** Congradulations, GenMatrixOp operations seem to check out. ***\n";
00853     else
00854       (*out)
00855         << "\n*** Oops, all of the tests for GenMatrixOp operations "
00856           "where not successful. ***\n";
00857   }
00858 
00859 
00860   return success;
00861 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines