AbstractLinAlgPack_VectorSpaceTester.cpp

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 
00033 #include "AbstractLinAlgPack_VectorSpaceTester.hpp"
00034 #include "AbstractLinAlgPack_VectorSpace.hpp"
00035 #include "AbstractLinAlgPack_VectorMutable.hpp"
00036 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00037 #include "AbstractLinAlgPack_VectorOut.hpp"
00038 #include "TestingHelperPack_update_success.hpp"
00039 #include "Teuchos_TestForException.hpp"
00040 
00041 namespace {
00042 template< class T >
00043 inline
00044 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00045 template< class T >
00046 inline
00047 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00048 } // end namespace
00049 
00050 namespace AbstractLinAlgPack {
00051 
00052 VectorSpaceTester::VectorSpaceTester(
00053   bool         print_all_tests
00054   ,bool        print_vectors
00055   ,bool        throw_exception
00056   ,size_type   num_random_tests
00057   ,value_type  warning_tol
00058   ,value_type  error_tol
00059   )
00060   :print_all_tests_(print_all_tests)
00061   ,print_vectors_(print_vectors)
00062   ,throw_exception_(throw_exception)
00063   ,num_random_tests_(num_random_tests)
00064   ,warning_tol_(warning_tol)
00065   ,error_tol_(error_tol)
00066 {}
00067 
00068 bool VectorSpaceTester::check_vector_space(
00069   const VectorSpace &space
00070   ,std::ostream     *out
00071   ) const
00072 {
00073   namespace rcp = MemMngPack;
00074 
00075   bool success = true, result = false;
00076 
00077   try {
00078 
00079   // Adapted from test_my_vector_library(...)
00080 
00081   const value_type 
00082     rand_l = -10.0,
00083     rand_u = +10.0;
00084 
00085   typedef VectorMutable::vec_ptr_t       vec_ptr_t;
00086   typedef VectorMutable::vec_mut_ptr_t   vec_mut_ptr_t;
00087 
00088   // Create three random non-mutable vectors
00089   vec_ptr_t            v_array[3];
00090   const Vector*  v[3];
00091   {for(int k = 0; k < 3; ++k) {
00092     vec_mut_ptr_t  r = space.create_member();
00093     random_vector( rand_l, rand_u, r.get() );
00094     v_array[k] = r;
00095     v[k] = v_array[k].get();
00096   }}
00097 
00098   // Create six random mutable vectors
00099   vec_mut_ptr_t          z_array[6];
00100   VectorMutable*   z[6];
00101   {for(int k = 0; k < 6; ++k) {
00102     random_vector( rand_l, rand_u, (z_array[k]= space.create_member()).get() );
00103     z[k] = z_array[k].get();
00104   }}
00105 
00106   if(out && print_all_tests())
00107     *out << std::boolalpha
00108        << "\n**************************************************"
00109        << "\n*** VectorSpaceTester::check_vector_space(...) ***"
00110        << "\n**************************************************\n";
00111 
00112   RTOp_index_type
00113     n = space.dim();
00114   RTOp_value_type
00115     err     = 0.0,
00116     max_err = 0.0,
00117     sum_err = 0.0;
00118   char z_name[20], v_name[20];
00119 
00120   // Print the vector space dimension
00121   if(out && print_all_tests())
00122     *out << "\nspace->dim() = " << n << std::endl;
00123 
00124   // Print the initial vectors
00125   if(out && print_vectors()) {
00126     *out << "\n*** Printing the immutable vectors\n";
00127     {for(int j = 0; j < 3; ++j) {
00128       sprintf( v_name, "v[%d]", j );
00129       *out << std::endl << v_name << " : " << typeName(*v[j]) << std::endl
00130          << *v[j];
00131     }}
00132     *out << "\n*** Printing the mutable vectors\n";
00133     {for(int k = 0; k < 6; ++k) {
00134       sprintf( z_name, "z[%d]", k );
00135       *out << std::endl << z_name << " : " << typeName(*z[k]) << std::endl
00136          << *z[k];
00137     }}
00138   }
00139 
00141   if(out && print_all_tests())
00142     *out << "\n*** Testing the obvious assertions\n";
00143   {
00144     {for(int k = 0; k < 6; ++k) {
00145       const Vector *vec = NULL;
00146       if( k < 3 ) {
00147         sprintf( v_name, "v[%d]", k );
00148         vec = v[k];
00149       }
00150       else {
00151         sprintf( v_name, "z[%d]", k-3 );
00152         vec = z[k-3];
00153       }
00154       
00155       result = vec->space().is_compatible(space);
00156       if(out && (print_all_tests() || !result))
00157         *out << "\ncheck: " << v_name << "->space().is_compatible(space) : " << result << std::endl;
00158       check_test( result ? 0.0 : -10.0 , out, &success );
00159       
00160       result = space.is_compatible(vec->space());
00161       if(out && (print_all_tests() || !result))
00162         *out << "check: space.is_compatible(" << v_name << "->space()) : " << result << std::endl;
00163       check_test( result ? 0.0 : -10.0 , out, &success );
00164       
00165       err = vec->dim() - space.dim();
00166       if(out && (print_all_tests() || !result))
00167         *out << "check: " << v_name << "->dim() - space.dim() = " << vec->dim() << " - "
00168            << space.dim() << " = " << err << std::endl;
00169       check_test( err , out, &success );
00170 
00171       result = vec->nz() <= space.dim();
00172       if(out && (print_all_tests() || !result))
00173         *out << "check: " << v_name << "->nz() <= space.dim() = " << vec->nz() << " <= " << space.dim()
00174            << " : " << result << std::endl;
00175       check_test( result ? 0.0 : -10.0 , out, &success );
00176 
00177     }}
00178   }
00179 
00181   if(out && print_all_tests())
00182     *out << "\n*** Testing scalar assignment and element access methods\n";
00183   {
00184     const index_type k = 0;
00185     sprintf( z_name, "z[%d]", k );
00186     if(out && print_all_tests())
00187       *out << "\n0.0 -> " << z_name << std::endl;
00188     *z[k] = 0.0;
00189     if(out && print_vectors())
00190       *out << std::endl << z_name << " =\n" << *z[k];
00191     {for(int r = 0; r < num_random_tests(); ++r) {
00192       std::srand( n / (1+r) + r ); // This is very important in a parallel program!
00193       const RTOp_index_type
00194         i = my_min( n, my_max( (RTOp_index_type)(((double)std::rand() / RAND_MAX) * n + 1.0), 1 ) );
00195       const RTOp_value_type
00196         val = 10.0;
00197       
00198       if(out && print_all_tests())
00199         *out << std::endl << z_name << ".set_ele("<<i<<","<<val<<")\n";
00200       z[k]->set_ele(i,val);
00201       if(out && print_vectors())
00202         *out << std::endl << z_name << " =\n" << *z[k];
00203       RTOp_value_type
00204         val_get = z[k]->get_ele(i);
00205       err = val - val_get;
00206       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00207         *out << "check: " << val << " - " << z_name << ".get_ele(" << i << ") = "
00208            << val << " - " << val_get << " = " << err << std::endl;
00209       check_test(err,out,&success);
00210       
00211       RTOp_value_type
00212         z_k_sum = sum(*z[k]);
00213       err = val - z_k_sum;
00214       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00215         *out << "check: " << val << " - sum(" << z_name << ") = "
00216            << val << " - " << z_k_sum << " = " << err << std::endl;
00217       check_test(err,out,&success);
00218 
00219       if(out && print_all_tests())
00220         *out << z_name << ".set_ele("<<i<<",0.0)\n";
00221       z[k]->set_ele(i,0.0);
00222       if(out && print_vectors())
00223         *out << std::endl << z_name << " =\n" << *z[k];
00224       z_k_sum = sum(*z[k]);
00225       err = z_k_sum;
00226       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00227         *out << "check: sum(" << z_name << ") = " << z_k_sum << std::endl;
00228       check_test(err,out,&success);
00229     }}
00230   }
00231 
00233   if(out && print_all_tests())
00234     *out << "\n*** Testing vector assignment\n";
00235   {
00236     {for( int k = 0; k < 3; ++k ) {
00237       sprintf( z_name, "z[%d]", k );
00238       const AbstractLinAlgPack::Vector *vec = NULL;
00239       sprintf( v_name, "v[%d]", k );
00240       if(out && print_all_tests())
00241         *out << "\n" << v_name << " -> " << z_name << "\n";
00242       *z[k] = *v[k];
00243       const RTOp_value_type
00244         sum_z_k = sum(*z[k]),
00245         sum_v_k = sum(*v[k]);
00246       err = (sum_z_k - sum_v_k)/n;
00247       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00248         *out << "check: (sum(" << z_name << ") - sum(" << v_name << "))/n = ("
00249            << sum_z_k << " - " << sum_v_k << ")/" << n << " = " << err << std::endl;
00250       check_test(err,out,&success);
00251     }}
00252   }
00253 
00254   /*
00255 
00257   if(out && print_all_tests())
00258     *out << "\n*** Testing sub-vector and sub-space access\n";
00259   {
00260     const index_type k = 0;
00261     sprintf( z_name, "z[%d]", k );
00262     if(out && print_all_tests())
00263       *out << "\n0.0 -> " << z_name << std::endl;
00264     *z[k] = 0.0;
00265     if(out && print_vectors())
00266       *out << std::endl << z_name << " =\n" << *z[k];
00267     {for(int r = 0; r < num_random_tests(); ++r) {
00268       RTOp_index_type
00269         i1 = my_min( n, (RTOp_index_type)(((double)rand() / RAND_MAX) * n + 1) ),
00270         i2 = my_min( n, (RTOp_index_type)(((double)rand() / RAND_MAX) * n + 1) );
00271       if( i1 > i2 ) std::swap( i1, i2 );
00272       RTOp_index_type
00273         sub_vec_dim = i2-i1+1;
00274       const RTOp_value_type
00275         val = 10.0;
00276       
00277       if(out && print_all_tests())
00278         *out << "\nsub_vec_mut = " << z_name
00279            << ".sub_view(" << i1 << "," << i2 << ")\n";
00280       VectorMutable::vec_mut_ptr_t
00281         sub_vec_mut = z[k]->sub_view(i1,i2);
00282       if(out && print_all_tests())
00283         *out << "sub_vec_mut = " << val << std::endl;
00284       *sub_vec_mut = val;
00285       if(out && print_vectors())
00286         *out << std::endl << z_name << " =\n" << *z[k];
00287 
00288       err = sub_vec_dim - sub_vec_mut->dim();
00289       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()))
00290         *out << "check: ("<<i2<<"-"<<i1<<"+1) - sub_vec_mut.dim() = "
00291            << sub_vec_dim <<" - " << sub_vec_mut->dim() << " = " << err << std::endl;
00292       check_test(err,out,&success);
00293 
00294       if(out && print_all_tests())
00295         *out << "\nsub_space = space.sub_space(" << i1 << "," << i2 << ")\n";
00296       VectorSpace::space_ptr_t
00297         sub_space = space.sub_space(i1,i2);
00298 
00299       result = sub_vec_mut->space().is_compatible(*sub_space);
00300       if(out && (print_all_tests() || !result))
00301         *out << "check: sub_vec_mut->space().is_compatible(*sub_space) : " << result << std::endl;
00302       check_test( result ? 0.0 : -10.0 , out, &success );
00303       
00304       result = sub_space->is_compatible(sub_vec_mut->space());
00305       if(out && (print_all_tests() || !result))
00306         *out << "check: sub_space->is_compatible(*sub_vec_mut->space()) : " << result << std::endl;
00307       check_test( result ? 0.0 : -10.0 , out, &success );
00308         
00309       RTOp_value_type
00310         expected_sum = val*sub_vec_dim,
00311         z_k_sum = sum( *z[k] );
00312       err = (expected_sum - z_k_sum)/sub_vec_mut->dim();
00313       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()))
00314         *out << "check: ("<<val<<"*("<<i2<<"-"<<i1<<"+1) - sum("<<z_name<<"))/"<<sub_vec_dim
00315            << " = ("<<expected_sum<<" - "<<z_k_sum<<")/"<<sub_vec_dim<<" = "<<err<<std::endl;
00316       check_test(err,out,&success);
00317          
00318       if(out && print_all_tests())
00319         *out << "sub_vec = "<<z_name<<"{const}.sub_view("<<i1<<","<<i2<<")\n";
00320       Vector::vec_ptr_t
00321         sub_vec = static_cast<const Vector*>(z[k])->sub_view(i1,i2);
00322 
00323       err = sub_vec_dim - sub_vec->dim();
00324       if(out && print_all_tests())
00325         *out << "check: ("<<i2<<"-"<<i1<<"+1) - sub_vec.dim() = "
00326            << sub_vec_dim << " - " << sub_vec->dim() << " = " << err << std::endl;
00327       check_test(err,out,&success);
00328         
00329       expected_sum = val*sub_vec_dim;
00330       z_k_sum = sum(*sub_vec);
00331       err = (expected_sum - z_k_sum)/sub_vec_mut->dim();
00332       if(out && print_all_tests())
00333         *out << "check: ("<<val<<"*("<<i2<<"-"<<i1<<"+1) - sum(sub_vec))/"<<sub_vec_dim
00334            << " = ("<<expected_sum<<" - "<<z_k_sum<<")/"<<sub_vec_dim<<" = "<<err<<std::endl;
00335 
00336       if(out && print_all_tests())
00337         *out << "sub_vec_mut = 0.0\n";
00338       *sub_vec_mut = 0.0;
00339       if(out && print_vectors())
00340         *out << std::endl << z_name << " =\n" << *z[k];
00341     }}
00342   }
00343 
00345   if(out && print_all_tests())
00346     *out << "\n*** Testing explicit sub-vector access\n";
00347   {
00348     const index_type k = 0;
00349     const Vector
00350       &v_from = *v[k];
00351     sprintf( v_name, "v[%d]", k );
00352     VectorMutable
00353       &z_to = *z[k];
00354     sprintf( z_name, "z[%d]", k );
00355     
00356     if(out && print_all_tests())
00357       *out << "\n0.0 -> " << z_name << std::endl;
00358     *z[k] = 0.0;
00359     if(out && print_vectors())
00360       *out << std::endl << z_name << " =\n" << *z[k];
00361 
00362     {for(int r = 0; r < num_random_tests(); ++r) {
00363       const RTOp_index_type // Get random small sub-vectors so parallel efficiency will be good
00364         i1 = my_min( n, (RTOp_index_type)(((double)rand() / RAND_MAX) * n + 1) ),
00365         i2 = my_min( (RTOp_index_type)(i1 + ((double)rand() / RAND_MAX) * 9), n );
00366       const RTOp_index_type
00367         sub_vec_dim = i2-i1+1;
00368 
00369       if(out && print_all_tests())
00370         *out << std::endl << v_name << ".get_sub_vector(Rang1D("<<i1<<","<<i2<<"),SPARSE,&sub_vec)\n";
00371       RTOpPack::SubVector sub_vec;
00372       v_from.get_sub_vector(Range1D(i1,i2),&sub_vec); 
00373 
00374       err = sub_vec_dim - sub_vec.subDim();
00375       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()))
00376         *out << "check: ("<<i2<<"-"<<i1<<"+1) - sub_vec.subDim() = "
00377            << sub_vec_dim << " - " << sub_vec.subDim() << "  = " << err << std::endl;
00378       check_test(err,out,&success);
00379 
00380       if(out && print_all_tests())
00381         *out << z_name << ".set_sub_vector(sub_vec)\n";
00382       RTOpPack::SparseSubVector spc_sub_vec( sub_vec );
00383       z_to.set_sub_vector(spc_sub_vec);
00384       if(out && print_vectors())
00385         *out << std::endl << z_name << " =\n" << z_to;
00386       
00387       const RTOp_value_type
00388         v_sub_vec_sum = sum(*v_from.sub_view(i1,i2)),
00389         z_sum         = sum(z_to);
00390       err = (v_sub_vec_sum - z_sum)/sub_vec_dim;
00391       if(out && (print_all_tests() || ::fabs(err) >= warning_tol()))
00392         *out << "check: (sum(*"<<v_name<<".sub_view("<<i1<<","<<i2<<"))-sum("<<z_name<<"))/sub_vec.subDim()"
00393           " = ("<<v_sub_vec_sum<<"-"<<z_sum<<")/"<<sub_vec_dim<<" = "<<err<<std::endl;
00394 
00395       if(out && print_all_tests())
00396         *out << v_name << ".free_sub_vector(&sub_vec)\n";
00397       v_from.free_sub_vector(&sub_vec); 
00398         
00399       if(out && print_all_tests())
00400         *out << "*" << z_name<<".sub_view("<<i1<<","<<i2<<") = 0.0\n";
00401       *z_to.sub_view(i1,i2) = 0.0;
00402       if(out && print_vectors())
00403         *out << std::endl << z_name << " =\n" << z_to;
00404 
00405     }}
00406   }
00407 
00408   */
00409 
00411   if(out && print_all_tests())
00412     *out << "\n*** Testing norms\n";
00413   if(n > 1) {
00414     const index_type k = 0;
00415     sprintf( z_name, "z[%d]", k );
00416 
00417     const value_type val1 = -2.0, val2 = 3.0;
00418     const index_type i_mid = n/2;
00419     
00420     if(out && print_all_tests())
00421       *out << std::endl << val1 << " -> *" << z_name << ".sub_view(1,"<<i_mid<<")\n";
00422     *z[k]->sub_view(1,i_mid) = val1;
00423     if(out && print_all_tests())
00424       *out << val2 << " -> *" << z_name << ".sub_view("<<i_mid+1<<","<<n<<")\n";
00425     *z[k]->sub_view(i_mid+1,n) = val2;
00426     if(out && print_vectors())
00427       *out << std::endl << z_name << " =\n" << *z[k] << std::endl;
00428 
00429     value_type
00430       norm_1        = z[k]->norm_1(),
00431       expect_norm_1 = (::fabs(val1)*(i_mid) + ::fabs(val2)*(n-i_mid));
00432     err = (norm_1 - expect_norm_1)/n;
00433     if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00434       *out << "check: (" << z_name << "->norm_1() - |"<<val1<<"|*("<<i_mid<<")+"
00435          << "|"<<val2<<"|*("<<n<<"-"<<i_mid<<"))/"<<n
00436          <<" = ("<<norm_1<<" - "<<expect_norm_1<<")/"<<n<<" = "<<err<<std::endl;
00437     check_test(err,out,&success);
00438 
00439     value_type
00440       norm_2        = z[k]->norm_2(),
00441       expect_norm_2 = ::sqrt(val1*val1*(i_mid) + val2*val2*(n-i_mid));
00442     err = (norm_2 - expect_norm_2)/n;
00443     if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00444       *out << "check: (" << z_name << "->norm_2() - ("<<val1<<")^2*("<<i_mid<<")+"
00445          << "("<<val2<<")^2*("<<n<<"-"<<i_mid<<"))/"<<n
00446          <<" = ("<<norm_2<<" - "<<expect_norm_2<<")/"<<n<<" = "<<err<<std::endl;
00447     check_test(err,out,&success);
00448 
00449     value_type
00450       norm_inf        = z[k]->norm_inf(),
00451       expect_norm_inf = my_max(::fabs(val1),::fabs(val2));
00452     err = (norm_inf - expect_norm_inf)/n;
00453     if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00454       *out << "check: (" << z_name << "->norm_inf() - max(|"<<val1<<"|,"
00455          << "|"<<val2<<"|)/"<<n<<" = ("<<norm_inf<<" - "<<expect_norm_inf<<")/"<<n
00456          <<" = "<<err<<std::endl;
00457     check_test(err,out,&success);
00458 
00459   }
00460   else {
00461     if(out && print_all_tests())
00462       *out << "space.dim() <= 1, can't test the norms...\n";
00463   }
00464   
00466   if(out && print_all_tests())
00467     *out << "\n*** Testing clone() method\n";
00468   {
00469     if(out && print_all_tests())
00470       *out << "\n(*vec = space.create_member()) = v[0]\n";
00471     VectorSpace::vec_mut_ptr_t
00472       vec = space.create_member();
00473     *vec = *v[0];
00474     if(out && print_all_tests())
00475       *out << "vec_clone = vec->clone()\n";
00476     VectorSpace::vec_mut_ptr_t
00477       vec_clone = vec->clone();
00478     if(out && print_all_tests())
00479       *out << "vec = NULL\n";
00480     vec = Teuchos::null;
00481     const value_type
00482       sum_vec       = sum(*v[0]),
00483       sum_vec_clone = sum(*vec_clone);
00484     err = (sum_vec - sum_vec_clone)/n;
00485     if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) )
00486       *out << "check: (sum(v[0]) - sum(vec_clone))/n = ("
00487          << sum_vec << " - " << sum_vec_clone << ")/" << n
00488          << " = " << err << std::endl;
00489     check_test(err,out,&success);
00490   }
00491 
00492   } // end try
00493   catch(const std::exception& except) {
00494     if(out)
00495       *out << "Caught a std::exception: " << except.what() << std::endl;
00496     success = false;
00497     if(throw_exception())
00498       throw;
00499   }
00500   catch(...) {
00501     if(out)
00502       *out << "Caught an unknown exception!\n";
00503     success = false;
00504     if(throw_exception())
00505       throw;
00506   }
00507 
00508   return success;
00509 }
00510 
00511 void VectorSpaceTester::check_test(value_type err, std::ostream* out, bool* success) const
00512 {
00513   if( ::fabs(err) >= error_tol() ) *success = false;
00514   if(out && (print_all_tests() || ::fabs(err) >= warning_tol()) ) {
00515     if( ::fabs(err) >= error_tol() )
00516       *out << "Error!  |" << err << "| = " << ::fabs(err) << " >= error_tol = "
00517          << error_tol() << std::endl;
00518     else if( ::fabs(err) >= warning_tol() )
00519       *out << "Warning!  |" << err << "| = " << ::fabs(err) << " >= warning_tol = "
00520              << warning_tol() << std::endl;
00521   }
00522   TEST_FOR_EXCEPTION(
00523     !*success && this->throw_exception(), std::logic_error
00524     ,"VectorSpaceTester::check_test(...): Error!  |" << err << "| = " << ::fabs(err) << " >= error_tol = "
00525     << error_tol() << std::endl );
00526 }
00527 
00528 } // end namespace AbstractLinAlgPack

Generated on Wed May 12 21:50:44 2010 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.4.7