MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_BasisSystemTester.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <math.h>
43 
44 #include <ostream>
45 #include <limits>
46 
58 
59 namespace AbstractLinAlgPack {
60 
62  EPrintTestLevel print_tests
63  ,bool dump_all
64  ,bool throw_exception
65  ,size_type num_random_tests
66  ,value_type warning_tol
67  ,value_type error_tol
68  )
69  :print_tests_(print_tests)
70  ,dump_all_(dump_all)
71  ,throw_exception_(throw_exception)
72  ,num_random_tests_(num_random_tests)
73  ,warning_tol_(warning_tol)
74  ,error_tol_(error_tol)
75 {}
76 
78  const BasisSystem &bs
79  ,const MatrixOp *Gc
80  ,const MatrixOpNonsing *C
81  ,const MatrixOp *N_in
82  ,const MatrixOp *D
83  ,const MatrixOp *GcUP
84  ,std::ostream *out
85  )
86 {
87  namespace rcp = MemMngPack;
88  using BLAS_Cpp::no_trans;
89  using BLAS_Cpp::trans;
96  using LinAlgOpPack::V_MtV;
97  using LinAlgOpPack::V_StV;
98  using LinAlgOpPack::V_VpV;
99  using LinAlgOpPack::Vp_V;
100 
101  // ToDo: Check the preconditions
102 
103  bool success = true, result, lresult, llresult;
104  const value_type
105  rand_y_l = -1.0,
106  rand_y_u = 1.0,
107  small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
108  alpha = 2.0;
109 
111  print_tests = ( this->print_tests() == PRINT_NOT_SELECTED ? PRINT_NONE : this->print_tests() );
112 
114 
115  // Print the input?
116  if( out && print_tests != PRINT_NONE ) {
117  if( print_tests >= PRINT_BASIC ) {
118  *out << "\n*************************************************"
119  << "\n*** BasisSystemTester::test_basis_system(...) ***"
120  << "\n*************************************************\n";
121  if(Gc)
122  *out << "\n||Gc||inf = " << Gc->calc_norm(mat_nrm_inf).value;
123  if(C) {
124  *out << "\n||C||inf = " << C->calc_norm(mat_nrm_inf).value;
125  *out << "\ncond_inf(C) = " << C->calc_cond_num(mat_nrm_inf).value;
126  }
127  if(N_in)
128  *out << "\n||N||inf = " << N_in->calc_norm(mat_nrm_inf).value;
129  if(D)
130  *out << "\n||D||inf = " << D->calc_norm(mat_nrm_inf).value;
131  if(GcUP)
132  *out << "\n||GcUP||inf = " << GcUP->calc_norm(mat_nrm_inf).value;
133  }
134  if(dump_all()) {
135  if(Gc)
136  *out << "\nGc =\n" << *Gc;
137  if(C)
138  *out << "\nC =\n" << *C;
139  if(N_in)
140  *out << "\nN =\n" << *N_in;
141  if(D)
142  *out << "\nD =\n" << *D;
143  if(GcUP)
144  *out << "\nGcUP =\n" << *GcUP;
145  }
146  }
147 
148  //
149  // Check the dimensions of everything
150  //
151 
152  const Range1D
153  var_dep = bs.var_dep(),
154  var_indep = bs.var_indep(),
155  equ_decomp = bs.equ_decomp(),
156  equ_undecomp = bs.equ_undecomp();
157 
158  if( out && print_tests >= PRINT_MORE ) {
159  *out
160  << "\nbs.var_dep() = ["<<var_dep.lbound()<<","<<var_dep.ubound()<<"]"
161  << "\nbs.var_indep( ) = ["<<var_indep.lbound()<<","<<var_indep.ubound()<<"]"
162  << "\nbs.equ_decomp() = ["<<equ_decomp.lbound()<<","<<equ_decomp.ubound()<<"]"
163  << "\nbs.equ_undecomp() = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
164  << std::endl;
165  }
166 
167  if( out && print_tests >= PRINT_BASIC )
168  *out << "\n1) Check the partitioning ranges ...";
169  lresult = true;
170 
171  if( out && print_tests >= PRINT_MORE )
172  *out << "\n\n1.a) check: var_dep.size() != equ_decomp.size() : ";
173  result = var_dep.size() == equ_decomp.size();
174  if(out && print_tests >= PRINT_MORE)
175  *out << ( result ? "passed" : "failed" );
176  if(!result) lresult = false;
177 
178  if(Gc) {
179  if( out && print_tests >= PRINT_MORE )
180  *out << "\n1.b) check: var_dep.size() + var_indep.size() == Gc->rows() : ";
181  result = var_dep.size() + var_indep.size() == Gc->rows();
182  if(out && print_tests >= PRINT_MORE)
183  *out << ( result ? "passed" : "failed" );
184  if(!result) lresult = false;
185  }
186 
187  if(Gc) {
188  if( out && print_tests >= PRINT_MORE )
189  *out << "\n1.d) check: equ_decomp.size() + equ_undecomp.size() == Gc->cols() : ";
190  result = equ_decomp.size() + equ_undecomp.size() == Gc->cols();
191  if(out && print_tests >= PRINT_MORE)
192  *out << ( result ? "passed" : "failed" );
193  if(!result) lresult = false;
194  }
195 
196  if(out && print_tests >= PRINT_MORE)
197  *out << std::endl;
198 
199  if(!lresult) success = false;
200  if( out && print_tests == PRINT_BASIC )
201  *out << " : " << ( lresult ? "passed" : "failed" );
202 
203  // Create the N matrix if not input
205  N = Teuchos::rcp(N_in,false);
206  if( Gc && C && N.get() == NULL ) {
207  if(out && print_tests >= PRINT_BASIC)
208  *out
209  << "\nCreating the matrix N since it was not input by the client ...";
210  if(out && print_tests >= PRINT_MORE)
211  *out
212  << std::endl;
214  N_comp = Teuchos::rcp(new AbstractLinAlgPack::MatrixComposite(var_dep.size(),var_indep.size()));
215  if( equ_decomp.size() )
216  N_comp->add_matrix( 0, 0, 1.0, equ_decomp, Gc, Teuchos::null, BLAS_Cpp::trans, var_indep );
217  N_comp->finish_construction(
218  Gc->space_rows().sub_space(equ_decomp)->clone()
219  ,Gc->space_cols().sub_space(var_indep)->clone()
220  );
221  if( out && dump_all() )
222  *out << "\nN =\n" << *N_comp;
223  N = N_comp;
224  }
225 
226  // Create the other auxillary matrix objects
227  if( equ_undecomp.size() ) {
228  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Create matrix objects for Gc(var_dep,equ_undecomp) and Gc(var_indep,equ_undecomp)
229  }
230 
231  //
232  // Perform the tests
233  //
234 
235  if( C && N.get() ) {
236 
237  if(out && print_tests >= PRINT_BASIC)
238  *out
239  << "\n2) Check the compatibility of the vector spaces for C, N, D and Gc ...";
240  lresult = true;
241 
242  if(out && print_tests >= PRINT_MORE)
243  *out
244  << "\n2.a) Check consistency of the vector spaces for:"
245  << "\n C.space_cols() == N.space_cols()";
246  llresult = true;
247  if(out && print_tests >= PRINT_ALL)
248  *out << "\n\n2.a.1) C->space_cols().is_compatible(N->space_cols()) == true : ";
249  result = C->space_cols().is_compatible(N->space_cols());
250  if(out && print_tests >= PRINT_ALL)
251  *out << ( result ? "passed" : "failed" )
252  << std::endl;
253  if(!result) llresult = false;
254  if(!llresult) lresult = false;
255  if( out && print_tests == PRINT_MORE )
256  *out << " : " << ( llresult ? "passed" : "failed" );
257 
258  if(D) {
259  if(out && print_tests >= PRINT_MORE)
260  *out
261  << "\n2.b) Check consistency of the vector spaces for:"
262  << "\n D.space_cols() == C.space_cols() and D.space_rows() == N.space_rows()";
263  llresult = true;
264  if(out && print_tests >= PRINT_ALL)
265  *out << "\n2.b.1) D->space_cols().is_compatible(C->space_cols()) == true : ";
266  result = D->space_cols().is_compatible(C->space_cols());
267  if(out && print_tests >= PRINT_ALL)
268  *out << ( result ? "passed" : "failed" );
269  if(!result) llresult = false;
270  if(out && print_tests >= PRINT_ALL)
271  *out << "\n2.b.2) D->space_rows().is_compatible(N->space_rows()) == true : ";
272  result = D->space_rows().is_compatible(N->space_rows());
273  if(out && print_tests >= PRINT_ALL)
274  *out << ( result ? "passed" : "failed" )
275  << std::endl;
276  if(!result) llresult = false;
277  if(!llresult) lresult = false;
278  if( out && print_tests == PRINT_MORE )
279  *out << " : " << ( llresult ? "passed" : "failed" );
280  }
281 
282  if(out && print_tests >= PRINT_MORE)
283  *out
284  << "\n2.c) Check consistency of the vector spaces for:"
285  << "\n Gc'(equ_decomp, var_dep) == C";
286  llresult = true;
287  if( equ_decomp.size() ) {
288  if(out && print_tests >= PRINT_ALL)
289  *out << "\n2.c.1) Gc->space_rows().sub_space(equ_decomp)->is_compatible(*C->space_cols().sub_space(equ_decomp)) == true : ";
290  result = Gc->space_rows().sub_space(equ_decomp)->is_compatible(*C->space_cols().sub_space(equ_decomp));
291  if(out && print_tests >= PRINT_ALL)
292  *out << ( result ? "passed" : "failed" );
293  if(!result) llresult = false;
294  if(out && print_tests >= PRINT_ALL)
295  *out << "\n2.c.2) Gc->space_cols().sub_space(var_dep)->is_compatible(C->space_rows()) == true : ";
296  result = Gc->space_cols().sub_space(var_dep)->is_compatible(C->space_rows());
297  if(out && print_tests >= PRINT_ALL)
298  *out << ( result ? "passed" : "failed" );
299  if(!result) llresult = false;
300  }
301  if(out && print_tests >= PRINT_ALL)
302  *out << std::endl;
303  if(!llresult) lresult = false;
304  if( out && print_tests == PRINT_MORE )
305  *out << " : " << ( llresult ? "passed" : "failed" );
306 
307  if(out && print_tests >= PRINT_MORE)
308  *out
309  << "\n2.d) Check consistency of the vector spaces for:"
310  << "\n Gc'(equ_decomp, var_indep) == N";
311  llresult = true;
312  if( equ_decomp.size() ) {
313  if(out && print_tests >= PRINT_ALL)
314  *out << "\n2.d.1) Gc->space_rows().sub_space(equ_decomp)->is_compatible(*N->space_cols().sub_space(equ_decomp)) == true : ";
315  result = Gc->space_rows().sub_space(equ_decomp)->is_compatible(*N->space_cols().sub_space(equ_decomp));
316  if(out && print_tests >= PRINT_ALL)
317  *out << ( result ? "passed" : "failed" );
318  if(!result) llresult = false;
319  if(out && print_tests >= PRINT_ALL)
320  *out << "\n2.d.2) Gc->space_cols().sub_space(var_indep)->is_compatible(N->space_rows()) == true : ";
321  result = Gc->space_cols().sub_space(var_indep)->is_compatible(N->space_rows());
322  if(out && print_tests >= PRINT_ALL)
323  *out << ( result ? "passed" : "failed" );
324  if(!result) llresult = false;
325  }
326  if(out && print_tests >= PRINT_ALL)
327  *out << std::endl;
328  if(!llresult) lresult = false;
329  if( out && print_tests == PRINT_MORE )
330  *out << " : " << ( llresult ? "passed" : "failed" )
331  << std::endl;
332 
333  if(!lresult) success = false;
334  if( out && print_tests == PRINT_BASIC )
335  *out << " : " << ( lresult ? "passed" : "failed" );
336 
337  if(out && print_tests >= PRINT_BASIC)
338  *out
339  << "\n3) Check the compatibility of the matrices C, N, D and Gc numerically ...";
340 
341  if(out && print_tests >= PRINT_MORE)
342  *out
343  << std::endl
344  << "\n3.a) Check consistency of:"
345  << "\n "
346  << "\n op ( alpha* [ Gc'(equ_decomp, var_dep) Gc'(equ_decomp, var_indep) ] ) * v"
347  << "\n \\______________________________________________________________/"
348  << "\n A"
349  << "\n == op( alpha*[ C N ] ) * v"
350  << "\n \\____________/"
351  << "\n B"
352  << "\nfor random vectors v ...";
353  {
354 
356  Gc_v_x = Gc->space_cols().create_member(),
357  Gc_v_c = Gc->space_rows().create_member(),
358  C_v_xD = C->space_rows().create_member(),
359  C_v_chD = C->space_cols().create_member(),
360  N_v_xI = N->space_rows().create_member(),
361  N_v_chD = N->space_cols().create_member(),
362  v_x = Gc->space_cols().create_member(),
363  v_x_tmp = v_x->space().create_member(),
364  v_chD = C_v_xD->space().create_member(),
365  v_chD_tmp = v_chD->space().create_member();
366 
367  if(out && print_tests >= PRINT_MORE)
368  *out << "\n\n3.a.1) Testing non-transposed A*v == B*v ...";
369  if(out && print_tests > PRINT_MORE)
370  *out << std::endl;
371  llresult = true;
372  {for( int k = 1; k <= num_random_tests(); ++k ) {
373  random_vector( rand_y_l, rand_y_u, v_x.get() );
374  if(out && print_tests >= PRINT_ALL) {
375  *out
376  << "\n3.a.1."<<k<<") random vector " << k << " ( ||v_x||_1 / n = " << (v_x->norm_1() / v_x->dim()) << " )\n";
377  if(dump_all() && print_tests >= PRINT_ALL)
378  *out << "\nv_x =\n" << *v_x;
379  }
380  if(Gc && equ_decomp.size()) {
381  V_StMtV( Gc_v_c.get(), alpha, *Gc, trans, *v_x );
382  *v_chD_tmp->sub_view(equ_decomp)
383  = *Gc_v_c->sub_view(equ_decomp);
384  }
385  V_StMtV( C_v_chD.get(), alpha, *C, no_trans, *v_x->sub_view(var_dep) );
386  V_StMtV( N_v_chD.get(), alpha, *N, no_trans, *v_x->sub_view(var_indep) );
387  V_VpV( v_chD.get(), *C_v_chD, *N_v_chD );
388  const value_type
389  sum_Bv = sum(*v_chD),
390  sum_Av = sum(*v_chD_tmp);
391  assert_print_nan_inf(sum_Bv, "sum(B*v_x)",true,out);
392  assert_print_nan_inf(sum_Av, "sum(A*v_x)",true,out);
393  const value_type
394  calc_err = ::fabs( ( sum_Av - sum_Bv )
395  /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
396  if(out && print_tests >= PRINT_ALL)
397  *out
398  << "\nrel_err(sum(A*v_x),sum(B*v_x)) = "
399  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
400  << calc_err << std::endl;
401  if( calc_err >= warning_tol() ) {
402  if(out && print_tests >= PRINT_ALL)
403  *out
404  << std::endl
405  << ( calc_err >= error_tol() ? "Error" : "Warning" )
406  << ", rel_err(sum(A*v_x),sum(B*v_x)) = "
407  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
408  << calc_err
409  << " exceeded "
410  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
411  << " = "
412  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
413  << std::endl;
414  if(calc_err >= error_tol()) {
415  if(dump_all() && print_tests >= PRINT_ALL) {
416  *out << "\nalpha = " << alpha << std::endl;
417  *out << "\nv_x =\n" << *v_x;
418  *out << "\nalpha*Gc*v_x =\n" << *Gc_v_c;
419  *out << "A*v =\n" << *v_chD_tmp;
420  *out << "\nalpha*C*v_x =\n" << *C_v_chD;
421  *out << "\nalpha*N*v_x =\n" << *N_v_chD;
422  *out << "\nB*v =\n" << *v_chD;
423  }
424  llresult = false;
425  }
426  }
427  }}
428  if(!llresult) lresult = false;
429  if( out && print_tests == PRINT_MORE )
430  *out << " : " << ( llresult ? "passed" : "failed" )
431  << std::endl;
432 
433  if(out && print_tests >= PRINT_MORE)
434  *out << "\n3.a.2) Testing transposed A'*v == B'*v ...";
435  if(out && print_tests > PRINT_MORE)
436  *out << std::endl;
437  llresult = true;
438  {for( int k = 1; k <= num_random_tests(); ++k ) {
439  random_vector( rand_y_l, rand_y_u, v_chD.get() );
440  if(out && print_tests >= PRINT_ALL) {
441  *out
442  << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_chD||_1 / n = " << (v_chD->norm_1() / v_chD->dim()) << " )\n";
443  if(dump_all() && print_tests >= PRINT_ALL)
444  *out << "\nv_chD =\n" << *v_chD;
445  }
446  *v_x_tmp = 0.0;
447  if(Gc && equ_decomp.size()) {
448  *Gc_v_c->sub_view(equ_decomp) = *v_chD->sub_view(equ_decomp);
449  if(equ_undecomp.size())
450  *Gc_v_c->sub_view(equ_undecomp) = 0.0;
451  V_StMtV( Gc_v_x.get(), alpha, *Gc, no_trans, *Gc_v_c );
452  Vp_V( v_x_tmp.get(), *Gc_v_x );
453  }
454  V_StMtV( C_v_xD.get(), alpha, *C, trans, *v_chD );
455  *v_x->sub_view(var_dep) = *C_v_xD;
456  V_StMtV( N_v_xI.get(), alpha, *N, trans, *v_chD );
457  *v_x->sub_view(var_indep) = *N_v_xI;
458  const value_type
459  sum_BTv = sum(*v_x),
460  sum_ATv = sum(*v_x_tmp);
461  assert_print_nan_inf(sum_BTv, "sum(B'*v_chD)",true,out);
462  assert_print_nan_inf(sum_ATv, "sum(A'*v_chD)",true,out);
463  const value_type
464  calc_err = ::fabs( ( sum_ATv - sum_BTv )
465  /( ::fabs(sum_ATv) + ::fabs(sum_BTv) + small_num ) );
466  if(out && print_tests >= PRINT_ALL)
467  *out
468  << "\nrel_err(sum(A'*v_chD),sum(B'*v_chD)) = "
469  << "rel_err(" << sum_ATv << "," << sum_BTv << ") = "
470  << calc_err << std::endl;
471  if( calc_err >= warning_tol() ) {
472  if(out && print_tests >= PRINT_ALL)
473  *out
474  << std::endl
475  << ( calc_err >= error_tol() ? "Error" : "Warning" )
476  << ", rel_err(sum(A'*v_chD),sum(B'*v_chD)) = "
477  << "rel_err(" << sum_ATv << "," << sum_BTv << ") = "
478  << calc_err << std::endl
479  << " exceeded "
480  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
481  << " = "
482  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
483  << std::endl;
484  if(calc_err >= error_tol()) {
485  if(dump_all() && print_tests >= PRINT_ALL) {
486  *out << "\nalpha = " << alpha << std::endl;
487  *out << "\nv_chD =\n" << *v_chD;
488  if(Gc_v_x.get() && equ_decomp.size()) {
489  *out << "\nGc_v_c =\n" << *Gc_v_c;
490  *out << "\nalpha*Gc'*[v_chD(equ_decomp); 0] =\n"
491  << *Gc_v_x;
492  }
493  *out << "A'*v =\n" << *v_x_tmp;
494  *out << "\nalpha*C*v_chD =\n" << *C_v_xD;
495  *out << "\nalpha*N*v_chD =\n" << *N_v_xI;
496  *out << "\nB'*v =\n" << *v_x;
497  }
498  llresult = false;
499  }
500  }
501  }}
502  if(!llresult) lresult = false;
503  if( out && print_tests == PRINT_MORE )
504  *out << " : " << ( llresult ? "passed" : "failed" )
505  << std::endl;
506  }
507 
508  if(out && print_tests >= PRINT_MORE)
509  *out
510  << "\n3.b) Check consistency of:"
511  << "\n alpha*op(C)*(op(inv(C)) * v) == alpha*v"
512  << "\nfor random vectors v ...";
513  {
515  v_xD = C->space_rows().create_member(),
516  v_xD_tmp = C->space_rows().create_member(),
517  v_chD = C->space_cols().create_member(),
518  v_chD_tmp = C->space_cols().create_member();
519 
520  if(out && print_tests >= PRINT_MORE)
521  *out << "\n\n3.b.1) Testing non-transposed: alpha*C*(inv(C)*v) == alpha*v ...";
522  if(out && print_tests > PRINT_MORE)
523  *out << std::endl;
524  llresult = true;
525  {for( int k = 1; k <= num_random_tests(); ++k ) {
526  random_vector( rand_y_l, rand_y_u, v_chD.get() );
527  if(out && print_tests >= PRINT_ALL) {
528  *out
529  << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_chD||_1 / n = " << (v_chD->norm_1() / v_chD->dim()) << " )\n";
530  if(dump_all() && print_tests >= PRINT_ALL)
531  *out << "\nv_chD =\n" << *v_chD;
532  }
533  V_InvMtV( v_xD_tmp.get(), *C, no_trans, *v_chD );
534  V_StMtV( v_chD_tmp.get(), alpha, *C, no_trans, *v_xD_tmp );
535  const value_type
536  sum_aCICv = sum(*v_chD_tmp),
537  sum_av = alpha * sum(*v_chD);
538  assert_print_nan_inf(sum_aCICv, "sum(alpha*C*(inv(C)*v)",true,out);
539  assert_print_nan_inf(sum_av, "sum(alpha*v)",true,out);
540  const value_type
541  calc_err = ::fabs( ( sum_aCICv - sum_av )
542  /( ::fabs(sum_aCICv) + ::fabs(sum_av) + small_num ) );
543  if(out && print_tests >= PRINT_ALL)
544  *out
545  << "\nrel_err(sum(alpha*C*(inv(C)*v),sum(alpha*v)) = "
546  << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
547  << calc_err << std::endl;
548  if( calc_err >= warning_tol() ) {
549  if(out && print_tests >= PRINT_ALL)
550  *out
551  << std::endl
552  << ( calc_err >= error_tol() ? "Error" : "Warning" )
553  << ", rel_err(sum(alpha*C*(inv(C)*v)),sum(alpha*v)) = "
554  << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
555  << calc_err
556  << " exceeded "
557  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
558  << " = "
559  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
560  << std::endl;
561  if(calc_err >= error_tol()) {
562  if(dump_all() && print_tests >= PRINT_ALL) {
563  *out << "\nalpha = " << alpha << std::endl;
564  *out << "\nv_chD =\n" << *v_chD;
565  *out << "\ninv(C)*v_chD =\n" << *v_xD_tmp;
566  *out << "\nalpha*C*inv(C)*v_chD =\n" << *v_chD_tmp;
567  }
568  llresult = false;
569  }
570  }
571  }}
572  if(!llresult) lresult = false;
573  if( out && print_tests == PRINT_MORE )
574  *out << " : " << ( llresult ? "passed" : "failed" )
575  << std::endl;
576 
577  if(out && print_tests >= PRINT_MORE)
578  *out << "\n3.b.2) Testing transposed: alpha*C'*(inv(C')*v) == alpha*v ...";
579  if(out && print_tests > PRINT_MORE)
580  *out << std::endl;
581  llresult = true;
582  {for( int k = 1; k <= num_random_tests(); ++k ) {
583  random_vector( rand_y_l, rand_y_u, v_xD.get() );
584  if(out && print_tests >= PRINT_ALL) {
585  *out
586  << "\n3.b.2."<<k<<") random vector " << k << " ( ||v_xD||_1 / n = " << (v_xD->norm_1() / v_xD->dim()) << " )\n";
587  if(dump_all() && print_tests >= PRINT_ALL)
588  *out << "\nv_xD =\n" << *v_xD;
589  }
590  V_InvMtV( v_chD_tmp.get(), *C, trans, *v_xD );
591  V_StMtV( v_xD_tmp.get(), alpha, *C, trans, *v_chD_tmp );
592  const value_type
593  sum_aCICv = sum(*v_xD_tmp),
594  sum_av = alpha * sum(*v_xD);
595  assert_print_nan_inf(sum_aCICv, "sum(alpha*C'*(inv(C')*v)",true,out);
596  assert_print_nan_inf(sum_av, "sum(alpha*v)",true,out);
597  const value_type
598  calc_err = ::fabs( ( sum_aCICv - sum_av )
599  /( ::fabs(sum_aCICv) + ::fabs(sum_av) + small_num ) );
600  if(out && print_tests >= PRINT_ALL)
601  *out
602  << "\nrel_err(sum(alpha*C'*(inv(C')*v)),sum(alpha*v)) = "
603  << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
604  << calc_err << std::endl;
605  if( calc_err >= warning_tol() ) {
606  if(out && print_tests >= PRINT_ALL)
607  *out
608  << std::endl
609  << ( calc_err >= error_tol() ? "Error" : "Warning" )
610  << ", rel_err(sum(alpha*C'*(inv(C')*v)),sum(alpha*v)) = "
611  << "rel_err(" << sum_aCICv << "," << sum_av << ") = "
612  << calc_err
613  << " exceeded "
614  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
615  << " = "
616  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
617  << std::endl;
618  if(calc_err >= error_tol()) {
619  if(dump_all() && print_tests >= PRINT_ALL) {
620  *out << "\nalpha = " << alpha << std::endl;
621  *out << "\nv_xD =\n" << *v_xD;
622  *out << "\ninv(C')*v_xD =\n" << *v_chD_tmp;
623  *out << "\nalpha*C'*inv(C')*v_xD =\n" << *v_xD_tmp;
624  }
625  llresult = false;
626  }
627  }
628  }}
629  if(!llresult) lresult = false;
630  if( out && print_tests == PRINT_MORE )
631  *out << " : " << ( llresult ? "passed" : "failed" )
632  << std::endl;
633  }
634 
635  if(D) {
636  if(out && print_tests >= PRINT_MORE)
637  *out
638  << "\n3.c) Check consistency of:"
639  << "\n alpha * op(-inv(C) * N) * v == alpha * op(D) * v"
640  << "\nfor random vectors v ...";
641 
642  {
644  v_xD = C->space_rows().create_member(),
645  v_xI = N->space_rows().create_member(),
646  v_xD_tmp = C->space_rows().create_member(),
647  v_xI_tmp = N->space_rows().create_member(),
648  v_chD_tmp = C->space_cols().create_member();
649 
650  if(out && print_tests >= PRINT_MORE)
651  *out << "\n\n3.b.1) Testing non-transposed: inv(C)*(-alpha*N*v) == alpha*D*v ...";
652  if(out && print_tests > PRINT_MORE)
653  *out << std::endl;
654  llresult = true;
655  {for( int k = 1; k <= num_random_tests(); ++k ) {
656  random_vector( rand_y_l, rand_y_u, v_xI.get() );
657  if(out && print_tests >= PRINT_ALL) {
658  *out
659  << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_xI||_1 / n = " << (v_xI->norm_1() / v_xI->dim()) << " )\n";
660  if(dump_all() && print_tests >= PRINT_ALL)
661  *out << "\nv_xI =\n" << *v_xI;
662  }
663  V_StMtV( v_chD_tmp.get(), -alpha, *N, no_trans, *v_xI );
664  V_InvMtV( v_xD_tmp.get(), *C, no_trans, *v_chD_tmp );
665  V_StMtV( v_xD.get(), alpha, *D, no_trans, *v_xI );
666  const value_type
667  sum_ICaNv = sum(*v_xD_tmp),
668  sum_aDv = sum(*v_xD);
669  assert_print_nan_inf(sum_ICaNv, "sum(inv(C)*(-alpha*N*v))",true,out);
670  assert_print_nan_inf(sum_aDv, "sum(alpha*D*v)",true,out);
671  const value_type
672  calc_err = ::fabs( ( sum_ICaNv - sum_aDv )
673  /( ::fabs(sum_ICaNv) + ::fabs(sum_aDv) + small_num ) );
674  if(out && print_tests >= PRINT_ALL)
675  *out
676  << "\nrel_err(sum(inv(C)*(-alpha*N*v)),sum(alpha*D*v)) = "
677  << "rel_err(" << sum_ICaNv << "," << sum_aDv << ") = "
678  << calc_err << std::endl;
679  if( calc_err >= warning_tol() ) {
680  if(out && print_tests >= PRINT_ALL)
681  *out
682  << std::endl
683  << ( calc_err >= error_tol() ? "Error" : "Warning" )
684  << ", rel_err(sum(inv(C)*(-alpha*N*v))),sum(alpha*D*v)) = "
685  << "rel_err(" << sum_ICaNv << "," << sum_aDv << ") = "
686  << calc_err
687  << " exceeded "
688  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
689  << " = "
690  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
691  << std::endl;
692  if(calc_err >= error_tol()) {
693  if(dump_all() && print_tests >= PRINT_ALL) {
694  *out << "\nalpha = " << alpha << std::endl;
695  *out << "\nv_xI =\n" << *v_xI;
696  *out << "\n-alpha*N*v_xI =\n" << *v_chD_tmp;
697  *out << "\ninv(C)*(-alpha*N*v_xI) =\n" << *v_xD_tmp;
698  *out << "\nalpha*D*v_xI =\n" << *v_xD;
699  }
700  llresult = false;
701  }
702  }
703  }}
704  if(!llresult) lresult = false;
705  if( out && print_tests == PRINT_MORE )
706  *out << " : " << ( llresult ? "passed" : "failed" )
707  << std::endl;
708 
709  if(out && print_tests >= PRINT_MORE)
710  *out << "\n3.b.1) Testing transposed: -alpha*N'*(inv(C')*v) == alpha*D'*v ...";
711  if(out && print_tests > PRINT_MORE)
712  *out << std::endl;
713  llresult = true;
714  {for( int k = 1; k <= num_random_tests(); ++k ) {
715  random_vector( rand_y_l, rand_y_u, v_xD.get() );
716  if(out && print_tests >= PRINT_ALL) {
717  *out
718  << "\n\n3.b.1."<<k<<") random vector " << k << " ( ||v_xD||_1 / n = " << (v_xD->norm_1() / v_xD->dim()) << " )\n";
719  if(dump_all() && print_tests >= PRINT_ALL)
720  *out << "\nv_xD =\n" << *v_xD;
721  }
722  V_InvMtV( v_chD_tmp.get(), *C, trans, *v_xD );
723  V_StMtV( v_xI_tmp.get(), -alpha, *N, trans, *v_chD_tmp );
724  V_StMtV( v_xI.get(), alpha, *D, trans, *v_xD );
725  const value_type
726  sum_aNTICTv = sum(*v_xI_tmp),
727  sum_aDTv = sum(*v_xI);
728  assert_print_nan_inf(sum_aNTICTv, "sum(-alpha*N'*(inv(C')*v))",true,out);
729  assert_print_nan_inf(sum_aDTv, "sum(alpha*D'*v)",true,out);
730  const value_type
731  calc_err = ::fabs( ( sum_aNTICTv - sum_aDTv )
732  /( ::fabs(sum_aNTICTv) + ::fabs(sum_aDTv) + small_num ) );
733  if(out && print_tests >= PRINT_ALL)
734  *out
735  << "\nrel_err(sum(-alpha*N'*(inv(C')*v)),sum(alpha*D'*v)) = "
736  << "rel_err(" << sum_aNTICTv << "," << sum_aDTv << ") = "
737  << calc_err << std::endl;
738  if( calc_err >= warning_tol() ) {
739  if(out && print_tests >= PRINT_ALL)
740  *out
741  << std::endl
742  << ( calc_err >= error_tol() ? "Error" : "Warning" )
743  << ", rel_err(sum(-alpha*N'*(inv(C')*v))),sum(alpha*D'*v)) = "
744  << "rel_err(" << sum_aNTICTv << "," << sum_aDTv << ") = "
745  << calc_err
746  << " exceeded "
747  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
748  << " = "
749  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
750  << std::endl;
751  if(calc_err >= error_tol()) {
752  if(dump_all() && print_tests >= PRINT_ALL) {
753  *out << "\nalpha = " << alpha << std::endl;
754  *out << "\nv_xD =\n" << *v_xD;
755  *out << "\ninv(C')**v_xD =\n" << *v_chD_tmp;
756  *out << "\n-alpha*N'*(inv(C')**v_xD) =\n" << *v_xI_tmp;
757  *out << "\nalpha*D*'v_xD =\n" << *v_xI;
758  }
759  llresult = false;
760  }
761  }
762  }}
763  if(!llresult) lresult = false;
764  if( out && print_tests == PRINT_MORE )
765  *out << " : " << ( llresult ? "passed" : "failed" )
766  << std::endl;
767 
768  }
769  }
770 
771  if( GcUP ) {
772  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Validate GcUP and the related matrices
773  }
774 
775  if(!lresult) success = false;
776  if( out && print_tests == PRINT_BASIC )
777  *out << " : " << ( lresult ? "passed" : "failed" );
778  }
779 
780  if(out && print_tests != PRINT_NONE ) {
781  if(success)
782  *out << "\nCongradulations! The BasisSystem object and its associated matrix objects seem to check out!\n";
783  else
784  *out << "\nOops! At last one of the tests did not check out!\n";
785  if(print_tests >= PRINT_BASIC)
786  *out << "\nEnd BasisSystemTester::test_basis_system(...)\n";
787  }
788 
789  return success;
790 }
791 
792 } // end namespace AbstractLinAlgPack
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
COOMatrixPartitionedViewUtilityPack::TransposedPartition< T_Indice, T_Value > trans(COOMatrixPartitionedViewUtilityPack::Partition< T_Indice, T_Value > &part)
Create a transposed view of a partition object.
Interface for the creation and maintainance of a basis matrix for a decomposition of linearlized cons...
void pow(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
vs_lhs = pow(vs_rhs1,vs_rhs2)
virtual Range1D equ_undecomp() const
Range of undecomposed general equality constriants.
void V_StMtV(SpVector *sv_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const DVectorSlice &vs_rhs2)
sv_lhs = alpha * op(P_rhs1) * vs_rhs2.
bool test_basis_system(const BasisSystem &basis_sys, const MatrixOp *Gc, const MatrixOpNonsing *C, const MatrixOp *N, const MatrixOp *D, const MatrixOp *GcUP, std::ostream *out)
Test a BasisSystem object after BasisSystem::update_basis() is called.
Index size() const
Return the size of the range (ubound() - lbound() + 1)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void add_matrix(size_type row_offset, size_type col_offset, value_type alpha, const GenPermMatrixSlice *P, const release_resource_ptr_t &P_release, BLAS_Cpp::Transp P_trans, const MatrixOp *A, const release_resource_ptr_t &A_release, BLAS_Cpp::Transp A_trans, const GenPermMatrixSlice *Q, const release_resource_ptr_t &Q_release, BLAS_Cpp::Transp Q_trans)
Add a sub-matrix alpha*op(P)*op(A)*op(Q).
Index ubound() const
Return upper bound of the range.
const MatNorm calc_norm(EMatNormType requested_norm_type=MAT_NORM_1, bool allow_replacement=false) const
Compute a norm of this matrix.
The induced infinity norm ||M||inf, i.e. max abs row sum.
T * get() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
void V_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = alpha * op(M_rhs1) * V_rhs2.
The print option has not been selected (will default to PRINT_NONE if not set)
Transposed.
virtual space_ptr_t sub_space(const Range1D &rng) const
Create a transient sub-space of the current vector space.
virtual void finish_construction(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows)
Call to finish the construction process.
virtual Range1D var_indep() const =0
Range of independnet (nonbasic) variables.
Not transposed.
virtual size_type cols() const
Return the number of columns in the matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
. One-based subregion index range class.
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
std::ostream * out
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
BasisSystemTester(EPrintTestLevel print_tests=PRINT_NOT_SELECTED, bool dump_all=false, bool throw_exception=true, size_type num_random_tests=1, value_type warning_tol=1e-14, value_type error_tol=1e-8)
Constructor (default options)
const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec LAPACK_C_Decl::f_dbl_prec * C
virtual bool is_compatible(const VectorSpace &vec_spc) const =0
Compare the compatibility of two vector spaces.
bool assert_print_nan_inf(const value_type &val, const char name[], bool throw_excpt, std::ostream *out)
This function asserts if a value_type scalare is a NaN or Inf and optionally prints out these entires...
Base class for all matrices that support basic matrix operations.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Index lbound() const
Return lower bound of the range.
virtual Range1D var_dep() const =0
Range of dependent (basic) variables.
virtual Range1D equ_decomp() const
Range of decomposed general equality constraints.
const MatNorm calc_cond_num(EMatNormType requested_norm_type=MAT_NORM_1, bool allow_replacement=false) const
Compute an estimate of the condition number of this matrix.
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
virtual const VectorSpace & space_cols() const =0
Vector space for vectors that are compatible with the columns of the matrix.
value_type sum(const Vector &v_rhs)
result = sum( v_rhs(i), i = 1,,,dim )
Abstract base class for all nonsingular polymorphic matrices that can be used to compute matrix-vecto...
const f_int const f_int & N
virtual vec_mut_ptr_t create_member() const =0
Create a vector member from the vector space.
virtual size_type rows() const
Return the number of rows in the matrix.
Print everything all the tests in great detail but output is independent of problem size...
Matrix class for matrices composed out of a set of other matrices and vectors.
void random_vector(value_type l, value_type u, VectorMutable *v)
Generate a random vector with elements uniformly distrubuted elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)
v_lhs += V_rhs.