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_MatrixOpNonsingTester.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_MatrixOpNonsingTester.hpp"
48 #include "AbstractLinAlgPack_MatrixOpNonsing.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 
62  ETestLevel test_level
63  ,EPrintTestLevel print_tests
64  ,bool dump_all
65  ,bool throw_exception
66  ,size_type num_random_tests
67  ,value_type warning_tol
68  ,value_type error_tol
69  )
70  :test_level_(test_level)
71  ,print_tests_(print_tests)
72  ,dump_all_(dump_all)
73  ,throw_exception_(throw_exception)
74  ,num_random_tests_(num_random_tests)
75  ,warning_tol_(warning_tol)
76  ,error_tol_(error_tol)
77 {}
78 
80  const MatrixOpNonsing &M
81  ,const char M_name[]
82  ,std::ostream *out
83  )
84 {
85  namespace rcp = MemMngPack;
86  using BLAS_Cpp::no_trans;
87  using BLAS_Cpp::trans;
88  using BLAS_Cpp::left;
89  using BLAS_Cpp::right;
93  using AbstractLinAlgPack::assert_print_nan_inf;
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;
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  //
111  // Perform the tests
112  //
113 
114  if(out && print_tests() >= PRINT_BASIC)
115  *out
116  << "\nCheck: alpha*op(op(inv("<<M_name<<"))*op("<<M_name<<"))*v == alpha*v ...";
117  if(out && print_tests() > PRINT_BASIC)
118  *out << std::endl;
119 
121  v_c1 = M.space_cols().create_member(),
122  v_c2 = M.space_cols().create_member(),
123  v_r1 = M.space_rows().create_member(),
124  v_r2 = M.space_rows().create_member();
125 
126  // Side of the matrix inverse
127  const BLAS_Cpp::Side a_side[2] = { BLAS_Cpp::left, BLAS_Cpp::right };
128  // If the matrices are transposed or not
129  const BLAS_Cpp::Transp a_trans[2] = { BLAS_Cpp::no_trans, BLAS_Cpp::trans };
130 
131  for( int side_i = 0; side_i < 2; ++side_i ) {
132  for( int trans_i = 0; trans_i < 2; ++trans_i ) {
133  const BLAS_Cpp::Side t_side = a_side[side_i];
134  const BLAS_Cpp::Transp t_trans = a_trans[trans_i];
135  if(out && print_tests() >= PRINT_MORE)
136  *out
137  << "\n" << side_i+1 << "." << trans_i+1 << ") "
138  << "Check: (t2 = "<<(t_side==left?"inv(":"alpha * ")<< M_name<<(t_trans==trans?"\'":"")<<(t_side==left?")":"")
139  << " * (t1 = "<<(t_side==right?"inv(":"alpha * ")<<M_name<<(t_trans==trans?"\'":"")<<(t_side==right?")":"")
140  << " * v)) == alpha * v ...";
141  if(out && print_tests() > PRINT_MORE)
142  *out << std::endl;
143  result = true;
145  *v = NULL,
146  *t1 = NULL,
147  *t2 = NULL;
148  if( (t_side == left && t_trans == no_trans) || (t_side == right && t_trans == trans) ) {
149  // (inv(R)*R*v or R'*inv(R')*v
150  v = v_r1.get();
151  t1 = v_c1.get();
152  t2 = v_r2.get();
153  }
154  else {
155  // (inv(R')*R'*v or R*inv(R)*v
156  v = v_c1.get();
157  t1 = v_r1.get();
158  t2 = v_c2.get();
159  }
160  for( int k = 1; k <= num_random_tests(); ++k ) {
161  lresult = true;
162  random_vector( rand_y_l, rand_y_u, v );
163  if(out && print_tests() >= PRINT_ALL) {
164  *out
165  << "\n"<<side_i+1<<"."<<trans_i+1<<"."<<k<<") random vector " << k
166  << " ( ||v||_1 / n = " << (v->norm_1() / v->dim()) << " )\n";
167  if(dump_all() && print_tests() >= PRINT_ALL)
168  *out << "\nv =\n" << *v;
169  }
170  // t1
171  if( t_side == right ) {
172  // t1 = inv(op(M))*v
173  V_InvMtV( t1, M, t_trans, *v );
174  }
175  else {
176  // t1 = alpha*op(M)*v
177  V_StMtV( t1, alpha, M, t_trans, *v );
178  }
179  // t2
180  if( t_side == left ) {
181  // t2 = inv(op(M))*t1
182  V_InvMtV( t2, M, t_trans, *t1 );
183  }
184  else {
185  // t2 = alpha*op(M)*t1
186  V_StMtV( t2, alpha, M, t_trans, *t1 );
187  }
188  const value_type
189  sum_t2 = sum(*t2),
190  sum_av = alpha*sum(*v);
191  assert_print_nan_inf(sum_t2, "sum(t2)",true,out);
192  assert_print_nan_inf(sum_av, "sum(alpha*t1)",true,out);
193  const value_type
194  calc_err = ::fabs( ( sum_av - sum_t2 )
195  /( ::fabs(sum_av) + ::fabs(sum_t2) + small_num ) );
196  if(out && print_tests() >= PRINT_ALL)
197  *out
198  << "\nrel_err(sum(alpha*v),sum(t2)) = "
199  << "rel_err(" << sum_av << "," << sum_t2 << ") = "
200  << calc_err << std::endl;
201  if( calc_err >= warning_tol() ) {
202  if(out && print_tests() >= PRINT_ALL)
203  *out
204  << std::endl
205  << ( calc_err >= error_tol() ? "Error" : "Warning" )
206  << ", rel_err(sum(alpha*v),sum(t2)) = "
207  << "rel_err(" << sum_av << "," << sum_t2 << ") = "
208  << calc_err
209  << " exceeded "
210  << ( calc_err >= error_tol() ? "error_tol" : "warning_tol" )
211  << " = "
212  << ( calc_err >= error_tol() ? error_tol() : warning_tol() )
213  << std::endl;
214  if(calc_err >= error_tol()) {
215  if(dump_all() && print_tests() >= PRINT_ALL) {
216  *out << "\nalpha = " << alpha << std::endl;
217  *out << "\nv =\n" << *v;
218  *out << "\nt1 =\n" << *t2;
219  *out << "\nt2 =\n" << *t2;
220  }
221  lresult = false;
222  }
223  }
224  if(!lresult) result = false;
225  }
226  if(!result) success = false;
227  if( out && print_tests() == PRINT_MORE )
228  *out << " : " << ( result ? "passed" : "failed" )
229  << std::endl;
230  }
231  }
232 
233  if( out && print_tests() == PRINT_BASIC )
234  *out << " : " << ( success ? "passed" : "failed" );
235 
236  return success;
237 }
238 
239 } // end namespace AbstractLinAlgPack
virtual const VectorSpace & space_rows() const =0
Vector space for vectors that are compatible with the rows of the matrix.
virtual value_type norm_1() const
One norm. ||v||_1 = sum( |v(i)|, i = 1,,,this->dim() )
bool test_matrix(const MatrixOpNonsing &M, const char M_name[], std::ostream *out)
Test a MatrixOpNonsing object.
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
T * get() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
Side
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
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
MatrixOpNonsingTester(ETestLevel test_level=TEST_LEVEL_2_BLAS, EPrintTestLevel print_tests=PRINT_NONE, 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)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
virtual index_type dim() const
Return the dimension of this vector.
virtual const VectorSpace & space_cols() const =0
Vector space for vectors that are compatible with the columns of the matrix.
Abstract interface for mutable coordinate vectors {abstract}.
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...
Print everything all the tests in great detail but output is independent of problem size...
virtual vec_mut_ptr_t create_member() const =0
Create a vector member from the vector space.
Transp
void random_vector(value_type l, value_type u, VectorMutable *v)
Generate a random vector with elements uniformly distrubuted elements.
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)
v_lhs += V_rhs.