ConstrainedOptPack_DecompositionSystemTester.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <math.h>
00030 
00031 #include <limits>
00032 #include <ostream>
00033 
00034 #include "ConstrainedOptPack_DecompositionSystemTester.hpp"
00035 #include "ConstrainedOptPack_DecompositionSystem.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_MatrixOpNonsingTester.hpp"
00042 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00043 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00044 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00045 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00046 
00047 namespace ConstrainedOptPack {
00048 
00049 DecompositionSystemTester::DecompositionSystemTester(
00050   EPrintTestLevel  print_tests
00051   ,bool            dump_all
00052   ,bool            throw_exception
00053   ,size_type       num_random_tests
00054   ,value_type      mult_warning_tol
00055   ,value_type      mult_error_tol
00056   ,value_type      solve_warning_tol
00057   ,value_type      solve_error_tol
00058   )
00059   :print_tests_(print_tests)
00060   ,dump_all_(dump_all)
00061   ,throw_exception_(throw_exception)
00062   ,num_random_tests_(num_random_tests)
00063   ,mult_warning_tol_(mult_warning_tol)
00064   ,mult_error_tol_(mult_error_tol)
00065   ,solve_warning_tol_(solve_warning_tol)
00066   ,solve_error_tol_(solve_error_tol)
00067 {}
00068  
00069 bool DecompositionSystemTester::test_decomp_system(
00070   const DecompositionSystem   &ds
00071   ,const MatrixOp             &Gc
00072   ,const MatrixOp             *Z
00073   ,const MatrixOp             *Y
00074   ,const MatrixOpNonsing      *R
00075   ,const MatrixOp             *Uz
00076   ,const MatrixOp             *Uy
00077   ,std::ostream               *out
00078   )
00079 {
00080   namespace rcp = MemMngPack;
00081   using BLAS_Cpp::no_trans;
00082   using BLAS_Cpp::trans;
00083   using AbstractLinAlgPack::sum;
00084   using AbstractLinAlgPack::dot;
00085   using AbstractLinAlgPack::Vp_StV;
00086   using AbstractLinAlgPack::Vp_StMtV;
00087   using AbstractLinAlgPack::assert_print_nan_inf;
00088   using AbstractLinAlgPack::random_vector;
00089   using LinAlgOpPack::V_StMtV;
00090   using LinAlgOpPack::V_MtV;
00091   using LinAlgOpPack::V_StV;
00092   using LinAlgOpPack::V_VpV;
00093   using LinAlgOpPack::Vp_V;
00094 
00095   bool success = true, result, lresult, llresult;
00096   const value_type
00097     rand_y_l  = -1.0,
00098     rand_y_u  = 1.0,
00099     small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
00100     alpha     = 2.0,
00101     beta      = 3.0;
00102 
00103   EPrintTestLevel
00104     print_tests = ( this->print_tests() == PRINT_NOT_SELECTED ? PRINT_NONE : this->print_tests() );
00105 
00106   // Print the input?
00107   if( out && print_tests != PRINT_NONE ) {
00108     if( print_tests >= PRINT_BASIC )
00109       *out << "\n**********************************************************"
00110          << "\n*** DecompositionSystemTester::test_decomp_system(...) ***"
00111          << "\n**********************************************************\n";
00112   }
00113 
00114   const size_type
00115     n = ds.n(),
00116     m = ds.m(),
00117     r = ds.r();
00118   const Range1D
00119     equ_decomp       = ds.equ_decomp(),
00120     equ_undecomp     = ds.equ_undecomp();
00121 
00122   // print dimensions, ranges
00123   if( out && print_tests >= PRINT_MORE ) {
00124     *out
00125       << "\nds.n()                  = " << n
00126       << "\nds.m()                  = " << m
00127       << "\nds.r()                  = " << r
00128       << "\nds.equ_decomp()         = ["<<equ_decomp.lbound()<<","<<equ_decomp.ubound()<<"]"
00129       << "\nds.equ_undecomp()       = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
00130       << "\nds.space_range()->dim() = " << ds.space_range()->dim()
00131       << "\nds.space_null()->dim()  = " << ds.space_null()->dim()
00132       << std::endl;
00133   }
00134 
00135   // validate input matrices
00136   TEST_FOR_EXCEPTION(
00137       Z==NULL&&Y==NULL&&R==NULL&&Uz==NULL&&Uy==NULL
00138     , std::invalid_argument
00139     ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
00140     "at least one of Z, Y, R, Uz or Uy can not be NULL!" );
00141   TEST_FOR_EXCEPTION(
00142     m == r && Uz != NULL, std::invalid_argument
00143     ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
00144     "Uz must be NULL if m==r is NULL!" );
00145   TEST_FOR_EXCEPTION(
00146     m == r && Uy != NULL, std::invalid_argument
00147     ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
00148     "Uy must be NULL if m==r is NULL!" );
00149 
00150   // Print the input?
00151   if( out && print_tests != PRINT_NONE ) {
00152     if(dump_all()) {
00153       *out << "\nGc =\n"       << Gc;
00154       if(Z)
00155         *out << "\nZ =\n"    << *Z;
00156       if(Y)
00157         *out << "\nY =\n"    << *Y;
00158       if(R)
00159         *out << "\nR =\n"    << *R;
00160       if(Uz)
00161         *out << "\nUz =\n"   << *Uz;
00162       if(Uy)
00163         *out << "\nUy =\n"   << *Uy;
00164     }
00165   }
00166 
00167   //
00168   // Check the dimensions of everything
00169   //
00170 
00171   if( out && print_tests >= PRINT_BASIC )
00172     *out << "\n1) Check the partitioning ranges and vector space dimensions ...";
00173   lresult = true;
00174 
00175   if( out && print_tests >= PRINT_MORE )
00176     *out << "\n\n1.a) check: equ_decomp.size() + equ_undecomp.size() == ds.m() : ";
00177   result = equ_decomp.size() + equ_undecomp.size() == ds.m();
00178   if(out && print_tests >= PRINT_MORE)
00179     *out << ( result ? "passed" : "failed" );
00180   if(!result) lresult = false;
00181 
00182   if( out && print_tests >= PRINT_MORE )
00183     *out << "\n\n1.b) check: equ_decomp.size() == ds.r() : ";
00184   result = equ_decomp.size() == ds.r();
00185   if(out && print_tests >= PRINT_MORE)
00186     *out << ( result ? "passed" : "failed" );
00187   if(!result) lresult = false;
00188 
00189   if( out && print_tests >= PRINT_MORE )
00190     *out << "\n\n1.c) check: ds.space_range()->dim() == ds.r() : ";
00191   result = ds.space_range()->dim() == ds.r();
00192   if(out && print_tests >= PRINT_MORE)
00193     *out << ( result ? "passed" : "failed" );
00194   if(!result) lresult = false;
00195 
00196   if( out && print_tests >= PRINT_MORE )
00197     *out << "\n\n1.d) check: ds.space_null()->dim() == ds.n() - ds.r() : ";
00198   result = ds.space_null()->dim() == ds.n() - ds.r();
00199   if(out && print_tests >= PRINT_MORE)
00200     *out << ( result ? "passed" : "failed" );
00201   if(!result) lresult = false;
00202 
00203   if(out && print_tests >= PRINT_MORE)
00204     *out << std::endl;
00205 
00206   if(!lresult) success = false;
00207   if( out && print_tests == PRINT_BASIC )
00208     *out << " : " << ( lresult ? "passed" : "failed" );
00209 
00210   //
00211   // Perform the tests
00212   //
00213 
00214   if(out && print_tests >= PRINT_BASIC)
00215     *out
00216       << "\n2) Check the compatibility of the vector spaces for Gc, Z, Y, R, Uz and Uy  ...";
00217   lresult = true;
00218   
00219   if(Z) {
00220     if(out && print_tests >= PRINT_MORE)
00221       *out
00222         << "\n2.a) Check consistency of the vector spaces for:"
00223         << "\n    Z.space_cols() == Gc.space_cols() and Z.space_rows() == ds.space_null()";
00224     llresult = true;
00225     if(out && print_tests >= PRINT_ALL)
00226       *out << "\n\n2.a.1) Z->space_cols().is_compatible(Gc.space_cols()) == true : ";
00227     result = Z->space_cols().is_compatible(Gc.space_cols());  
00228     if(out && print_tests >= PRINT_ALL)
00229       *out << ( result ? "passed" : "failed" )
00230          << std::endl;
00231     if(!result) llresult = false;
00232     if(out && print_tests >= PRINT_ALL)
00233       *out << "\n\n2.a.2) Z->space_cols().is_compatible(*ds.space_null()) == true : ";
00234     result = Z->space_rows().is_compatible(*ds.space_null()); 
00235     if(out && print_tests >= PRINT_ALL)
00236       *out << ( result ? "passed" : "failed" )
00237          << std::endl;
00238     if(!result) llresult = false;
00239     if(!llresult) lresult = false;
00240     if( out && print_tests == PRINT_MORE )
00241       *out << " : " << ( llresult ? "passed" : "failed" );
00242   }
00243 
00244   if(Y) {
00245     if(out && print_tests >= PRINT_MORE)
00246       *out
00247         << "\n2.b) Check consistency of the vector spaces for:"
00248         << "\n    Y.space_cols() == Gc.space_cols() and Y.space_rows() == ds.space_range()";
00249     llresult = true;
00250     if(out && print_tests >= PRINT_ALL)
00251       *out << "\n\n2.b.1) Y->space_cols().is_compatible(Gc.space_cols()) == true : ";
00252     result = Y->space_cols().is_compatible(Gc.space_cols());  
00253     if(out && print_tests >= PRINT_ALL)
00254       *out << ( result ? "passed" : "failed" )
00255          << std::endl;
00256     if(!result) llresult = false;
00257     if(out && print_tests >= PRINT_ALL)
00258       *out << "\n\n2.b.2) Y->space_cols().is_compatible(*ds.space_range()) == true : ";
00259     result = Y->space_rows().is_compatible(*ds.space_range());  
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(R) {
00270     if(out && print_tests >= PRINT_MORE)
00271       *out
00272         << "\n2.c) Check consistency of the vector spaces for:"
00273         << "\n    R.space_cols() == Gc.space_cols()(equ_decomp) and R.space_rows() == ds.space_range()";
00274     llresult = true;
00275     if(out && print_tests >= PRINT_ALL)
00276       *out << "\n\n2.c.1) R->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_decomp)) == true : ";
00277     result = R->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_decomp)); 
00278     if(out && print_tests >= PRINT_ALL)
00279       *out << ( result ? "passed" : "failed" )
00280          << std::endl;
00281     if(!result) llresult = false;
00282     if(out && print_tests >= PRINT_ALL)
00283       *out << "\n\n2.c.2) R->space_cols().is_compatible(*ds.space_range()) == true : ";
00284     result = R->space_rows().is_compatible(*ds.space_range());  
00285     if(out && print_tests >= PRINT_ALL)
00286       *out << ( result ? "passed" : "failed" )
00287          << std::endl;
00288     if(!result) llresult = false;
00289     if(!llresult) lresult = false;
00290     if( out && print_tests == PRINT_MORE )
00291       *out << " : " << ( llresult ? "passed" : "failed" );
00292   }
00293 
00294   if(Uz) {
00295     if(out && print_tests >= PRINT_MORE)
00296       *out
00297         << "\n2.d) Check consistency of the vector spaces for:"
00298         << "\n    Uz.space_cols() == Gc.space_cols()(equ_undecomp) and Uz.space_rows() == ds.space_null()";
00299     llresult = true;
00300     if(out && print_tests >= PRINT_ALL)
00301       *out << "\n\n2.d.1) Uz->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp)) == true : ";
00302     result = Uz->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp));  
00303     if(out && print_tests >= PRINT_ALL)
00304       *out << ( result ? "passed" : "failed" )
00305          << std::endl;
00306     if(!result) llresult = false;
00307     if(out && print_tests >= PRINT_ALL)
00308       *out << "\n\n2.d.2) Uz->space_cols().is_compatible(*ds.space_null()) == true : ";
00309     result = Uz->space_rows().is_compatible(*ds.space_null());  
00310     if(out && print_tests >= PRINT_ALL)
00311       *out << ( result ? "passed" : "failed" )
00312          << std::endl;
00313     if(!result) llresult = false;
00314     if(!llresult) lresult = false;
00315     if( out && print_tests == PRINT_MORE )
00316       *out << " : " << ( llresult ? "passed" : "failed" );
00317   }
00318 
00319   if(Uy) {
00320     if(out && print_tests >= PRINT_MORE)
00321       *out
00322         << "\n2.e) Check consistency of the vector spaces for:"
00323         << "\n    Uy.space_cols() == Gc.space_cols()(equ_undecomp) and Uy.space_rows() == ds.space_range()";
00324     llresult = true;
00325     if(out && print_tests >= PRINT_ALL)
00326       *out << "\n\n2.e.1) Uy->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp)) == true : ";
00327     result = Uy->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp));  
00328     if(out && print_tests >= PRINT_ALL)
00329       *out << ( result ? "passed" : "failed" )
00330          << std::endl;
00331     if(!result) llresult = false;
00332     if(out && print_tests >= PRINT_ALL)
00333       *out << "\n\n2.e.2) Uy->space_cols().is_compatible(*ds.space_range()) == true : ";
00334     result = Uy->space_rows().is_compatible(*ds.space_range()); 
00335     if(out && print_tests >= PRINT_ALL)
00336       *out << ( result ? "passed" : "failed" )
00337          << std::endl;
00338     if(!result) llresult = false;
00339     if(!llresult) lresult = false;
00340     if( out && print_tests == PRINT_MORE )
00341       *out << " : " << ( llresult ? "passed" : "failed" );
00342   }
00343 
00344   if(!lresult) success = false;
00345   if( out && print_tests == PRINT_BASIC )
00346     *out << " : " << ( lresult ? "passed" : "failed" );
00347 
00348   if(out && print_tests >= PRINT_BASIC)
00349     *out
00350       << "\n3) Check the compatibility of the matrices Gc, Z, Y, R, Uz and Uy numerically ...";
00351   
00352   if(Z) {
00353 
00354     if(out && print_tests >= PRINT_MORE)
00355       *out
00356         << std::endl
00357         << "\n3.a) Check consistency of:"
00358         << "\n     op ( alpha* Gc(:,equ_decomp)' * beta*Z ) * v"
00359         << "\n         \\__________________________________/"
00360         << "\n                         A"
00361         << "\n    ==  op( alpha*beta*Uz * v"
00362         << "\n            \\___________/"
00363         << "\n                     B"
00364         << "\nfor random vectors v ...";
00365 
00366     VectorSpace::vec_mut_ptr_t
00367       v_c       = Gc.space_rows().create_member(),
00368       v_c_tmp   = v_c->space().create_member(),
00369       v_x       = Gc.space_cols().create_member(),
00370       v_x_tmp   = v_x->space().create_member(),
00371       v_z       = ds.space_null()->create_member(),
00372       v_z_tmp   = v_z->space().create_member();
00373     
00374     if(out && print_tests >= PRINT_MORE)
00375       *out << "\n\n3.a.1) Testing non-transposed A*v == B*v ...";
00376     if(out && print_tests > PRINT_MORE)
00377       *out << std::endl;
00378     llresult = true;
00379     {for( int k = 1; k <= num_random_tests(); ++k ) {
00380       random_vector( rand_y_l, rand_y_u, v_z.get() );
00381       if(out && print_tests >= PRINT_ALL) {
00382         *out
00383           << "\n3.a.1."<<k<<") random vector " << k << " ( ||v_z||_1 / n = " << (v_z->norm_1() / v_z->dim()) << " )\n";
00384         if(dump_all() && print_tests >= PRINT_ALL)
00385           *out << "\nv_z =\n" << *v_z;
00386       }
00387       V_StMtV( v_x.get(), beta, *Z, no_trans, *v_z );
00388       V_StMtV( v_c.get(), alpha, Gc, trans, *v_x );
00389       *v_c_tmp->sub_view(equ_decomp) = 0.0;
00390       if(equ_undecomp.size()) {
00391         if(Uz)
00392           V_StMtV( v_c_tmp->sub_view(equ_undecomp).get(), alpha*beta, *Uz, no_trans, *v_z );
00393         else
00394           *v_c_tmp->sub_view(equ_undecomp).get() = *v_c->sub_view(equ_undecomp);
00395       }
00396       const value_type
00397         sum_Bv  = sum(*v_c_tmp), // should be zero if equ_undecomp.size() == 0 so scale by 1.0
00398         sum_Av  = sum(*v_c);
00399       assert_print_nan_inf(sum_Bv, "sum(B*v_z)",true,out);
00400       assert_print_nan_inf(sum_Av, "sum(A*v_z)",true,out);
00401       const value_type
00402         calc_err = ::fabs( ( sum_Av - sum_Bv )
00403                    /( ::fabs(sum_Av) + ::fabs(sum_Bv) + (equ_undecomp.size() ? small_num : 1.0) ) );
00404       if(out && print_tests >= PRINT_ALL)
00405         *out
00406           << "\nrel_err(sum(A*v_z),sum(B*v_z)) = "
00407           << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00408           << calc_err << std::endl;
00409       if( calc_err >= mult_warning_tol() ) {
00410         if(out && print_tests >= PRINT_ALL)
00411           *out
00412             << std::endl
00413             << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
00414             << ", rel_err(sum(A*v_z),sum(B*v_z)) = "
00415             << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00416             << calc_err
00417             << " exceeded "
00418             << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
00419             << " = "
00420             << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
00421             << std::endl;
00422         if(calc_err >= mult_error_tol()) {
00423           if(dump_all() && print_tests >= PRINT_ALL) {
00424             *out << "\nalpha = " << alpha << std::endl;
00425             *out << "\nbeta  = " << beta  << std::endl;
00426             *out << "\nv_z =\n"           << *v_z;
00427             *out << "\nbeta*Z*v_z =\n"    << *v_x;
00428             *out << "\nA*v_z =\n"         << *v_c;
00429             *out << "\nB*v_z =\n"         << *v_c_tmp;
00430           }
00431           llresult = false;
00432         }
00433       }
00434     }}
00435     if(!llresult) lresult = false;
00436     if( out && print_tests == PRINT_MORE )
00437       *out << " : " << ( llresult ? "passed" : "failed" )
00438          << std::endl;
00439     
00440     if(out && print_tests >= PRINT_MORE)
00441       *out << "\n\n3.a.2) Testing transposed A'*v == B'*v ...";
00442     if(out && print_tests > PRINT_MORE)
00443       *out << std::endl;
00444     llresult = true;
00445     {for( int k = 1; k <= num_random_tests(); ++k ) {
00446       random_vector( rand_y_l, rand_y_u, v_c.get() );
00447       if(out && print_tests >= PRINT_ALL) {
00448         *out
00449           << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_c||_1 / n = " << (v_c->norm_1() / v_c->dim()) << " )\n";
00450         if(dump_all() && print_tests >= PRINT_ALL)
00451           *out << "\nv_c =\n" << *v_c;
00452       }
00453       V_StMtV( v_x.get(), alpha, Gc, no_trans, *v_c );
00454       V_StMtV( v_z.get(), beta,  *Z, trans,    *v_x );
00455       *v_z_tmp = 0.0;
00456       if(equ_undecomp.size()) {
00457         if(Uz)
00458           V_StMtV( v_z_tmp.get(), alpha*beta, *Uz, trans, *v_c->sub_view(equ_undecomp) );
00459         else
00460           *v_z_tmp = *v_z;
00461       }
00462       const value_type
00463         sum_Bv  = sum(*v_z_tmp), // should be zero so scale by 1.0
00464         sum_Av  = sum(*v_z);
00465       assert_print_nan_inf(sum_Bv, "sum(B'*v_c)",true,out);
00466       assert_print_nan_inf(sum_Av, "sum(A'*v_c)",true,out);
00467       const value_type
00468         calc_err = ::fabs( ( sum_Av - sum_Bv )
00469                    /( ::fabs(sum_Av) + ::fabs(sum_Bv) + (equ_undecomp.size() ? small_num : 1.0) ) );
00470       if(out && print_tests >= PRINT_ALL)
00471         *out
00472           << "\nrel_err(sum(A'*v_c),sum(B'*v_c)) = "
00473           << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00474           << calc_err << std::endl;
00475       if( calc_err >= mult_warning_tol() ) {
00476         if(out && print_tests >= PRINT_ALL)
00477           *out
00478             << std::endl
00479             << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
00480             << ", rel_err(sum(A'*v_c),sum(B'*v_c)) = "
00481             << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00482             << calc_err
00483             << " exceeded "
00484             << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
00485             << " = "
00486             << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
00487             << std::endl;
00488         if(calc_err >= mult_error_tol()) {
00489           if(dump_all() && print_tests >= PRINT_ALL) {
00490             *out << "\nalpha = " << alpha << std::endl;
00491             *out << "\nbeta  = " << beta  << std::endl;
00492             *out << "\nv_c =\n"           << *v_c;
00493             *out << "\nalpha*Gc*v_c =\n"  << *v_x;
00494             *out << "\nA'*v_c =\n"        << *v_z;
00495             *out << "\nB'*v_c =\n"        << *v_z_tmp;
00496           }
00497           llresult = false;
00498         }
00499       }
00500     }}
00501     if(!llresult) lresult = false;
00502     if( out && print_tests == PRINT_MORE )
00503       *out << " : " << ( llresult ? "passed" : "failed" )
00504          << std::endl;
00505 
00506   }
00507   else {
00508     if(out && print_tests >= PRINT_MORE)
00509       *out
00510         << std::endl
00511         << "\n3.a) Warning! Z ==NULL; Z, and Uz are not checked numerically ...\n";
00512   }
00513 
00514   if(Y) {
00515 
00516     if(out && print_tests >= PRINT_MORE)
00517       *out
00518         << std::endl
00519         << "\n3.b) Check consistency of:"
00520         << "\n     op ( alpha*[ Gc(:,equ_decomp)'   ]"
00521         << "\n                [ Gc(:,equ_undecomp)' ] * beta*Y ) * v"
00522         << "\n         \\_____________________________________/"
00523         << "\n                         A"
00524         << "\n    ==  op( alpha*beta*[ R  ]"
00525         << "\n                       [ Uy ] ) * v"
00526         << "\n            \\_______________/"
00527         << "\n                     B"
00528         << "\nfor random vectors v ...";
00529 
00530     VectorSpace::vec_mut_ptr_t
00531       v_c       = Gc.space_rows().create_member(),
00532       v_c_tmp   = v_c->space().create_member(),
00533       v_x       = Gc.space_cols().create_member(),
00534       v_x_tmp   = v_x->space().create_member(),
00535       v_y       = ds.space_range()->create_member(),
00536       v_y_tmp   = v_y->space().create_member();
00537     
00538     if(out && print_tests >= PRINT_MORE)
00539       *out << "\n\n3.b.1) Testing non-transposed A*v == B*v ...";
00540     if(out && print_tests > PRINT_MORE)
00541       *out << std::endl;
00542     llresult = true;
00543     {for( int k = 1; k <= num_random_tests(); ++k ) {
00544       random_vector( rand_y_l, rand_y_u, v_y.get() );
00545       if(out && print_tests >= PRINT_ALL) {
00546         *out
00547           << "\n3.b.1."<<k<<") random vector " << k << " ( ||v_y||_1 / n = " << (v_y->norm_1() / v_y->dim()) << " )\n";
00548         if(dump_all() && print_tests >= PRINT_ALL)
00549           *out << "\nv_y =\n" << *v_y;
00550       }
00551       V_StMtV( v_x.get(), beta, *Y, no_trans, *v_y );
00552       V_StMtV( v_c.get(), alpha, Gc, trans, *v_x );
00553       V_StMtV( v_c_tmp->sub_view(equ_decomp).get(), alpha*beta, *R, no_trans, *v_y );
00554       if(equ_undecomp.size()) {
00555         if(Uy)
00556           V_StMtV( v_c_tmp->sub_view(equ_undecomp).get(), alpha*beta, *Uy, no_trans, *v_y );
00557         else
00558           *v_c_tmp->sub_view(equ_undecomp) = *v_c->sub_view(equ_undecomp);
00559       }
00560       const value_type
00561         sum_Bv  = sum(*v_c_tmp),
00562         sum_Av  = sum(*v_c);
00563       assert_print_nan_inf(sum_Bv, "sum(B*v_y)",true,out);
00564       assert_print_nan_inf(sum_Av, "sum(A*v_y)",true,out);
00565       const value_type
00566         calc_err = ::fabs( ( sum_Av - sum_Bv )
00567                    /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
00568       if(out && print_tests >= PRINT_ALL)
00569         *out
00570           << "\nrel_err(sum(A*v_y),sum(B*v_y)) = "
00571           << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00572           << calc_err << std::endl;
00573       if( calc_err >= mult_warning_tol() ) {
00574         if(out && print_tests >= PRINT_ALL)
00575           *out
00576             << std::endl
00577             << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
00578             << ", rel_err(sum(A*v_y),sum(B*v_y)) = "
00579             << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00580             << calc_err
00581             << " exceeded "
00582             << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
00583             << " = "
00584             << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
00585             << std::endl;
00586         if(calc_err >= mult_error_tol()) {
00587           if(dump_all() && print_tests >= PRINT_ALL) {
00588             *out << "\nalpha = " << alpha << std::endl;
00589             *out << "\nbeta  = " << beta  << std::endl;
00590             *out << "\nv_y =\n"           << *v_y;
00591             *out << "\nbeta*Y*v_y =\n"    << *v_x;
00592             *out << "\nA*v_y =\n"         << *v_c;
00593             *out << "\nB*v_y =\n"         << *v_c_tmp;
00594           }
00595           llresult = false;
00596         }
00597       }
00598     }}
00599     if(!llresult) lresult = false;
00600     if( out && print_tests == PRINT_MORE )
00601       *out << " : " << ( llresult ? "passed" : "failed" )
00602          << std::endl;
00603     
00604     if(out && print_tests >= PRINT_MORE)
00605       *out << "\n\n3.b.2) Testing transposed A'*v == B'*v ...";
00606     if(out && print_tests > PRINT_MORE)
00607       *out << std::endl;
00608     llresult = true;
00609     {for( int k = 1; k <= num_random_tests(); ++k ) {
00610       random_vector( rand_y_l, rand_y_u, v_c.get() );
00611       if(out && print_tests >= PRINT_ALL) {
00612         *out
00613           << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_c||_1 / n = " << (v_c->norm_1() / v_c->dim()) << " )\n";
00614         if(dump_all() && print_tests >= PRINT_ALL)
00615           *out << "\nv_c =\n" << *v_c;
00616       }
00617       V_StMtV( v_x.get(), alpha, Gc, no_trans, *v_c );
00618       V_StMtV( v_y.get(), beta,  *Y, trans,    *v_x );
00619       V_StMtV( v_y_tmp.get(), alpha*beta, *R, trans, *v_c->sub_view(equ_decomp) );
00620       if(equ_undecomp.size()) {
00621         if(Uy)
00622           Vp_StMtV( v_y_tmp.get(), alpha*beta, *Uy, trans, *v_c->sub_view(equ_undecomp) );
00623         else
00624           Vp_V( v_y_tmp.get(), *v_y );
00625       }
00626       const value_type
00627         sum_Bv  = sum(*v_y_tmp), // should be zero so scale by 1.0
00628         sum_Av  = sum(*v_y);
00629       assert_print_nan_inf(sum_Bv, "sum(B'*v_c)",true,out);
00630       assert_print_nan_inf(sum_Av, "sum(A'*v_c)",true,out);
00631       const value_type
00632         calc_err = ::fabs( ( sum_Av - sum_Bv )
00633                    /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
00634       if(out && print_tests >= PRINT_ALL)
00635         *out
00636           << "\nrel_err(sum(A'*v_c),sum(B'*v_c)) = "
00637           << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00638           << calc_err << std::endl;
00639       if( calc_err >= mult_warning_tol() ) {
00640         if(out && print_tests >= PRINT_ALL)
00641           *out
00642             << std::endl
00643             << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
00644             << ", rel_err(sum(A'*v_c),sum(B'*v_c)) = "
00645             << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
00646             << calc_err
00647             << " exceeded "
00648             << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
00649             << " = "
00650             << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
00651             << std::endl;
00652         if(calc_err >= mult_error_tol()) {
00653           if(dump_all() && print_tests >= PRINT_ALL) {
00654             *out << "\nalpha = " << alpha << std::endl;
00655             *out << "\nbeta  = " << beta  << std::endl;
00656             *out << "\nv_c =\n"           << *v_c;
00657             *out << "\nalpha*Gc*v_c =\n"  << *v_x;
00658             *out << "\nA'*v_c =\n"        << *v_y;
00659             *out << "\nB'*v_c =\n"        << *v_y_tmp;
00660           }
00661           llresult = false;
00662         }
00663       }
00664     }}
00665     if(!llresult) lresult = false;
00666     if( out && print_tests == PRINT_MORE )
00667       *out << " : " << ( llresult ? "passed" : "failed" )
00668          << std::endl;
00669 
00670   }
00671   else {
00672     if(out && print_tests >= PRINT_MORE)
00673       *out
00674         << std::endl
00675         << "\n3.b) Warning! Y ==NULL; Y, R and Uy are not checked numerically ...\n";
00676   }
00677 
00678   if(R) {
00679     if(out && print_tests >= PRINT_MORE)
00680       *out
00681         << std::endl
00682         << "\n3.b) Check consistency of: op(op(inv(R))*op(R)) == I ...\n";
00683     typedef MatrixOpNonsingTester  MWONST_t;
00684     MWONST_t::EPrintTestLevel
00685       olevel;
00686     switch(print_tests) {
00687       case PRINT_NONE:
00688       case PRINT_BASIC:
00689         olevel = MWONST_t::PRINT_NONE;
00690         break;
00691       case PRINT_MORE:
00692         olevel = MWONST_t::PRINT_MORE;
00693         break;
00694       case PRINT_ALL:
00695         olevel = MWONST_t::PRINT_ALL;
00696         break;
00697       default:
00698         TEST_FOR_EXCEPT(true); // Should not get here
00699     }
00700     MWONST_t
00701       R_tester(
00702         MWONST_t::TEST_LEVEL_2_BLAS
00703         ,olevel
00704         ,dump_all()
00705         ,throw_exception()
00706         ,num_random_tests()
00707         ,solve_warning_tol()
00708         ,solve_error_tol()
00709         );
00710     lresult = R_tester.test_matrix(*R,"R",out);
00711   }
00712 
00713   if(!lresult) success = false;
00714   if( out && print_tests == PRINT_BASIC )
00715     *out << " : " << ( lresult ? "passed" : "failed" );
00716   
00717   if( out && print_tests != PRINT_NONE ) {
00718     if(success)
00719       *out << "\nCongradulations! The DecompositionSystem object and its associated matrix objects seem to check out!\n";
00720     else
00721       *out << "\nOops! At least one of the tests did not check out!\n";
00722     if( print_tests >= PRINT_BASIC )
00723       *out << "\nEnd DecompositionSystemTester::test_decomp_system(...)\n";
00724   }
00725 
00726   return success;
00727 }
00728 
00729 } // end namespace ConstrainedOptPack

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