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