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