AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_BasisSystemTester.cpp
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 
47 #include "AbstractLinAlgPack_BasisSystemTester.hpp"
48 #include "AbstractLinAlgPack_BasisSystem.hpp"
49 #include "AbstractLinAlgPack_VectorSpace.hpp"
50 #include "AbstractLinAlgPack_VectorMutable.hpp"
51 #include "AbstractLinAlgPack_VectorStdOps.hpp"
52 #include "AbstractLinAlgPack_VectorOut.hpp"
53 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
54 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
55 #include "AbstractLinAlgPack_MatrixComposite.hpp"
56 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
57 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
58 
59 namespace AbstractLinAlgPack {
60 
61 BasisSystemTester::BasisSystemTester(
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 
77 bool BasisSystemTester::test_basis_system(
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;
90  using AbstractLinAlgPack::sum;
91  using AbstractLinAlgPack::dot;
92  using AbstractLinAlgPack::Vp_StV;
93  using AbstractLinAlgPack::assert_print_nan_inf;
94  using AbstractLinAlgPack::random_vector;
95  using LinAlgOpPack::V_StMtV;
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 
110  EPrintTestLevel
111  print_tests = ( this->print_tests() == PRINT_NOT_SELECTED ? PRINT_NONE : this->print_tests() );
112 
113  MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
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
204  Teuchos::RCP<const AbstractLinAlgPack::MatrixOp>
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;
213  Teuchos::RCP<AbstractLinAlgPack::MatrixComposite>
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 
355  VectorSpace::vec_mut_ptr_t
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  {
514  VectorSpace::vec_mut_ptr_t
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  {
643  VectorSpace::vec_mut_ptr_t
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
AbstractLinAlgPack::MatrixBase::space_rows
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
AbstractLinAlgPack::BasisSystem
Interface for the creation and maintainance of a basis matrix for a decomposition of linearlized cons...
Definition: AbstractLinAlgPack_BasisSystem.hpp:228
AbstractLinAlgPack::BasisSystem::equ_undecomp
virtual Range1D equ_undecomp() const
Range of undecomposed general equality constriants.
Definition: AbstractLinAlgPack_BasisSystem.cpp:70
AbstractLinAlgPack::BasisSystemTester::test_basis_system
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.
Definition: AbstractLinAlgPack_BasisSystemTester.cpp:77
RangePack::Range1D::size
Index size() const
AbstractLinAlgPack::Vp_StV
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
Definition: AbstractLinAlgPack_VectorStdOps.cpp:191
AbstractLinAlgPack::MatrixComposite::add_matrix
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).
Definition: AbstractLinAlgPack_MatrixComposite.cpp:344
RangePack::Range1D::ubound
Index ubound() const
AbstractLinAlgPack::MatrixOp::MAT_NORM_INF
The induced infinity norm ||M||inf, i.e. max abs row sum.
Definition: AbstractLinAlgPack_MatrixOp.hpp:245
Teuchos::RCP::get
T * get() const
LinAlgOpPack::V_StV
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:450
AbstractLinAlgPack::BasisSystemTester::PRINT_NOT_SELECTED
The print option has not been selected (will default to PRINT_NONE if not set)
Definition: AbstractLinAlgPack_BasisSystemTester.hpp:77
trans
AbstractLinAlgPack::VectorSpace::sub_space
virtual space_ptr_t sub_space(const Range1D &rng) const
Create a transient sub-space of the current vector space.
Definition: AbstractLinAlgPack_VectorSpace.cpp:104
AbstractLinAlgPack::MatrixComposite::finish_construction
virtual void finish_construction(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows)
Call to finish the construction process.
Definition: AbstractLinAlgPack_MatrixComposite.cpp:554
AbstractLinAlgPack::BasisSystem::var_indep
virtual Range1D var_indep() const =0
Range of independnet (nonbasic) variables.
no_trans
AbstractLinAlgPack::MatrixBase::cols
virtual size_type cols() const
Return the number of columns in the matrix.
Definition: AbstractLinAlgPack_MatrixBase.cpp:52
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
LinAlgOpPack::V_VpV
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:458
AbstractLinAlgPack::V_InvMtV
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
Definition: AbstractLinAlgPack_MatrixNonsing.hpp:249
AbstractLinAlgPack::BasisSystemTester::BasisSystemTester
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)
Definition: AbstractLinAlgPack_BasisSystemTester.cpp:61
AbstractLinAlgPack::VectorSpace::is_compatible
virtual bool is_compatible(const VectorSpace &vec_spc) const =0
Compare the compatibility of two vector spaces.
AbstractLinAlgPack::MatrixOp
Base class for all matrices that support basic matrix operations.
Definition: AbstractLinAlgPack_MatrixOp.hpp:129
AbstractLinAlgPack::MatrixOp::calc_norm
const MatNorm calc_norm(EMatNormType requested_norm_type=MAT_NORM_1, bool allow_replacement=false) const
Compute a norm of this matrix.
Definition: AbstractLinAlgPack_MatrixOp.cpp:123
AbstractLinAlgPack::dot
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Definition: AbstractLinAlgPack_VectorStdOps.cpp:130
RangePack::Range1D::lbound
Index lbound() const
AbstractLinAlgPack::BasisSystem::var_dep
virtual Range1D var_dep() const =0
Range of dependent (basic) variables.
AbstractLinAlgPack::BasisSystem::equ_decomp
virtual Range1D equ_decomp() const
Range of decomposed general equality constraints.
Definition: AbstractLinAlgPack_BasisSystem.cpp:64
AbstractLinAlgPack::MatrixOpNonsing::calc_cond_num
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.
Definition: AbstractLinAlgPack_MatrixOpNonsing.cpp:64
AbstractLinAlgPack::MatrixBase::space_cols
virtual const VectorSpace & space_cols() const =0
Vector space for vectors that are compatible with the columns of the matrix.
AbstractLinAlgPack::sum
value_type sum(const Vector &v_rhs)
result = sum( v_rhs(i), i = 1,,,dim )
Definition: AbstractLinAlgPack_VectorStdOps.cpp:121
AbstractLinAlgPack::MatrixOpNonsing
Abstract base class for all nonsingular polymorphic matrices that can be used to compute matrix-vecto...
Definition: AbstractLinAlgPack_MatrixOpNonsing.hpp:54
AbstractLinAlgPack::VectorSpace::create_member
virtual vec_mut_ptr_t create_member() const =0
Create a vector member from the vector space.
AbstractLinAlgPack::MatrixBase::rows
virtual size_type rows() const
Return the number of rows in the matrix.
Definition: AbstractLinAlgPack_MatrixBase.cpp:47
AbstractLinAlgPack::BasisSystemTester::PRINT_ALL
Print everything all the tests in great detail but output is independent of problem size...
Definition: AbstractLinAlgPack_BasisSystemTester.hpp:81
AbstractLinAlgPack::MatrixComposite
Matrix class for matrices composed out of a set of other matrices and vectors.
Definition: AbstractLinAlgPack_MatrixComposite.hpp:126
AbstractLinAlgPack::random_vector
void random_vector(value_type l, value_type u, VectorMutable *v)
Generate a random vector with elements uniformly distrubuted elements.
Definition: AbstractLinAlgPack_VectorStdOps.cpp:247
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
LinAlgOpPack::Vp_V
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)
v_lhs += V_rhs.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:342
Generated on Wed Dec 2 2015 08:49:46 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by   doxygen 1.8.6