|
ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
|
|
ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
|
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
1.7.4
