MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConstrainedOptPack_DecompositionSystemTester.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 <limits>
45 #include <ostream>
46 
47 #include "ConstrainedOptPack_DecompositionSystemTester.hpp"
48 #include "ConstrainedOptPack_DecompositionSystem.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_MatrixOpNonsingTester.hpp"
55 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
56 #include "AbstractLinAlgPack_MatrixComposite.hpp"
57 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
58 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
59 
60 namespace ConstrainedOptPack {
61 
62 DecompositionSystemTester::DecompositionSystemTester(
63  EPrintTestLevel print_tests
64  ,bool dump_all
65  ,bool throw_exception
66  ,size_type num_random_tests
67  ,value_type mult_warning_tol
68  ,value_type mult_error_tol
69  ,value_type solve_warning_tol
70  ,value_type solve_error_tol
71  )
72  :print_tests_(print_tests)
73  ,dump_all_(dump_all)
74  ,throw_exception_(throw_exception)
75  ,num_random_tests_(num_random_tests)
76  ,mult_warning_tol_(mult_warning_tol)
77  ,mult_error_tol_(mult_error_tol)
78  ,solve_warning_tol_(solve_warning_tol)
79  ,solve_error_tol_(solve_error_tol)
80 {}
81 
82 bool DecompositionSystemTester::test_decomp_system(
83  const DecompositionSystem &ds
84  ,const MatrixOp &Gc
85  ,const MatrixOp *Z
86  ,const MatrixOp *Y
87  ,const MatrixOpNonsing *R
88  ,const MatrixOp *Uz
89  ,const MatrixOp *Uy
90  ,std::ostream *out
91  )
92 {
93  namespace rcp = MemMngPack;
94  using BLAS_Cpp::no_trans;
95  using BLAS_Cpp::trans;
96  using AbstractLinAlgPack::sum;
97  using AbstractLinAlgPack::dot;
98  using AbstractLinAlgPack::Vp_StV;
99  using AbstractLinAlgPack::Vp_StMtV;
100  using AbstractLinAlgPack::assert_print_nan_inf;
101  using AbstractLinAlgPack::random_vector;
102  using LinAlgOpPack::V_StMtV;
103  using LinAlgOpPack::V_MtV;
104  using LinAlgOpPack::V_StV;
105  using LinAlgOpPack::V_VpV;
106  using LinAlgOpPack::Vp_V;
107 
108  bool success = true, result, lresult, llresult;
109  const value_type
110  rand_y_l = -1.0,
111  rand_y_u = 1.0,
112  small_num = ::pow(std::numeric_limits<value_type>::epsilon(),0.25),
113  alpha = 2.0,
114  beta = 3.0;
115 
116  EPrintTestLevel
117  print_tests = ( this->print_tests() == PRINT_NOT_SELECTED ? PRINT_NONE : this->print_tests() );
118 
119  // Print the input?
120  if( out && print_tests != PRINT_NONE ) {
121  if( print_tests >= PRINT_BASIC )
122  *out << "\n**********************************************************"
123  << "\n*** DecompositionSystemTester::test_decomp_system(...) ***"
124  << "\n**********************************************************\n";
125  }
126 
127  const size_type
128  n = ds.n(),
129  m = ds.m(),
130  r = ds.r();
131  const Range1D
132  equ_decomp = ds.equ_decomp(),
133  equ_undecomp = ds.equ_undecomp();
134 
135  // print dimensions, ranges
136  if( out && print_tests >= PRINT_MORE ) {
137  *out
138  << "\nds.n() = " << n
139  << "\nds.m() = " << m
140  << "\nds.r() = " << r
141  << "\nds.equ_decomp() = ["<<equ_decomp.lbound()<<","<<equ_decomp.ubound()<<"]"
142  << "\nds.equ_undecomp() = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
143  << "\nds.space_range()->dim() = " << ds.space_range()->dim()
144  << "\nds.space_null()->dim() = " << ds.space_null()->dim()
145  << std::endl;
146  }
147 
148  // validate input matrices
149  TEUCHOS_TEST_FOR_EXCEPTION(
150  Z==NULL&&Y==NULL&&R==NULL&&Uz==NULL&&Uy==NULL
151  , std::invalid_argument
152  ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
153  "at least one of Z, Y, R, Uz or Uy can not be NULL!" );
154  TEUCHOS_TEST_FOR_EXCEPTION(
155  m == r && Uz != NULL, std::invalid_argument
156  ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
157  "Uz must be NULL if m==r is NULL!" );
158  TEUCHOS_TEST_FOR_EXCEPTION(
159  m == r && Uy != NULL, std::invalid_argument
160  ,"DecompositionSystemTester::test_decomp_system(...) : Error, "
161  "Uy must be NULL if m==r is NULL!" );
162 
163  // Print the input?
164  if( out && print_tests != PRINT_NONE ) {
165  if(dump_all()) {
166  *out << "\nGc =\n" << Gc;
167  if(Z)
168  *out << "\nZ =\n" << *Z;
169  if(Y)
170  *out << "\nY =\n" << *Y;
171  if(R)
172  *out << "\nR =\n" << *R;
173  if(Uz)
174  *out << "\nUz =\n" << *Uz;
175  if(Uy)
176  *out << "\nUy =\n" << *Uy;
177  }
178  }
179 
180  //
181  // Check the dimensions of everything
182  //
183 
184  if( out && print_tests >= PRINT_BASIC )
185  *out << "\n1) Check the partitioning ranges and vector space dimensions ...";
186  lresult = true;
187 
188  if( out && print_tests >= PRINT_MORE )
189  *out << "\n\n1.a) check: equ_decomp.size() + equ_undecomp.size() == ds.m() : ";
190  result = equ_decomp.size() + equ_undecomp.size() == ds.m();
191  if(out && print_tests >= PRINT_MORE)
192  *out << ( result ? "passed" : "failed" );
193  if(!result) lresult = false;
194 
195  if( out && print_tests >= PRINT_MORE )
196  *out << "\n\n1.b) check: equ_decomp.size() == ds.r() : ";
197  result = equ_decomp.size() == ds.r();
198  if(out && print_tests >= PRINT_MORE)
199  *out << ( result ? "passed" : "failed" );
200  if(!result) lresult = false;
201 
202  if( out && print_tests >= PRINT_MORE )
203  *out << "\n\n1.c) check: ds.space_range()->dim() == ds.r() : ";
204  result = ds.space_range()->dim() == ds.r();
205  if(out && print_tests >= PRINT_MORE)
206  *out << ( result ? "passed" : "failed" );
207  if(!result) lresult = false;
208 
209  if( out && print_tests >= PRINT_MORE )
210  *out << "\n\n1.d) check: ds.space_null()->dim() == ds.n() - ds.r() : ";
211  result = ds.space_null()->dim() == ds.n() - ds.r();
212  if(out && print_tests >= PRINT_MORE)
213  *out << ( result ? "passed" : "failed" );
214  if(!result) lresult = false;
215 
216  if(out && print_tests >= PRINT_MORE)
217  *out << std::endl;
218 
219  if(!lresult) success = false;
220  if( out && print_tests == PRINT_BASIC )
221  *out << " : " << ( lresult ? "passed" : "failed" );
222 
223  //
224  // Perform the tests
225  //
226 
227  if(out && print_tests >= PRINT_BASIC)
228  *out
229  << "\n2) Check the compatibility of the vector spaces for Gc, Z, Y, R, Uz and Uy ...";
230  lresult = true;
231 
232  if(Z) {
233  if(out && print_tests >= PRINT_MORE)
234  *out
235  << "\n2.a) Check consistency of the vector spaces for:"
236  << "\n Z.space_cols() == Gc.space_cols() and Z.space_rows() == ds.space_null()";
237  llresult = true;
238  if(out && print_tests >= PRINT_ALL)
239  *out << "\n\n2.a.1) Z->space_cols().is_compatible(Gc.space_cols()) == true : ";
240  result = Z->space_cols().is_compatible(Gc.space_cols());
241  if(out && print_tests >= PRINT_ALL)
242  *out << ( result ? "passed" : "failed" )
243  << std::endl;
244  if(!result) llresult = false;
245  if(out && print_tests >= PRINT_ALL)
246  *out << "\n\n2.a.2) Z->space_cols().is_compatible(*ds.space_null()) == true : ";
247  result = Z->space_rows().is_compatible(*ds.space_null());
248  if(out && print_tests >= PRINT_ALL)
249  *out << ( result ? "passed" : "failed" )
250  << std::endl;
251  if(!result) llresult = false;
252  if(!llresult) lresult = false;
253  if( out && print_tests == PRINT_MORE )
254  *out << " : " << ( llresult ? "passed" : "failed" );
255  }
256 
257  if(Y) {
258  if(out && print_tests >= PRINT_MORE)
259  *out
260  << "\n2.b) Check consistency of the vector spaces for:"
261  << "\n Y.space_cols() == Gc.space_cols() and Y.space_rows() == ds.space_range()";
262  llresult = true;
263  if(out && print_tests >= PRINT_ALL)
264  *out << "\n\n2.b.1) Y->space_cols().is_compatible(Gc.space_cols()) == true : ";
265  result = Y->space_cols().is_compatible(Gc.space_cols());
266  if(out && print_tests >= PRINT_ALL)
267  *out << ( result ? "passed" : "failed" )
268  << std::endl;
269  if(!result) llresult = false;
270  if(out && print_tests >= PRINT_ALL)
271  *out << "\n\n2.b.2) Y->space_cols().is_compatible(*ds.space_range()) == true : ";
272  result = Y->space_rows().is_compatible(*ds.space_range());
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(R) {
283  if(out && print_tests >= PRINT_MORE)
284  *out
285  << "\n2.c) Check consistency of the vector spaces for:"
286  << "\n R.space_cols() == Gc.space_cols()(equ_decomp) and R.space_rows() == ds.space_range()";
287  llresult = true;
288  if(out && print_tests >= PRINT_ALL)
289  *out << "\n\n2.c.1) R->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_decomp)) == true : ";
290  result = R->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_decomp));
291  if(out && print_tests >= PRINT_ALL)
292  *out << ( result ? "passed" : "failed" )
293  << std::endl;
294  if(!result) llresult = false;
295  if(out && print_tests >= PRINT_ALL)
296  *out << "\n\n2.c.2) R->space_cols().is_compatible(*ds.space_range()) == true : ";
297  result = R->space_rows().is_compatible(*ds.space_range());
298  if(out && print_tests >= PRINT_ALL)
299  *out << ( result ? "passed" : "failed" )
300  << std::endl;
301  if(!result) llresult = false;
302  if(!llresult) lresult = false;
303  if( out && print_tests == PRINT_MORE )
304  *out << " : " << ( llresult ? "passed" : "failed" );
305  }
306 
307  if(Uz) {
308  if(out && print_tests >= PRINT_MORE)
309  *out
310  << "\n2.d) Check consistency of the vector spaces for:"
311  << "\n Uz.space_cols() == Gc.space_cols()(equ_undecomp) and Uz.space_rows() == ds.space_null()";
312  llresult = true;
313  if(out && print_tests >= PRINT_ALL)
314  *out << "\n\n2.d.1) Uz->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp)) == true : ";
315  result = Uz->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp));
316  if(out && print_tests >= PRINT_ALL)
317  *out << ( result ? "passed" : "failed" )
318  << std::endl;
319  if(!result) llresult = false;
320  if(out && print_tests >= PRINT_ALL)
321  *out << "\n\n2.d.2) Uz->space_cols().is_compatible(*ds.space_null()) == true : ";
322  result = Uz->space_rows().is_compatible(*ds.space_null());
323  if(out && print_tests >= PRINT_ALL)
324  *out << ( result ? "passed" : "failed" )
325  << std::endl;
326  if(!result) llresult = false;
327  if(!llresult) lresult = false;
328  if( out && print_tests == PRINT_MORE )
329  *out << " : " << ( llresult ? "passed" : "failed" );
330  }
331 
332  if(Uy) {
333  if(out && print_tests >= PRINT_MORE)
334  *out
335  << "\n2.e) Check consistency of the vector spaces for:"
336  << "\n Uy.space_cols() == Gc.space_cols()(equ_undecomp) and Uy.space_rows() == ds.space_range()";
337  llresult = true;
338  if(out && print_tests >= PRINT_ALL)
339  *out << "\n\n2.e.1) Uy->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp)) == true : ";
340  result = Uy->space_cols().is_compatible(*Gc.space_cols().sub_space(equ_undecomp));
341  if(out && print_tests >= PRINT_ALL)
342  *out << ( result ? "passed" : "failed" )
343  << std::endl;
344  if(!result) llresult = false;
345  if(out && print_tests >= PRINT_ALL)
346  *out << "\n\n2.e.2) Uy->space_cols().is_compatible(*ds.space_range()) == true : ";
347  result = Uy->space_rows().is_compatible(*ds.space_range());
348  if(out && print_tests >= PRINT_ALL)
349  *out << ( result ? "passed" : "failed" )
350  << std::endl;
351  if(!result) llresult = false;
352  if(!llresult) lresult = false;
353  if( out && print_tests == PRINT_MORE )
354  *out << " : " << ( llresult ? "passed" : "failed" );
355  }
356 
357  if(!lresult) success = false;
358  if( out && print_tests == PRINT_BASIC )
359  *out << " : " << ( lresult ? "passed" : "failed" );
360 
361  if(out && print_tests >= PRINT_BASIC)
362  *out
363  << "\n3) Check the compatibility of the matrices Gc, Z, Y, R, Uz and Uy numerically ...";
364 
365  if(Z) {
366 
367  if(out && print_tests >= PRINT_MORE)
368  *out
369  << std::endl
370  << "\n3.a) Check consistency of:"
371  << "\n op ( alpha* Gc(:,equ_decomp)' * beta*Z ) * v"
372  << "\n \\__________________________________/"
373  << "\n A"
374  << "\n == op( alpha*beta*Uz * v"
375  << "\n \\___________/"
376  << "\n B"
377  << "\nfor random vectors v ...";
378 
379  VectorSpace::vec_mut_ptr_t
380  v_c = Gc.space_rows().create_member(),
381  v_c_tmp = v_c->space().create_member(),
382  v_x = Gc.space_cols().create_member(),
383  v_x_tmp = v_x->space().create_member(),
384  v_z = ds.space_null()->create_member(),
385  v_z_tmp = v_z->space().create_member();
386 
387  if(out && print_tests >= PRINT_MORE)
388  *out << "\n\n3.a.1) Testing non-transposed A*v == B*v ...";
389  if(out && print_tests > PRINT_MORE)
390  *out << std::endl;
391  llresult = true;
392  {for( int k = 1; k <= num_random_tests(); ++k ) {
393  random_vector( rand_y_l, rand_y_u, v_z.get() );
394  if(out && print_tests >= PRINT_ALL) {
395  *out
396  << "\n3.a.1."<<k<<") random vector " << k << " ( ||v_z||_1 / n = " << (v_z->norm_1() / v_z->dim()) << " )\n";
397  if(dump_all() && print_tests >= PRINT_ALL)
398  *out << "\nv_z =\n" << *v_z;
399  }
400  V_StMtV( v_x.get(), beta, *Z, no_trans, *v_z );
401  V_StMtV( v_c.get(), alpha, Gc, trans, *v_x );
402  *v_c_tmp->sub_view(equ_decomp) = 0.0;
403  if(equ_undecomp.size()) {
404  if(Uz)
405  V_StMtV( v_c_tmp->sub_view(equ_undecomp).get(), alpha*beta, *Uz, no_trans, *v_z );
406  else
407  *v_c_tmp->sub_view(equ_undecomp).get() = *v_c->sub_view(equ_undecomp);
408  }
409  const value_type
410  sum_Bv = sum(*v_c_tmp), // should be zero if equ_undecomp.size() == 0 so scale by 1.0
411  sum_Av = sum(*v_c);
412  assert_print_nan_inf(sum_Bv, "sum(B*v_z)",true,out);
413  assert_print_nan_inf(sum_Av, "sum(A*v_z)",true,out);
414  const value_type
415  calc_err = ::fabs( ( sum_Av - sum_Bv )
416  /( ::fabs(sum_Av) + ::fabs(sum_Bv) + (equ_undecomp.size() ? small_num : 1.0) ) );
417  if(out && print_tests >= PRINT_ALL)
418  *out
419  << "\nrel_err(sum(A*v_z),sum(B*v_z)) = "
420  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
421  << calc_err << std::endl;
422  if( calc_err >= mult_warning_tol() ) {
423  if(out && print_tests >= PRINT_ALL)
424  *out
425  << std::endl
426  << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
427  << ", rel_err(sum(A*v_z),sum(B*v_z)) = "
428  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
429  << calc_err
430  << " exceeded "
431  << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
432  << " = "
433  << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
434  << std::endl;
435  if(calc_err >= mult_error_tol()) {
436  if(dump_all() && print_tests >= PRINT_ALL) {
437  *out << "\nalpha = " << alpha << std::endl;
438  *out << "\nbeta = " << beta << std::endl;
439  *out << "\nv_z =\n" << *v_z;
440  *out << "\nbeta*Z*v_z =\n" << *v_x;
441  *out << "\nA*v_z =\n" << *v_c;
442  *out << "\nB*v_z =\n" << *v_c_tmp;
443  }
444  llresult = false;
445  }
446  }
447  }}
448  if(!llresult) lresult = false;
449  if( out && print_tests == PRINT_MORE )
450  *out << " : " << ( llresult ? "passed" : "failed" )
451  << std::endl;
452 
453  if(out && print_tests >= PRINT_MORE)
454  *out << "\n\n3.a.2) Testing transposed A'*v == B'*v ...";
455  if(out && print_tests > PRINT_MORE)
456  *out << std::endl;
457  llresult = true;
458  {for( int k = 1; k <= num_random_tests(); ++k ) {
459  random_vector( rand_y_l, rand_y_u, v_c.get() );
460  if(out && print_tests >= PRINT_ALL) {
461  *out
462  << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_c||_1 / n = " << (v_c->norm_1() / v_c->dim()) << " )\n";
463  if(dump_all() && print_tests >= PRINT_ALL)
464  *out << "\nv_c =\n" << *v_c;
465  }
466  V_StMtV( v_x.get(), alpha, Gc, no_trans, *v_c );
467  V_StMtV( v_z.get(), beta, *Z, trans, *v_x );
468  *v_z_tmp = 0.0;
469  if(equ_undecomp.size()) {
470  if(Uz)
471  V_StMtV( v_z_tmp.get(), alpha*beta, *Uz, trans, *v_c->sub_view(equ_undecomp) );
472  else
473  *v_z_tmp = *v_z;
474  }
475  const value_type
476  sum_Bv = sum(*v_z_tmp), // should be zero so scale by 1.0
477  sum_Av = sum(*v_z);
478  assert_print_nan_inf(sum_Bv, "sum(B'*v_c)",true,out);
479  assert_print_nan_inf(sum_Av, "sum(A'*v_c)",true,out);
480  const value_type
481  calc_err = ::fabs( ( sum_Av - sum_Bv )
482  /( ::fabs(sum_Av) + ::fabs(sum_Bv) + (equ_undecomp.size() ? small_num : 1.0) ) );
483  if(out && print_tests >= PRINT_ALL)
484  *out
485  << "\nrel_err(sum(A'*v_c),sum(B'*v_c)) = "
486  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
487  << calc_err << std::endl;
488  if( calc_err >= mult_warning_tol() ) {
489  if(out && print_tests >= PRINT_ALL)
490  *out
491  << std::endl
492  << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
493  << ", rel_err(sum(A'*v_c),sum(B'*v_c)) = "
494  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
495  << calc_err
496  << " exceeded "
497  << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
498  << " = "
499  << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
500  << std::endl;
501  if(calc_err >= mult_error_tol()) {
502  if(dump_all() && print_tests >= PRINT_ALL) {
503  *out << "\nalpha = " << alpha << std::endl;
504  *out << "\nbeta = " << beta << std::endl;
505  *out << "\nv_c =\n" << *v_c;
506  *out << "\nalpha*Gc*v_c =\n" << *v_x;
507  *out << "\nA'*v_c =\n" << *v_z;
508  *out << "\nB'*v_c =\n" << *v_z_tmp;
509  }
510  llresult = false;
511  }
512  }
513  }}
514  if(!llresult) lresult = false;
515  if( out && print_tests == PRINT_MORE )
516  *out << " : " << ( llresult ? "passed" : "failed" )
517  << std::endl;
518 
519  }
520  else {
521  if(out && print_tests >= PRINT_MORE)
522  *out
523  << std::endl
524  << "\n3.a) Warning! Z ==NULL; Z, and Uz are not checked numerically ...\n";
525  }
526 
527  if(Y) {
528 
529  if(out && print_tests >= PRINT_MORE)
530  *out
531  << std::endl
532  << "\n3.b) Check consistency of:"
533  << "\n op ( alpha*[ Gc(:,equ_decomp)' ]"
534  << "\n [ Gc(:,equ_undecomp)' ] * beta*Y ) * v"
535  << "\n \\_____________________________________/"
536  << "\n A"
537  << "\n == op( alpha*beta*[ R ]"
538  << "\n [ Uy ] ) * v"
539  << "\n \\_______________/"
540  << "\n B"
541  << "\nfor random vectors v ...";
542 
543  VectorSpace::vec_mut_ptr_t
544  v_c = Gc.space_rows().create_member(),
545  v_c_tmp = v_c->space().create_member(),
546  v_x = Gc.space_cols().create_member(),
547  v_x_tmp = v_x->space().create_member(),
548  v_y = ds.space_range()->create_member(),
549  v_y_tmp = v_y->space().create_member();
550 
551  if(out && print_tests >= PRINT_MORE)
552  *out << "\n\n3.b.1) Testing non-transposed A*v == B*v ...";
553  if(out && print_tests > PRINT_MORE)
554  *out << std::endl;
555  llresult = true;
556  {for( int k = 1; k <= num_random_tests(); ++k ) {
557  random_vector( rand_y_l, rand_y_u, v_y.get() );
558  if(out && print_tests >= PRINT_ALL) {
559  *out
560  << "\n3.b.1."<<k<<") random vector " << k << " ( ||v_y||_1 / n = " << (v_y->norm_1() / v_y->dim()) << " )\n";
561  if(dump_all() && print_tests >= PRINT_ALL)
562  *out << "\nv_y =\n" << *v_y;
563  }
564  V_StMtV( v_x.get(), beta, *Y, no_trans, *v_y );
565  V_StMtV( v_c.get(), alpha, Gc, trans, *v_x );
566  V_StMtV( v_c_tmp->sub_view(equ_decomp).get(), alpha*beta, *R, no_trans, *v_y );
567  if(equ_undecomp.size()) {
568  if(Uy)
569  V_StMtV( v_c_tmp->sub_view(equ_undecomp).get(), alpha*beta, *Uy, no_trans, *v_y );
570  else
571  *v_c_tmp->sub_view(equ_undecomp) = *v_c->sub_view(equ_undecomp);
572  }
573  const value_type
574  sum_Bv = sum(*v_c_tmp),
575  sum_Av = sum(*v_c);
576  assert_print_nan_inf(sum_Bv, "sum(B*v_y)",true,out);
577  assert_print_nan_inf(sum_Av, "sum(A*v_y)",true,out);
578  const value_type
579  calc_err = ::fabs( ( sum_Av - sum_Bv )
580  /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
581  if(out && print_tests >= PRINT_ALL)
582  *out
583  << "\nrel_err(sum(A*v_y),sum(B*v_y)) = "
584  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
585  << calc_err << std::endl;
586  if( calc_err >= mult_warning_tol() ) {
587  if(out && print_tests >= PRINT_ALL)
588  *out
589  << std::endl
590  << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
591  << ", rel_err(sum(A*v_y),sum(B*v_y)) = "
592  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
593  << calc_err
594  << " exceeded "
595  << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
596  << " = "
597  << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
598  << std::endl;
599  if(calc_err >= mult_error_tol()) {
600  if(dump_all() && print_tests >= PRINT_ALL) {
601  *out << "\nalpha = " << alpha << std::endl;
602  *out << "\nbeta = " << beta << std::endl;
603  *out << "\nv_y =\n" << *v_y;
604  *out << "\nbeta*Y*v_y =\n" << *v_x;
605  *out << "\nA*v_y =\n" << *v_c;
606  *out << "\nB*v_y =\n" << *v_c_tmp;
607  }
608  llresult = false;
609  }
610  }
611  }}
612  if(!llresult) lresult = false;
613  if( out && print_tests == PRINT_MORE )
614  *out << " : " << ( llresult ? "passed" : "failed" )
615  << std::endl;
616 
617  if(out && print_tests >= PRINT_MORE)
618  *out << "\n\n3.b.2) Testing transposed A'*v == B'*v ...";
619  if(out && print_tests > PRINT_MORE)
620  *out << std::endl;
621  llresult = true;
622  {for( int k = 1; k <= num_random_tests(); ++k ) {
623  random_vector( rand_y_l, rand_y_u, v_c.get() );
624  if(out && print_tests >= PRINT_ALL) {
625  *out
626  << "\n3.a.2."<<k<<") random vector " << k << " ( ||v_c||_1 / n = " << (v_c->norm_1() / v_c->dim()) << " )\n";
627  if(dump_all() && print_tests >= PRINT_ALL)
628  *out << "\nv_c =\n" << *v_c;
629  }
630  V_StMtV( v_x.get(), alpha, Gc, no_trans, *v_c );
631  V_StMtV( v_y.get(), beta, *Y, trans, *v_x );
632  V_StMtV( v_y_tmp.get(), alpha*beta, *R, trans, *v_c->sub_view(equ_decomp) );
633  if(equ_undecomp.size()) {
634  if(Uy)
635  Vp_StMtV( v_y_tmp.get(), alpha*beta, *Uy, trans, *v_c->sub_view(equ_undecomp) );
636  else
637  Vp_V( v_y_tmp.get(), *v_y );
638  }
639  const value_type
640  sum_Bv = sum(*v_y_tmp), // should be zero so scale by 1.0
641  sum_Av = sum(*v_y);
642  assert_print_nan_inf(sum_Bv, "sum(B'*v_c)",true,out);
643  assert_print_nan_inf(sum_Av, "sum(A'*v_c)",true,out);
644  const value_type
645  calc_err = ::fabs( ( sum_Av - sum_Bv )
646  /( ::fabs(sum_Av) + ::fabs(sum_Bv) + small_num ) );
647  if(out && print_tests >= PRINT_ALL)
648  *out
649  << "\nrel_err(sum(A'*v_c),sum(B'*v_c)) = "
650  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
651  << calc_err << std::endl;
652  if( calc_err >= mult_warning_tol() ) {
653  if(out && print_tests >= PRINT_ALL)
654  *out
655  << std::endl
656  << ( calc_err >= mult_error_tol() ? "Error" : "Warning" )
657  << ", rel_err(sum(A'*v_c),sum(B'*v_c)) = "
658  << "rel_err(" << sum_Av << "," << sum_Bv << ") = "
659  << calc_err
660  << " exceeded "
661  << ( calc_err >= mult_error_tol() ? "mult_error_tol" : "mult_warning_tol" )
662  << " = "
663  << ( calc_err >= mult_error_tol() ? mult_error_tol() : mult_warning_tol() )
664  << std::endl;
665  if(calc_err >= mult_error_tol()) {
666  if(dump_all() && print_tests >= PRINT_ALL) {
667  *out << "\nalpha = " << alpha << std::endl;
668  *out << "\nbeta = " << beta << std::endl;
669  *out << "\nv_c =\n" << *v_c;
670  *out << "\nalpha*Gc*v_c =\n" << *v_x;
671  *out << "\nA'*v_c =\n" << *v_y;
672  *out << "\nB'*v_c =\n" << *v_y_tmp;
673  }
674  llresult = false;
675  }
676  }
677  }}
678  if(!llresult) lresult = false;
679  if( out && print_tests == PRINT_MORE )
680  *out << " : " << ( llresult ? "passed" : "failed" )
681  << std::endl;
682 
683  }
684  else {
685  if(out && print_tests >= PRINT_MORE)
686  *out
687  << std::endl
688  << "\n3.b) Warning! Y ==NULL; Y, R and Uy are not checked numerically ...\n";
689  }
690 
691  if(R) {
692  if(out && print_tests >= PRINT_MORE)
693  *out
694  << std::endl
695  << "\n3.b) Check consistency of: op(op(inv(R))*op(R)) == I ...\n";
696  typedef MatrixOpNonsingTester MWONST_t;
697  MWONST_t::EPrintTestLevel
698  olevel;
699  switch(print_tests) {
700  case PRINT_NONE:
701  case PRINT_BASIC:
702  olevel = MWONST_t::PRINT_NONE;
703  break;
704  case PRINT_MORE:
705  olevel = MWONST_t::PRINT_MORE;
706  break;
707  case PRINT_ALL:
708  olevel = MWONST_t::PRINT_ALL;
709  break;
710  default:
711  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here
712  }
713  MWONST_t
714  R_tester(
715  MWONST_t::TEST_LEVEL_2_BLAS
716  ,olevel
717  ,dump_all()
718  ,throw_exception()
719  ,num_random_tests()
720  ,solve_warning_tol()
721  ,solve_error_tol()
722  );
723  lresult = R_tester.test_matrix(*R,"R",out);
724  }
725 
726  if(!lresult) success = false;
727  if( out && print_tests == PRINT_BASIC )
728  *out << " : " << ( lresult ? "passed" : "failed" );
729 
730  if( out && print_tests != PRINT_NONE ) {
731  if(success)
732  *out << "\nCongradulations! The DecompositionSystem object and its associated matrix objects seem to check out!\n";
733  else
734  *out << "\nOops! At least one of the tests did not check out!\n";
735  if( print_tests >= PRINT_BASIC )
736  *out << "\nEnd DecompositionSystemTester::test_decomp_system(...)\n";
737  }
738 
739  return success;
740 }
741 
742 } // end namespace ConstrainedOptPack
BLAS_C_Decl::Y
f_dbl_prec const f_int f_dbl_prec * Y
Definition: DenseLinAlgPack_BLAS_Cpp.cpp:64
LinAlgOpPack::size_type
AbstractLinAlgPack::size_type size_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:54
ConstrainedOptPack::DecompositionSystem::equ_undecomp
virtual Range1D equ_undecomp() const
Returns the range of the undecomposed equalities.
Definition: ConstrainedOptPack_DecompositionSystem.cpp:63
DenseLinAlgPack::pow
void pow(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
vs_lhs = pow(vs_rhs1,vs_rhs2)
Definition: DenseLinAlgPack_DVectorOpBLAS.cpp:250
ConstrainedOptPack_DecompositionSystemTester.hpp
ConstrainedOptPack_DecompositionSystem.hpp
AbstractLinAlgPack_MatrixComposite.hpp
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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
ConstrainedOptPack::DecompositionSystemTester::test_decomp_system
bool test_decomp_system(const DecompositionSystem &decomp_sys, const MatrixOp &Gc, const MatrixOp *Z, const MatrixOp *Y, const MatrixOpNonsing *R, const MatrixOp *Uz, const MatrixOp *Uy, std::ostream *out)
Test a DecompositionSystem object after DecompositionSystem::update_basis() is called.
Definition: ConstrainedOptPack_DecompositionSystemTester.cpp:82
ConstrainedOptPack::DecompositionSystem::r
virtual size_type r() const
Returns the rank of Gc(:,equ_decomp()).
Definition: ConstrainedOptPack_DecompositionSystem.cpp:51
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
ConstrainedOptPack::DecompositionSystemTester::DecompositionSystemTester
DecompositionSystemTester(EPrintTestLevel print_tests=PRINT_NOT_SELECTED, bool dump_all=false, bool throw_exception=true, size_type num_random_tests=1, value_type mult_warning_tol=1e-14, value_type mult_error_tol=1e-8, value_type solve_warning_tol=1e-14, value_type solve_error_tol=1e-8)
Constructor (default options)
Definition: ConstrainedOptPack_DecompositionSystemTester.cpp:62
LinAlgOpPack::V_StMtV
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.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:491
ConstrainedOptPack::DecompositionSystemTester::PRINT_ALL
Print everything all the tests in great detail but output is independent of problem size...
Definition: ConstrainedOptPack_DecompositionSystemTester.hpp:81
ConstrainedOptPack::DecompositionSystem::space_range
virtual const VectorSpace::space_ptr_t space_range() const =0
Return a VectorSpace object for the range space.
BLAS_Cpp::trans
Transposed.
Definition: BLAS_Cpp_Types.hpp:75
ConstrainedOptPack::DecompositionSystem::equ_decomp
virtual Range1D equ_decomp() const
Returns the range of the decomposed equalities.
Definition: ConstrainedOptPack_DecompositionSystem.cpp:56
BLAS_Cpp::no_trans
Not transposed.
Definition: BLAS_Cpp_Types.hpp:74
AbstractLinAlgPack_assert_print_nan_inf.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
AbstractLinAlgPack_LinAlgOpPack.hpp
AbstractLinAlgPack_VectorMutable.hpp
LinAlgOpPack::V_VpV
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:458
QPKWIK_Output::out
std::ostream * out
AbstractLinAlgPack_VectorStdOps.hpp
ConstrainedOptPack::DecompositionSystem
This class abstracts a decomposition choice for the quasi-range space Y and null space Z matrices for...
Definition: ConstrainedOptPack_DecompositionSystem.hpp:89
AbstractLinAlgPack::Vp_StMtV
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
Definition: AbstractLinAlgPack_MatrixOp.hpp:857
AbstractLinAlgPack::assert_print_nan_inf
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...
Definition: AbstractLinAlgPack_assert_print_nan_inf.cpp:70
AbstractLinAlgPack::dot
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Definition: AbstractLinAlgPack_VectorStdOps.cpp:130
ConstrainedOptPack::DecompositionSystem::n
virtual size_type n() const
Return the number of rows in Gc.
Definition: ConstrainedOptPack_DecompositionSystem.cpp:46
LinAlgOpPack::V_MtV
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.
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:500
LinAlgOpPack::value_type
AbstractLinAlgPack::value_type value_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:56
AbstractLinAlgPack_VectorSpace.hpp
ConstrainedOptPack::DecompositionSystem::space_null
virtual const VectorSpace::space_ptr_t space_null() const =0
Return a VectorSpace object for the range space.
AbstractLinAlgPack::sum
value_type sum(const Vector &v_rhs)
result = sum( v_rhs(i), i = 1,,,dim )
Definition: AbstractLinAlgPack_VectorStdOps.cpp:121
ConstrainedOptPack::DecompositionSystemTester::PRINT_NOT_SELECTED
The print option has not been selected (will default to PRINT_NONE if not set)
Definition: ConstrainedOptPack_DecompositionSystemTester.hpp:77
ConstrainedOptPack::DecompositionSystem::m
virtual size_type m() const =0
Return the number of columns in Gc.
AbstractLinAlgPack_MatrixOpNonsingTester.hpp
AbstractLinAlgPack_MatrixOpNonsing.hpp
n
int n
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)
AbstractLinAlgPack_VectorOut.hpp
AbstractLinAlgPack_MatrixOpOut.hpp
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:50:24 for MOOCHO (Single Doxygen Collection) by   doxygen 1.8.6