AbstractLinAlgPack_MatrixOpNonsingTester.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <math.h>
00030 
00031 #include <ostream>
00032 #include <limits>
00033 
00034 #include "AbstractLinAlgPack_MatrixOpNonsingTester.hpp"
00035 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00036 #include "AbstractLinAlgPack_VectorSpace.hpp"
00037 #include "AbstractLinAlgPack_VectorMutable.hpp"
00038 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00039 #include "AbstractLinAlgPack_VectorOut.hpp"
00040 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00041 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00042 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00043 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00044 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00045 
00046 namespace AbstractLinAlgPack {
00047 
00048 MatrixOpNonsingTester::MatrixOpNonsingTester(
00049   ETestLevel       test_level
00050   ,EPrintTestLevel print_tests
00051   ,bool            dump_all
00052   ,bool            throw_exception
00053   ,size_type       num_random_tests
00054   ,value_type      warning_tol
00055   ,value_type      error_tol
00056   )
00057   :test_level_(test_level)
00058   ,print_tests_(print_tests)
00059   ,dump_all_(dump_all)
00060   ,throw_exception_(throw_exception)
00061   ,num_random_tests_(num_random_tests)
00062   ,warning_tol_(warning_tol)
00063   ,error_tol_(error_tol)
00064 {}
00065 
00066 bool MatrixOpNonsingTester::test_matrix(
00067   const MatrixOpNonsing   &M
00068   ,const char                     M_name[]
00069   ,std::ostream                   *out
00070   )
00071 {
00072   namespace rcp = MemMngPack;
00073   using BLAS_Cpp::no_trans;
00074   using BLAS_Cpp::trans;
00075   using BLAS_Cpp::left;
00076   using BLAS_Cpp::right;
00077   using AbstractLinAlgPack::sum;
00078   using AbstractLinAlgPack::dot;
00079   using AbstractLinAlgPack::Vp_StV;
00080   using AbstractLinAlgPack::assert_print_nan_inf;
00081   using AbstractLinAlgPack::random_vector;
00082   using LinAlgOpPack::V_StMtV;
00083   using LinAlgOpPack::V_MtV;
00084   using LinAlgOpPack::V_StV;
00085   using LinAlgOpPack::V_VpV;
00086   using LinAlgOpPack::Vp_V;
00087   
00088   // ToDo: Check the preconditions
00089   
00090   bool success = true, result, lresult;
00091   const value_type
00092     rand_y_l  = -1.0,
00093     rand_y_u  = 1.0,
00094     small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
00095     alpha     = 2.0;
00096   
00097   //
00098   // Perform the tests
00099   //
00100 
00101   if(out && print_tests() >= PRINT_BASIC)
00102     *out
00103       << "\nCheck: alpha*op(op(inv("<<M_name<<"))*op("<<M_name<<"))*v == alpha*v ...";
00104   if(out && print_tests() > PRINT_BASIC)
00105     *out << std::endl;
00106 
00107   VectorSpace::vec_mut_ptr_t
00108     v_c1 = M.space_cols().create_member(),
00109     v_c2 = M.space_cols().create_member(),
00110     v_r1 = M.space_rows().create_member(),
00111     v_r2 = M.space_rows().create_member();
00112 
00113   // Side of the matrix inverse 
00114   const BLAS_Cpp::Side    a_side[2]  = { BLAS_Cpp::left,     BLAS_Cpp::right };
00115   // If the matrices are transposed or not
00116   const BLAS_Cpp::Transp  a_trans[2] = { BLAS_Cpp::no_trans, BLAS_Cpp::trans };
00117 
00118   for( int side_i = 0; side_i < 2; ++side_i ) {
00119     for( int trans_i = 0; trans_i < 2; ++trans_i ) {
00120       const BLAS_Cpp::Side    t_side  = a_side[side_i];
00121       const BLAS_Cpp::Transp  t_trans = a_trans[trans_i];
00122       if(out && print_tests() >= PRINT_MORE)
00123         *out
00124           << "\n" << side_i+1 << "." << trans_i+1 << ") "
00125           << "Check: (t2 = "<<(t_side==left?"inv(":"alpha * ")<< M_name<<(t_trans==trans?"\'":"")<<(t_side==left?")":"")
00126           << " * (t1 = "<<(t_side==right?"inv(":"alpha * ")<<M_name<<(t_trans==trans?"\'":"")<<(t_side==right?")":"")
00127           << " * v)) == alpha * v ...";
00128       if(out && print_tests() > PRINT_MORE)
00129         *out << std::endl;
00130       result = true;
00131       VectorMutable
00132         *v  = NULL,
00133         *t1 = NULL,
00134         *t2 = NULL;
00135       if( (t_side == left && t_trans == no_trans) || (t_side == right && t_trans == trans) ) {
00136         // (inv(R)*R*v or R'*inv(R')*v
00137         v  = v_r1.get();
00138         t1 = v_c1.get();
00139         t2 = v_r2.get();
00140       }
00141       else {
00142         // (inv(R')*R'*v or R*inv(R)*v
00143         v  = v_c1.get();
00144         t1 = v_r1.get();
00145         t2 = v_c2.get();
00146       }
00147       for( int k = 1; k <= num_random_tests(); ++k ) {
00148         lresult = true;
00149         random_vector( rand_y_l, rand_y_u, v );
00150           if(out && print_tests() >= PRINT_ALL) {
00151           *out
00152             << "\n"<<side_i+1<<"."<<trans_i+1<<"."<<k<<") random vector " << k
00153             << " ( ||v||_1 / n = " << (v->norm_1() / v->dim()) << " )\n";
00154           if(dump_all() && print_tests() >= PRINT_ALL)
00155             *out << "\nv =\n" << *v;
00156         }
00157         // t1
00158         if( t_side == right ) {
00159           // t1 = inv(op(M))*v
00160           V_InvMtV( t1, M, t_trans, *v );
00161         }
00162         else {
00163           // t1 = alpha*op(M)*v
00164           V_StMtV( t1, alpha, M, t_trans, *v );
00165         }
00166         // t2
00167         if( t_side == left ) {
00168           // t2 = inv(op(M))*t1
00169           V_InvMtV( t2, M, t_trans, *t1 );
00170         }
00171         else {
00172           // t2 = alpha*op(M)*t1
00173           V_StMtV( t2, alpha, M, t_trans, *t1 );
00174         }
00175         const value_type
00176           sum_t2  = sum(*t2),
00177           sum_av  = alpha*sum(*v);
00178         assert_print_nan_inf(sum_t2, "sum(t2)",true,out);
00179         assert_print_nan_inf(sum_av, "sum(alpha*t1)",true,out);
00180         const value_type
00181           calc_err = ::fabs( ( sum_av - sum_t2 )
00182                      /( ::fabs(sum_av) + ::fabs(sum_t2) + small_num ) );
00183         if(out && print_tests() >= PRINT_ALL)
00184           *out
00185             << "\nrel_err(sum(alpha*v),sum(t2)) = "
00186             << "rel_err(" << sum_av << "," << sum_t2 << ") = "
00187             << calc_err << std::endl;
00188         if( calc_err >= warning_tol() ) {
00189           if(out && print_tests() >= PRINT_ALL)
00190             *out
00191               << std::endl
00192               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00193               << ", rel_err(sum(alpha*v),sum(t2)) = "
00194               << "rel_err(" << sum_av << "," << sum_t2 << ") = "
00195               << calc_err
00196               << " exceeded "
00197               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00198               << " = "
00199               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00200               << std::endl;
00201           if(calc_err >= error_tol()) {
00202             if(dump_all() && print_tests() >= PRINT_ALL) {
00203               *out << "\nalpha = " << alpha << std::endl;
00204               *out << "\nv =\n"    << *v;
00205               *out << "\nt1 =\n"   << *t2;
00206               *out << "\nt2 =\n"   << *t2;
00207             }
00208             lresult = false;
00209           }
00210         }
00211         if(!lresult) result = false;
00212       }
00213       if(!result) success = false;
00214       if( out && print_tests() == PRINT_MORE )
00215         *out << " : " << ( result ? "passed" : "failed" )
00216            << std::endl;
00217     }
00218   }
00219 
00220   if( out && print_tests() == PRINT_BASIC )
00221     *out << " : " << ( success ? "passed" : "failed" );
00222 
00223   return success;
00224 }
00225 
00226 } // end namespace AbstractLinAlgPack

Generated on Tue Oct 20 12:51:43 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7