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