AbstractLinAlgPack_BasisSystemTester.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 #include <limits>
00033 
00034 #include "AbstractLinAlgPack_BasisSystemTester.hpp"
00035 #include "AbstractLinAlgPack_BasisSystem.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 BasisSystemTester::BasisSystemTester(
00049   EPrintTestLevel  print_tests
00050   ,bool            dump_all
00051   ,bool            throw_exception
00052   ,size_type       num_random_tests
00053   ,value_type      warning_tol
00054   ,value_type      error_tol
00055   )
00056   :print_tests_(print_tests)
00057   ,dump_all_(dump_all)
00058   ,throw_exception_(throw_exception)
00059   ,num_random_tests_(num_random_tests)
00060   ,warning_tol_(warning_tol)
00061   ,error_tol_(error_tol)
00062 {}
00063 
00064 bool BasisSystemTester::test_basis_system(
00065   const BasisSystem           &bs
00066   ,const MatrixOp             *Gc
00067   ,const MatrixOpNonsing      *C
00068   ,const MatrixOp             *N_in
00069   ,const MatrixOp             *D
00070   ,const MatrixOp             *GcUP
00071   ,std::ostream               *out
00072   )
00073 {
00074   namespace rcp = MemMngPack;
00075   using BLAS_Cpp::no_trans;
00076   using BLAS_Cpp::trans;
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, llresult;
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   EPrintTestLevel
00098     print_tests = ( this->print_tests() == PRINT_NOT_SELECTED ? PRINT_NONE : this->print_tests() );
00099 
00100   MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
00101   
00102   // Print the input?
00103   if( out && print_tests != PRINT_NONE ) {
00104     if( print_tests >= PRINT_BASIC ) {
00105       *out << "\n*************************************************"
00106          << "\n*** BasisSystemTester::test_basis_system(...) ***"
00107          << "\n*************************************************\n";
00108       if(Gc)
00109         *out << "\n||Gc||inf   = " << Gc->calc_norm(mat_nrm_inf).value;
00110       if(C) {
00111         *out << "\n||C||inf    = " << C->calc_norm(mat_nrm_inf).value;
00112         *out << "\ncond_inf(C) = " << C->calc_cond_num(mat_nrm_inf).value;
00113       }
00114       if(N_in)
00115         *out << "\n||N||inf    = " << N_in->calc_norm(mat_nrm_inf).value;
00116       if(D)
00117         *out << "\n||D||inf    = " << D->calc_norm(mat_nrm_inf).value;
00118       if(GcUP)
00119         *out << "\n||GcUP||inf = " << GcUP->calc_norm(mat_nrm_inf).value;
00120     }
00121     if(dump_all()) {
00122       if(Gc)
00123         *out << "\nGc =\n"    << *Gc;
00124       if(C)
00125         *out << "\nC =\n"     << *C;
00126       if(N_in)
00127         *out << "\nN =\n"     << *N_in;
00128       if(D)
00129         *out << "\nD =\n"     << *D;
00130       if(GcUP)
00131         *out << "\nGcUP =\n"  << *GcUP;
00132     }
00133   }
00134 
00135   //
00136   // Check the dimensions of everything
00137   //
00138 
00139   const Range1D
00140     var_dep          = bs.var_dep(),
00141     var_indep        = bs.var_indep(),
00142     equ_decomp       = bs.equ_decomp(),
00143     equ_undecomp     = bs.equ_undecomp();
00144 
00145   if( out && print_tests >= PRINT_MORE ) {
00146     *out
00147       << "\nbs.var_dep()        = ["<<var_dep.lbound()<<","<<var_dep.ubound()<<"]"
00148       << "\nbs.var_indep( )     = ["<<var_indep.lbound()<<","<<var_indep.ubound()<<"]"
00149       << "\nbs.equ_decomp()     = ["<<equ_decomp.lbound()<<","<<equ_decomp.ubound()<<"]"
00150       << "\nbs.equ_undecomp()   = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
00151       << std::endl;
00152   }
00153 
00154   if( out && print_tests >= PRINT_BASIC )
00155     *out << "\n1) Check the partitioning ranges ...";
00156   lresult = true;
00157 
00158   if( out && print_tests >= PRINT_MORE )
00159     *out << "\n\n1.a) check: var_dep.size() != equ_decomp.size() : ";
00160   result = var_dep.size() == equ_decomp.size();
00161   if(out && print_tests >= PRINT_MORE)
00162     *out << ( result ? "passed" : "failed" );
00163   if(!result) lresult = false;
00164 
00165   if(Gc) {
00166     if( out && print_tests >= PRINT_MORE )
00167       *out << "\n1.b) check: var_dep.size() + var_indep.size() == Gc->rows() : ";
00168     result = var_dep.size() + var_indep.size() == Gc->rows();
00169     if(out && print_tests >= PRINT_MORE)
00170       *out << ( result ? "passed" : "failed" );
00171     if(!result) lresult = false;
00172   } 
00173 
00174   if(Gc) {
00175     if( out && print_tests >= PRINT_MORE )
00176       *out << "\n1.d) check: equ_decomp.size() + equ_undecomp.size() == Gc->cols() : ";
00177     result = equ_decomp.size() + equ_undecomp.size() == Gc->cols();
00178     if(out && print_tests >= PRINT_MORE)
00179       *out << ( result ? "passed" : "failed" );
00180     if(!result) lresult = false;
00181   } 
00182 
00183   if(out && print_tests >= PRINT_MORE)
00184     *out << std::endl;
00185 
00186   if(!lresult) success = false;
00187   if( out && print_tests == PRINT_BASIC )
00188     *out << " : " << ( lresult ? "passed" : "failed" );
00189 
00190   // Create the N matrix if not input
00191   Teuchos::RCP<const AbstractLinAlgPack::MatrixOp>
00192     N = Teuchos::rcp(N_in,false);
00193   if( Gc && C && N.get() == NULL ) {
00194     if(out && print_tests >= PRINT_BASIC)
00195       *out
00196         << "\nCreating the matrix N since it was not input by the client ...";
00197     if(out && print_tests >= PRINT_MORE)
00198       *out
00199         << std::endl;
00200     Teuchos::RCP<AbstractLinAlgPack::MatrixComposite>
00201       N_comp = Teuchos::rcp(new AbstractLinAlgPack::MatrixComposite(var_dep.size(),var_indep.size()));
00202     if( equ_decomp.size() )
00203       N_comp->add_matrix( 0, 0, 1.0, equ_decomp, Gc, Teuchos::null, BLAS_Cpp::trans, var_indep );
00204     N_comp->finish_construction(
00205       Gc->space_rows().sub_space(equ_decomp)->clone()
00206       ,Gc->space_cols().sub_space(var_indep)->clone()
00207       );
00208     if( out && dump_all() )
00209       *out << "\nN =\n" << *N_comp;
00210     N = N_comp;
00211   }
00212 
00213   // Create the other auxillary matrix objects
00214   if( equ_undecomp.size() ) {
00215     TEST_FOR_EXCEPT(true); // ToDo: Create matrix objects for Gc(var_dep,equ_undecomp) and Gc(var_indep,equ_undecomp)
00216   }
00217 
00218   //
00219   // Perform the tests
00220   //
00221 
00222   if( C && N.get() ) {
00223 
00224     if(out && print_tests >= PRINT_BASIC)
00225       *out
00226         << "\n2) Check the compatibility of the vector spaces for C, N, D and Gc ...";
00227     lresult = true;
00228 
00229     if(out && print_tests >= PRINT_MORE)
00230       *out
00231         << "\n2.a) Check consistency of the vector spaces for:"
00232         << "\n    C.space_cols() == N.space_cols()";
00233     llresult = true;
00234     if(out && print_tests >= PRINT_ALL)
00235        *out << "\n\n2.a.1) C->space_cols().is_compatible(N->space_cols()) == true : ";
00236     result = C->space_cols().is_compatible(N->space_cols());  
00237     if(out && print_tests >= PRINT_ALL)
00238        *out << ( result ? "passed" : "failed" )
00239          << std::endl;
00240     if(!result) llresult = false;
00241     if(!llresult) lresult = false;
00242     if( out && print_tests == PRINT_MORE )
00243     *out << " : " << ( llresult ? "passed" : "failed" );
00244 
00245     if(D) {
00246       if(out && print_tests >= PRINT_MORE)
00247         *out
00248           << "\n2.b) Check consistency of the vector spaces for:"
00249           << "\n    D.space_cols() == C.space_cols() and D.space_rows() == N.space_rows()";
00250       llresult = true;
00251       if(out && print_tests >= PRINT_ALL)
00252         *out << "\n2.b.1) D->space_cols().is_compatible(C->space_cols()) == true : ";
00253       result = D->space_cols().is_compatible(C->space_cols());
00254       if(out && print_tests >= PRINT_ALL)
00255         *out << ( result ? "passed" : "failed" );
00256       if(!result) llresult = false;
00257       if(out && print_tests >= PRINT_ALL)
00258         *out << "\n2.b.2) D->space_rows().is_compatible(N->space_rows()) == true : ";
00259       result = D->space_rows().is_compatible(N->space_rows());
00260       if(out && print_tests >= PRINT_ALL)
00261         *out << ( result ? "passed" : "failed" )
00262            << std::endl;
00263       if(!result) llresult = false;
00264       if(!llresult) lresult = false;
00265       if( out && print_tests == PRINT_MORE )
00266         *out << " : " << ( llresult ? "passed" : "failed" );
00267     }
00268 
00269     if(out && print_tests >= PRINT_MORE)
00270       *out
00271         << "\n2.c) Check consistency of the vector spaces for:"
00272         << "\n    Gc'(equ_decomp,  var_dep) == C";
00273     llresult = true;
00274     if( equ_decomp.size() ) {
00275       if(out && print_tests >= PRINT_ALL)
00276         *out << "\n2.c.1) Gc->space_rows().sub_space(equ_decomp)->is_compatible(*C->space_cols().sub_space(equ_decomp)) == true : ";
00277       result = Gc->space_rows().sub_space(equ_decomp)->is_compatible(*C->space_cols().sub_space(equ_decomp));
00278       if(out && print_tests >= PRINT_ALL)
00279         *out << ( result ? "passed" : "failed" );
00280       if(!result) llresult = false;
00281       if(out && print_tests >= PRINT_ALL)
00282         *out << "\n2.c.2) Gc->space_cols().sub_space(var_dep)->is_compatible(C->space_rows()) == true : ";
00283       result = Gc->space_cols().sub_space(var_dep)->is_compatible(C->space_rows());
00284       if(out && print_tests >= PRINT_ALL)
00285         *out << ( result ? "passed" : "failed" );
00286       if(!result) llresult = false;
00287     }
00288     if(out && print_tests >= PRINT_ALL)
00289       *out << std::endl;
00290     if(!llresult) lresult = false;
00291     if( out && print_tests == PRINT_MORE )
00292       *out << " : " << ( llresult ? "passed" : "failed" );
00293 
00294     if(out && print_tests >= PRINT_MORE)
00295       *out
00296         << "\n2.d) Check consistency of the vector spaces for:"
00297         << "\n    Gc'(equ_decomp, var_indep) == N";
00298     llresult = true;
00299     if( equ_decomp.size() ) {
00300       if(out && print_tests >= PRINT_ALL)
00301         *out << "\n2.d.1) Gc->space_rows().sub_space(equ_decomp)->is_compatible(*N->space_cols().sub_space(equ_decomp)) == true : ";
00302       result = Gc->space_rows().sub_space(equ_decomp)->is_compatible(*N->space_cols().sub_space(equ_decomp));
00303       if(out && print_tests >= PRINT_ALL)
00304         *out << ( result ? "passed" : "failed" );
00305       if(!result) llresult = false;
00306       if(out && print_tests >= PRINT_ALL)
00307         *out << "\n2.d.2) Gc->space_cols().sub_space(var_indep)->is_compatible(N->space_rows()) == true : ";
00308       result = Gc->space_cols().sub_space(var_indep)->is_compatible(N->space_rows());
00309       if(out && print_tests >= PRINT_ALL)
00310         *out << ( result ? "passed" : "failed" );
00311       if(!result) llresult = false;
00312     }
00313     if(out && print_tests >= PRINT_ALL)
00314       *out << std::endl;
00315     if(!llresult) lresult = false;
00316     if( out && print_tests == PRINT_MORE )
00317       *out << " : " << ( llresult ? "passed" : "failed" )
00318          << std::endl;
00319 
00320     if(!lresult) success = false;
00321     if( out && print_tests == PRINT_BASIC )
00322       *out << " : " << ( lresult ? "passed" : "failed" );
00323 
00324     if(out && print_tests >= PRINT_BASIC)
00325       *out
00326         << "\n3) Check the compatibility of the matrices C, N, D and Gc numerically ...";
00327 
00328     if(out && print_tests >= PRINT_MORE)
00329       *out
00330         << std::endl
00331         << "\n3.a) Check consistency of:"
00332         << "\n                "
00333         << "\n    op ( alpha* [ Gc'(equ_decomp,  var_dep) Gc'(equ_decomp,  var_indep) ] ) * v"
00334         << "\n         \\______________________________________________________________/"
00335         << "\n                                        A"
00336         << "\n    ==  op( alpha*[ C  N ] ) * v"
00337         << "\n            \\____________/"
00338         << "\n                   B"
00339         << "\nfor random vectors v ...";
00340     {
00341 
00342       VectorSpace::vec_mut_ptr_t
00343         Gc_v_x    = Gc->space_cols().create_member(),
00344         Gc_v_c    = Gc->space_rows().create_member(),
00345         C_v_xD    = C->space_rows().create_member(),
00346         C_v_chD   = C->space_cols().create_member(),
00347         N_v_xI    = N->space_rows().create_member(),
00348         N_v_chD   = N->space_cols().create_member(),
00349         v_x       = Gc->space_cols().create_member(),
00350         v_x_tmp   = v_x->space().create_member(),
00351         v_chD     = C_v_xD->space().create_member(),
00352         v_chD_tmp = v_chD->space().create_member();
00353 
00354       if(out && print_tests >= PRINT_MORE)
00355         *out << "\n\n3.a.1) Testing non-transposed A*v == B*v ...";
00356       if(out && print_tests > PRINT_MORE)
00357         *out << std::endl;
00358       llresult = true;
00359        {for( int k = 1; k <= num_random_tests(); ++k ) {
00360         random_vector( rand_y_l, rand_y_u, v_x.get() );
00361          if(out && print_tests >= PRINT_ALL) {
00362           *out
00363             << "\n3.a.1."<<k<<") random vector " << k << " ( ||v_x||_1 / n = " << (v_x->norm_1() / v_x->dim()) << " )\n";
00364           if(dump_all() && print_tests >= PRINT_ALL)
00365             *out << "\nv_x =\n" << *v_x;
00366         }
00367         if(Gc && equ_decomp.size()) {
00368           V_StMtV( Gc_v_c.get(), alpha, *Gc, trans, *v_x );
00369           *v_chD_tmp->sub_view(equ_decomp)
00370             = *Gc_v_c->sub_view(equ_decomp);
00371         }
00372         V_StMtV( C_v_chD.get(), alpha, *C, no_trans, *v_x->sub_view(var_dep) );
00373         V_StMtV( N_v_chD.get(), alpha, *N, no_trans, *v_x->sub_view(var_indep) );
00374         V_VpV( v_chD.get(), *C_v_chD, *N_v_chD );
00375         const value_type
00376           sum_Bv  = sum(*v_chD),
00377           sum_Av  = sum(*v_chD_tmp);
00378         assert_print_nan_inf(sum_Bv, "sum(B*v_x)",true,out);
00379         assert_print_nan_inf(sum_Av, "sum(A*v_x)",true,out);
00380         const value_type
00381           calc_err = ::fabs( ( sum_Av - sum_Bv )
00382                      /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
00383         if(out && print_tests >= PRINT_ALL)
00384           *out
00385             << "\nrel_err(sum(A*v_x),sum(B*v_x)) = "
00386             << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00387             << calc_err << std::endl;
00388         if( calc_err >= warning_tol() ) {
00389           if(out && print_tests >= PRINT_ALL)
00390             *out
00391               << std::endl
00392               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00393               << ", rel_err(sum(A*v_x),sum(B*v_x)) = "
00394               << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00395               << calc_err
00396               << " exceeded "
00397               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00398               << " = "
00399               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00400               << std::endl;
00401           if(calc_err >= error_tol()) {
00402             if(dump_all() && print_tests >= PRINT_ALL) {
00403               *out << "\nalpha = "             << alpha << std::endl;
00404               *out << "\nv_x =\n"              << *v_x;
00405               *out << "\nalpha*Gc*v_x =\n"     << *Gc_v_c;
00406               *out << "A*v =\n"                << *v_chD_tmp;
00407               *out << "\nalpha*C*v_x =\n"      << *C_v_chD;
00408               *out << "\nalpha*N*v_x =\n"      << *N_v_chD;
00409               *out << "\nB*v =\n"              << *v_chD;
00410             }
00411             llresult = false;
00412           }
00413         }
00414       }}
00415       if(!llresult) lresult = false;
00416       if( out && print_tests == PRINT_MORE )
00417         *out << " : " << ( llresult ? "passed" : "failed" )
00418            << std::endl;
00419 
00420       if(out && print_tests >= PRINT_MORE)
00421         *out << "\n3.a.2) Testing transposed A'*v == B'*v ...";
00422       if(out && print_tests > PRINT_MORE)
00423         *out << std::endl;
00424       llresult = true;
00425       {for( int k = 1; k <= num_random_tests(); ++k ) {
00426         random_vector( rand_y_l, rand_y_u, v_chD.get() );
00427          if(out && print_tests >= PRINT_ALL) {
00428           *out
00429             << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_chD||_1 / n = " << (v_chD->norm_1() / v_chD->dim()) << " )\n";
00430           if(dump_all() && print_tests >= PRINT_ALL)
00431             *out << "\nv_chD =\n" << *v_chD;
00432         }
00433         *v_x_tmp = 0.0;
00434         if(Gc && equ_decomp.size()) {
00435           *Gc_v_c->sub_view(equ_decomp) = *v_chD->sub_view(equ_decomp);
00436           if(equ_undecomp.size())
00437             *Gc_v_c->sub_view(equ_undecomp) = 0.0;
00438           V_StMtV( Gc_v_x.get(), alpha, *Gc, no_trans, *Gc_v_c );
00439           Vp_V( v_x_tmp.get(), *Gc_v_x );
00440         }
00441         V_StMtV( C_v_xD.get(), alpha, *C, trans, *v_chD );
00442         *v_x->sub_view(var_dep) = *C_v_xD;
00443         V_StMtV( N_v_xI.get(), alpha, *N, trans, *v_chD );
00444         *v_x->sub_view(var_indep) = *N_v_xI;
00445         const value_type
00446           sum_BTv  = sum(*v_x),
00447           sum_ATv  = sum(*v_x_tmp);
00448         assert_print_nan_inf(sum_BTv, "sum(B'*v_chD)",true,out);
00449         assert_print_nan_inf(sum_ATv, "sum(A'*v_chD)",true,out);
00450         const value_type
00451           calc_err = ::fabs( ( sum_ATv - sum_BTv )
00452                      /( ::fabs(sum_ATv) + ::fabs(sum_BTv) + small_num ) );
00453         if(out && print_tests >= PRINT_ALL)
00454           *out
00455             << "\nrel_err(sum(A'*v_chD),sum(B'*v_chD)) = "
00456             << "rel_err(" << sum_ATv << "," << sum_BTv << ") = "
00457             << calc_err << std::endl;
00458         if( calc_err >= warning_tol() ) {
00459           if(out && print_tests >= PRINT_ALL)
00460             *out
00461               << std::endl
00462               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00463               << ", rel_err(sum(A'*v_chD),sum(B'*v_chD)) = "
00464               << "rel_err(" << sum_ATv << "," << sum_BTv << ") = "
00465               << calc_err << std::endl
00466               << " exceeded "
00467               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00468               << " = "
00469               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00470               << std::endl;
00471           if(calc_err >= error_tol()) {
00472             if(dump_all() && print_tests >= PRINT_ALL) {
00473               *out << "\nalpha = " << alpha << std::endl;
00474               *out << "\nv_chD =\n"            << *v_chD;
00475               if(Gc_v_x.get() && equ_decomp.size()) {
00476                 *out << "\nGc_v_c =\n"       << *Gc_v_c;
00477                 *out << "\nalpha*Gc'*[v_chD(equ_decomp); 0] =\n"
00478                    << *Gc_v_x;
00479               }
00480               *out << "A'*v =\n"               << *v_x_tmp;
00481               *out << "\nalpha*C*v_chD =\n"    << *C_v_xD;
00482               *out << "\nalpha*N*v_chD =\n"    << *N_v_xI;
00483               *out << "\nB'*v =\n"             << *v_x;
00484             }
00485             llresult = false;
00486           }
00487         }
00488       }}
00489       if(!llresult) lresult = false;
00490       if( out && print_tests == PRINT_MORE )
00491         *out << " : " << ( llresult ? "passed" : "failed" )
00492            << std::endl;
00493     }
00494     
00495     if(out && print_tests >= PRINT_MORE)
00496       *out
00497         << "\n3.b) Check consistency of:"
00498         << "\n    alpha*op(C)*(op(inv(C)) * v) == alpha*v"
00499         << "\nfor random vectors v ...";
00500     {
00501       VectorSpace::vec_mut_ptr_t
00502         v_xD      = C->space_rows().create_member(),
00503         v_xD_tmp  = C->space_rows().create_member(),
00504         v_chD     = C->space_cols().create_member(),
00505         v_chD_tmp = C->space_cols().create_member();
00506 
00507       if(out && print_tests >= PRINT_MORE)
00508         *out << "\n\n3.b.1) Testing non-transposed: alpha*C*(inv(C)*v) == alpha*v ...";
00509       if(out && print_tests > PRINT_MORE)
00510         *out << std::endl;
00511       llresult = true;
00512        {for( int k = 1; k <= num_random_tests(); ++k ) {
00513         random_vector( rand_y_l, rand_y_u, v_chD.get() );
00514          if(out && print_tests >= PRINT_ALL) {
00515           *out
00516             << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_chD||_1 / n = " << (v_chD->norm_1() / v_chD->dim()) << " )\n";
00517           if(dump_all() && print_tests >= PRINT_ALL)
00518             *out << "\nv_chD =\n" << *v_chD;
00519         }
00520         V_InvMtV( v_xD_tmp.get(), *C, no_trans, *v_chD );
00521         V_StMtV( v_chD_tmp.get(), alpha, *C, no_trans, *v_xD_tmp );
00522         const value_type
00523           sum_aCICv  =         sum(*v_chD_tmp),
00524           sum_av     = alpha * sum(*v_chD);
00525         assert_print_nan_inf(sum_aCICv, "sum(alpha*C*(inv(C)*v)",true,out);
00526         assert_print_nan_inf(sum_av, "sum(alpha*v)",true,out);
00527         const value_type
00528           calc_err = ::fabs( ( sum_aCICv - sum_av )
00529                      /( ::fabs(sum_aCICv) + ::fabs(sum_av) + small_num ) );
00530         if(out && print_tests >= PRINT_ALL)
00531           *out
00532             << "\nrel_err(sum(alpha*C*(inv(C)*v),sum(alpha*v)) = "
00533             << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
00534             << calc_err << std::endl;
00535         if( calc_err >= warning_tol() ) {
00536           if(out && print_tests >= PRINT_ALL)
00537             *out
00538               << std::endl
00539               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00540               << ", rel_err(sum(alpha*C*(inv(C)*v)),sum(alpha*v)) = "
00541               << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
00542               << calc_err
00543               << " exceeded "
00544               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00545               << " = "
00546               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00547               << std::endl;
00548           if(calc_err >= error_tol()) {
00549             if(dump_all() && print_tests >= PRINT_ALL) {
00550               *out << "\nalpha = " << alpha << std::endl;
00551               *out << "\nv_chD =\n"                << *v_chD;
00552               *out << "\ninv(C)*v_chD =\n"         << *v_xD_tmp;
00553               *out << "\nalpha*C*inv(C)*v_chD =\n" << *v_chD_tmp;
00554             }
00555             llresult = false;
00556           }
00557         }
00558       }}
00559       if(!llresult) lresult = false;
00560       if( out && print_tests == PRINT_MORE )
00561         *out << " : " << ( llresult ? "passed" : "failed" )
00562            << std::endl;
00563 
00564       if(out && print_tests >= PRINT_MORE)
00565         *out << "\n3.b.2) Testing transposed: alpha*C'*(inv(C')*v) == alpha*v ...";
00566       if(out && print_tests > PRINT_MORE)
00567         *out << std::endl;
00568       llresult = true;
00569        {for( int k = 1; k <= num_random_tests(); ++k ) {
00570         random_vector( rand_y_l, rand_y_u, v_xD.get() );
00571          if(out && print_tests >= PRINT_ALL) {
00572           *out
00573             << "\n3.b.2."<<k<<") random vector " << k << " ( ||v_xD||_1 / n = " << (v_xD->norm_1() / v_xD->dim()) << " )\n";
00574           if(dump_all() && print_tests >= PRINT_ALL)
00575             *out << "\nv_xD =\n" << *v_xD;
00576         }
00577         V_InvMtV( v_chD_tmp.get(), *C, trans, *v_xD );
00578         V_StMtV( v_xD_tmp.get(), alpha, *C, trans, *v_chD_tmp );
00579         const value_type
00580           sum_aCICv  =         sum(*v_xD_tmp),
00581           sum_av     = alpha * sum(*v_xD);
00582         assert_print_nan_inf(sum_aCICv, "sum(alpha*C'*(inv(C')*v)",true,out);
00583         assert_print_nan_inf(sum_av, "sum(alpha*v)",true,out);
00584         const value_type
00585           calc_err = ::fabs( ( sum_aCICv - sum_av )
00586                      /( ::fabs(sum_aCICv) + ::fabs(sum_av) + small_num ) );
00587         if(out && print_tests >= PRINT_ALL)
00588           *out
00589             << "\nrel_err(sum(alpha*C'*(inv(C')*v)),sum(alpha*v)) = "
00590             << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
00591             << calc_err << std::endl;
00592         if( calc_err >= warning_tol() ) {
00593           if(out && print_tests >= PRINT_ALL)
00594             *out
00595               << std::endl
00596               << ( calc_err >= error_tol() ? "Error" : "Warning" )
00597               << ", rel_err(sum(alpha*C'*(inv(C')*v)),sum(alpha*v)) = "
00598               << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
00599               << calc_err
00600               << " exceeded "
00601               << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00602               << " = "
00603               << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00604               << std::endl;
00605           if(calc_err >= error_tol()) {
00606             if(dump_all() && print_tests >= PRINT_ALL) {
00607               *out << "\nalpha = " << alpha << std::endl;
00608               *out << "\nv_xD =\n"                   << *v_xD;
00609               *out << "\ninv(C')*v_xD =\n"           << *v_chD_tmp;
00610               *out << "\nalpha*C'*inv(C')*v_xD =\n"  << *v_xD_tmp;
00611             }
00612             llresult = false;
00613           }
00614         }
00615       }}
00616       if(!llresult) lresult = false;
00617       if( out && print_tests == PRINT_MORE )
00618         *out << " : " << ( llresult ? "passed" : "failed" )
00619            << std::endl;
00620     }
00621     
00622     if(D) {
00623       if(out && print_tests >= PRINT_MORE)
00624         *out
00625           << "\n3.c) Check consistency of:"
00626           << "\n    alpha * op(-inv(C) * N) * v == alpha * op(D) * v"
00627           << "\nfor random vectors v ...";
00628       
00629     {
00630         VectorSpace::vec_mut_ptr_t
00631           v_xD      = C->space_rows().create_member(),
00632           v_xI      = N->space_rows().create_member(),
00633           v_xD_tmp  = C->space_rows().create_member(),
00634           v_xI_tmp  = N->space_rows().create_member(),
00635           v_chD_tmp = C->space_cols().create_member();
00636 
00637         if(out && print_tests >= PRINT_MORE)
00638           *out << "\n\n3.b.1) Testing non-transposed: inv(C)*(-alpha*N*v) == alpha*D*v ...";
00639         if(out && print_tests > PRINT_MORE)
00640           *out << std::endl;
00641         llresult = true;
00642          {for( int k = 1; k <= num_random_tests(); ++k ) {
00643           random_vector( rand_y_l, rand_y_u, v_xI.get() );
00644            if(out && print_tests >= PRINT_ALL) {
00645             *out
00646               << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_xI||_1 / n = " << (v_xI->norm_1() / v_xI->dim()) << " )\n";
00647             if(dump_all() && print_tests >= PRINT_ALL)
00648               *out << "\nv_xI =\n" << *v_xI;
00649           }
00650           V_StMtV( v_chD_tmp.get(), -alpha, *N, no_trans, *v_xI );
00651           V_InvMtV( v_xD_tmp.get(), *C, no_trans, *v_chD_tmp );
00652           V_StMtV( v_xD.get(), alpha, *D, no_trans, *v_xI );
00653           const value_type
00654             sum_ICaNv  = sum(*v_xD_tmp),
00655             sum_aDv    = sum(*v_xD);
00656           assert_print_nan_inf(sum_ICaNv, "sum(inv(C)*(-alpha*N*v))",true,out);
00657           assert_print_nan_inf(sum_aDv, "sum(alpha*D*v)",true,out);
00658           const value_type
00659             calc_err = ::fabs( ( sum_ICaNv - sum_aDv )
00660                        /( ::fabs(sum_ICaNv) + ::fabs(sum_aDv) + small_num ) );
00661           if(out && print_tests >= PRINT_ALL)
00662             *out
00663               << "\nrel_err(sum(inv(C)*(-alpha*N*v)),sum(alpha*D*v)) = "
00664               << "rel_err(" << sum_ICaNv << "," << sum_aDv << ") = "
00665               << calc_err << std::endl;
00666           if( calc_err >= warning_tol() ) {
00667             if(out && print_tests >= PRINT_ALL)
00668               *out
00669                 << std::endl
00670                 << ( calc_err >= error_tol() ? "Error" : "Warning" )
00671                 << ", rel_err(sum(inv(C)*(-alpha*N*v))),sum(alpha*D*v)) = "
00672                 << "rel_err(" << sum_ICaNv << "," << sum_aDv << ") = "
00673                 << calc_err
00674                 << " exceeded "
00675                 << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00676                 << " = "
00677                 << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00678                 << std::endl;
00679             if(calc_err >= error_tol()) {
00680               if(dump_all() && print_tests >= PRINT_ALL) {
00681                 *out << "\nalpha = " << alpha << std::endl;
00682                 *out << "\nv_xI =\n"                   << *v_xI;
00683                 *out << "\n-alpha*N*v_xI =\n"          << *v_chD_tmp;
00684                 *out << "\ninv(C)*(-alpha*N*v_xI) =\n" << *v_xD_tmp;
00685                 *out << "\nalpha*D*v_xI =\n"           << *v_xD;
00686               }
00687               llresult = false;
00688             }
00689           }
00690         }}
00691         if(!llresult) lresult = false;
00692         if( out && print_tests == PRINT_MORE )
00693           *out << " : " << ( llresult ? "passed" : "failed" )
00694              << std::endl;
00695 
00696         if(out && print_tests >= PRINT_MORE)
00697           *out << "\n3.b.1) Testing transposed: -alpha*N'*(inv(C')*v) == alpha*D'*v ...";
00698         if(out && print_tests > PRINT_MORE)
00699           *out << std::endl;
00700         llresult = true;
00701          {for( int k = 1; k <= num_random_tests(); ++k ) {
00702           random_vector( rand_y_l, rand_y_u, v_xD.get() );
00703            if(out && print_tests >= PRINT_ALL) {
00704             *out
00705               << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_xD||_1 / n = " << (v_xD->norm_1() / v_xD->dim()) << " )\n";
00706             if(dump_all() && print_tests >= PRINT_ALL)
00707               *out << "\nv_xD =\n" << *v_xD;
00708           }
00709           V_InvMtV( v_chD_tmp.get(), *C, trans, *v_xD );
00710           V_StMtV( v_xI_tmp.get(), -alpha, *N, trans, *v_chD_tmp );
00711           V_StMtV( v_xI.get(), alpha, *D, trans, *v_xD );
00712           const value_type
00713             sum_aNTICTv  = sum(*v_xI_tmp),
00714             sum_aDTv     = sum(*v_xI);
00715           assert_print_nan_inf(sum_aNTICTv, "sum(-alpha*N'*(inv(C')*v))",true,out);
00716           assert_print_nan_inf(sum_aDTv, "sum(alpha*D'*v)",true,out);
00717           const value_type
00718             calc_err = ::fabs( ( sum_aNTICTv - sum_aDTv )
00719                        /( ::fabs(sum_aNTICTv) + ::fabs(sum_aDTv) + small_num ) );
00720           if(out && print_tests >= PRINT_ALL)
00721             *out
00722               << "\nrel_err(sum(-alpha*N'*(inv(C')*v)),sum(alpha*D'*v)) = "
00723               << "rel_err(" << sum_aNTICTv << "," << sum_aDTv << ") = "
00724               << calc_err << std::endl;
00725           if( calc_err >= warning_tol() ) {
00726             if(out && print_tests >= PRINT_ALL)
00727               *out
00728                 << std::endl
00729                 << ( calc_err >= error_tol() ? "Error" : "Warning" )
00730                 << ", rel_err(sum(-alpha*N'*(inv(C')*v))),sum(alpha*D'*v)) = "
00731                 << "rel_err(" << sum_aNTICTv << "," << sum_aDTv << ") = "
00732                 << calc_err
00733                 << " exceeded "
00734                 << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
00735                 << " = "
00736                 << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
00737                 << std::endl;
00738             if(calc_err >= error_tol()) {
00739               if(dump_all() && print_tests >= PRINT_ALL) {
00740                 *out << "\nalpha = " << alpha << std::endl;
00741                 *out << "\nv_xD =\n"                      << *v_xD;
00742                 *out << "\ninv(C')**v_xD =\n"             << *v_chD_tmp;
00743                 *out << "\n-alpha*N'*(inv(C')**v_xD) =\n" << *v_xI_tmp;
00744                 *out << "\nalpha*D*'v_xD =\n"             << *v_xI;
00745               }
00746               llresult = false;
00747             }
00748           }
00749         }}
00750         if(!llresult) lresult = false;
00751         if( out && print_tests == PRINT_MORE )
00752           *out << " : " << ( llresult ? "passed" : "failed" )
00753              << std::endl;
00754 
00755       }
00756     }
00757     
00758     if( GcUP ) {
00759         TEST_FOR_EXCEPT(true); // ToDo: Validate GcUP and the related matrices
00760     }
00761 
00762     if(!lresult) success = false;
00763     if( out && print_tests == PRINT_BASIC )
00764       *out << " : " << ( lresult ? "passed" : "failed" );
00765   }
00766 
00767   if(out && print_tests != PRINT_NONE ) {
00768     if(success)
00769       *out << "\nCongradulations! The BasisSystem object and its associated matrix objects seem to check out!\n";
00770     else
00771       *out << "\nOops! At last one of the tests did not check out!\n";
00772     if(print_tests >= PRINT_BASIC)
00773       *out << "\nEnd BasisSystemTester::test_basis_system(...)\n";
00774   }
00775   
00776   return success;
00777 }
00778 
00779 } // end namespace AbstractLinAlgPack

Generated on Tue Jul 13 09:28:50 2010 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.4.7