MOOCHO (Single Doxygen Collection) Version of the Day
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 // 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 <math.h>
00043 
00044 #include <ostream>
00045 #include <limits>
00046 
00047 #include "AbstractLinAlgPack_MatrixOpNonsingTester.hpp"
00048 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00049 #include "AbstractLinAlgPack_VectorSpace.hpp"
00050 #include "AbstractLinAlgPack_VectorMutable.hpp"
00051 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00052 #include "AbstractLinAlgPack_VectorOut.hpp"
00053 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00054 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00055 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00056 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00057 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00058 
00059 namespace AbstractLinAlgPack {
00060 
00061 MatrixOpNonsingTester::MatrixOpNonsingTester(
00062   ETestLevel       test_level
00063   ,EPrintTestLevel print_tests
00064   ,bool            dump_all
00065   ,bool            throw_exception
00066   ,size_type       num_random_tests
00067   ,value_type      warning_tol
00068   ,value_type      error_tol
00069   )
00070   :test_level_(test_level)
00071   ,print_tests_(print_tests)
00072   ,dump_all_(dump_all)
00073   ,throw_exception_(throw_exception)
00074   ,num_random_tests_(num_random_tests)
00075   ,warning_tol_(warning_tol)
00076   ,error_tol_(error_tol)
00077 {}
00078 
00079 bool MatrixOpNonsingTester::test_matrix(
00080   const MatrixOpNonsing   &M
00081   ,const char                     M_name[]
00082   ,std::ostream                   *out
00083   )
00084 {
00085   namespace rcp = MemMngPack;
00086   using BLAS_Cpp::no_trans;
00087   using BLAS_Cpp::trans;
00088   using BLAS_Cpp::left;
00089   using BLAS_Cpp::right;
00090   using AbstractLinAlgPack::sum;
00091   using AbstractLinAlgPack::dot;
00092   using AbstractLinAlgPack::Vp_StV;
00093   using AbstractLinAlgPack::assert_print_nan_inf;
00094   using AbstractLinAlgPack::random_vector;
00095   using LinAlgOpPack::V_StMtV;
00096   using LinAlgOpPack::V_MtV;
00097   using LinAlgOpPack::V_StV;
00098   using LinAlgOpPack::V_VpV;
00099   using LinAlgOpPack::Vp_V;
00100   
00101   // ToDo: Check the preconditions
00102   
00103   bool success = true, result, lresult;
00104   const value_type
00105     rand_y_l  = -1.0,
00106     rand_y_u  = 1.0,
00107     small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
00108     alpha     = 2.0;
00109   
00110   //
00111   // Perform the tests
00112   //
00113 
00114   if(out && print_tests() >= PRINT_BASIC)
00115     *out
00116       << "\nCheck: alpha*op(op(inv("<<M_name<<"))*op("<<M_name<<"))*v == alpha*v ...";
00117   if(out && print_tests() > PRINT_BASIC)
00118     *out << std::endl;
00119 
00120   VectorSpace::vec_mut_ptr_t
00121     v_c1 = M.space_cols().create_member(),
00122     v_c2 = M.space_cols().create_member(),
00123     v_r1 = M.space_rows().create_member(),
00124     v_r2 = M.space_rows().create_member();
00125 
00126   // Side of the matrix inverse 
00127   const BLAS_Cpp::Side    a_side[2]  = { BLAS_Cpp::left,     BLAS_Cpp::right };
00128   // If the matrices are transposed or not
00129   const BLAS_Cpp::Transp  a_trans[2] = { BLAS_Cpp::no_trans, BLAS_Cpp::trans };
00130 
00131   for( int side_i = 0; side_i < 2; ++side_i ) {
00132     for( int trans_i = 0; trans_i < 2; ++trans_i ) {
00133       const BLAS_Cpp::Side    t_side  = a_side[side_i];
00134       const BLAS_Cpp::Transp  t_trans = a_trans[trans_i];
00135       if(out && print_tests() >= PRINT_MORE)
00136         *out
00137           << "\n" << side_i+1 << "." << trans_i+1 << ") "
00138           << "Check: (t2 = "<<(t_side==left?"inv(":"alpha * ")<< M_name<<(t_trans==trans?"\'":"")<<(t_side==left?")":"")
00139           << " * (t1 = "<<(t_side==right?"inv(":"alpha * ")<<M_name<<(t_trans==trans?"\'":"")<<(t_side==right?")":"")
00140           << " * v)) == alpha * v ...";
00141       if(out && print_tests() > PRINT_MORE)
00142         *out << std::endl;
00143       result = true;
00144       VectorMutable
00145         *v  = NULL,
00146         *t1 = NULL,
00147         *t2 = NULL;
00148       if( (t_side == left && t_trans == no_trans) || (t_side == right && t_trans == trans) ) {
00149         // (inv(R)*R*v or R'*inv(R')*v
00150         v  = v_r1.get();
00151         t1 = v_c1.get();
00152         t2 = v_r2.get();
00153       }
00154       else {
00155         // (inv(R')*R'*v or R*inv(R)*v
00156         v  = v_c1.get();
00157         t1 = v_r1.get();
00158         t2 = v_c2.get();
00159       }
00160       for( int k = 1; k <= num_random_tests(); ++k ) {
00161         lresult = true;
00162         random_vector( rand_y_l, rand_y_u, v );
00163           if(out && print_tests() >= PRINT_ALL) {
00164           *out
00165             << "\n"<<side_i+1<<"."<<trans_i+1<<"."<<k<<") random vector " << k
00166             << " ( ||v||_1 / n = " << (v->norm_1() / v->dim()) << " )\n";
00167           if(dump_all() && print_tests() >= PRINT_ALL)
00168             *out << "\nv =\n" << *v;
00169         }
00170         // t1
00171         if( t_side == right ) {
00172           // t1 = inv(op(M))*v
00173           V_InvMtV( t1, M, t_trans, *v );
00174         }
00175         else {
00176           // t1 = alpha*op(M)*v
00177           V_StMtV( t1, alpha, M, t_trans, *v );
00178         }
00179         // t2
00180         if( t_side == left ) {
00181           // t2 = inv(op(M))*t1
00182           V_InvMtV( t2, M, t_trans, *t1 );
00183         }
00184         else {
00185           // t2 = alpha*op(M)*t1
00186           V_StMtV( t2, alpha, M, t_trans, *t1 );
00187         }
00188         const value_type
00189           sum_t2  = sum(*t2),
00190           sum_av  = alpha*sum(*v);
00191         assert_print_nan_inf(sum_t2, "sum(t2)",true,out);
00192         assert_print_nan_inf(sum_av, "sum(alpha*t1)",true,out);
00193         const value_type
00194           calc_err = ::fabs( ( sum_av - sum_t2 )
00195                      /( ::fabs(sum_av) + ::fabs(sum_t2) + small_num ) );
00196         if(out && print_tests() >= PRINT_ALL)
00197           *out
00198             << "\nrel_err(sum(alpha*v),sum(t2)) = "
00199             << "rel_err(" << sum_av << "," << sum_t2 << ") = "
00200             << calc_err << std::endl;
00201         if( calc_err >= warning_tol() ) {
00202           if(out && print_tests() >= PRINT_ALL)
00203             *out
00204               << std::endl
00205               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00206               << ", rel_err(sum(alpha*v),sum(t2)) = "
00207               << "rel_err(" << sum_av << "," << sum_t2 << ") = "
00208               << calc_err
00209               << " exceeded "
00210               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00211               << " = "
00212               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00213               << std::endl;
00214           if(calc_err >= error_tol()) {
00215             if(dump_all() && print_tests() >= PRINT_ALL) {
00216               *out << "\nalpha = " << alpha << std::endl;
00217               *out << "\nv =\n"    << *v;
00218               *out << "\nt1 =\n"   << *t2;
00219               *out << "\nt2 =\n"   << *t2;
00220             }
00221             lresult = false;
00222           }
00223         }
00224         if(!lresult) result = false;
00225       }
00226       if(!result) success = false;
00227       if( out && print_tests() == PRINT_MORE )
00228         *out << " : " << ( result ? "passed" : "failed" )
00229            << std::endl;
00230     }
00231   }
00232 
00233   if( out && print_tests() == PRINT_BASIC )
00234     *out << " : " << ( success ? "passed" : "failed" );
00235 
00236   return success;
00237 }
00238 
00239 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines