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