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