MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DenseLinAlgPack_TestGenMatrixClass.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 <iomanip>
43 #include <ostream>
44 #include <vector>
45 #include <typeinfo>
46 
47 #include "DenseLinAlgPack_TestDenseLinAlgPack.hpp"
48 #include "DenseLinAlgPack_DMatrixClass.hpp"
49 #include "DenseLinAlgPack_DVectorOut.hpp"
50 #include "DenseLinAlgPack_DMatrixOut.hpp"
51 #include "DenseLinAlgPack_MatVecCompare.hpp"
52 
53 namespace {
54 
55 using DenseLinAlgPack::sqrt_eps;
56 
57 // Check consistency of row(), col(), diag() and operator()().
58 template<class M_t>
59 void check_access( M_t& M, typename M_t::size_type row_offset, typename M_t::size_type col_offset
60  , std::ostream* out, bool* success )
61 {
62  if(out)
63  *out << "Checking M(i,j) == M.row(i)(j) == M.col(j)(i) == "
64  << "M.diag(...)(...) == "
65  << "(i + "<<row_offset<<") + 0.1*(j+"<<col_offset<<") : ";
66 
67  bool result = true;
68 
69  for( typename M_t::size_type i = 1; i <= M.rows(); ++i ) {
70  for( typename M_t::size_type j = 1; j <= M.rows(); ++j ) {
71  const typename M_t::value_type
72  Mij = M(i,j);
73  typename M_t::value_type
74  val = (i+row_offset)+0.1*(j+col_offset);
75  if( ::fabs(Mij-val) > sqrt_eps ) {
76  result = false;
77  if(out) *out << "(M("<<i<<","<<j<<") -> "<<Mij<<") != "<<val<<std::endl;
78  }
79  if( Mij != (val = M.row(i)(j)) ) {
80  result = false;
81  if(out) *out << "M("<<i<<","<<j<<") != (M.row("<<i<<")("<<j<<") -> "<<val<<")\n";
82  }
83  if( Mij != (val = M.col(j)(i)) ) {
84  result = false;
85  if(out) *out << "M("<<i<<","<<j<<") != (M.col("<<j<<")("<<i<<") -> "<<val<<")\n";
86  }
87  const int k = ( i > j ? -i + j : j - i );
88  const typename M_t::size_type k_i = ( i > j ? j : i );
89  if( Mij != (val = M.diag(k)(k_i) ) ) {
90  result = false;
91  if(out) *out << "M("<<i<<","<<j<<") != (M.diag("<<k<<")("<<k_i<<") -> "<<val<<")\n";
92  }
93  }
94  }
95  if(out) *out << result << std::endl;
96  if(!result) *success = false;
97 }
98 
99 // Print out a string for overlap
100 const char* overlap_str( DenseLinAlgPack::EOverLap overlap ) {
101  switch(overlap) {
102  case DenseLinAlgPack::NO_OVERLAP:
103  return "NO_OVERLAP";
104  case DenseLinAlgPack::SOME_OVERLAP:
105  return "SOME_OVERLAP";
106  case DenseLinAlgPack::SAME_MEM:
107  return "SAME_MEM";
108  }
109  return "Invalid value for EOverLap";
110 }
111 
112 } // end namespace
113 
114 bool DenseLinAlgPack::TestingPack::TestGenMatrixClass(std::ostream* out)
115 {
116 
117  using DenseLinAlgPack::comp;
118  using DenseLinAlgPack::sqrt_eps;
119 
120  bool success = true;
121  bool result;
122 
123  if(out)
124  *out << "\n****************************************************"
125  << "\n*** Testing DMatrix and DMatrixSlice classes ***"
126  << "\n****************************************************\n"
127  << std::boolalpha;
128 
129  try {
130 
131  const size_type
132  m = 6,
133  n = 8;
134 
135  const value_type
136  ptr[m*n] =
137  { 1.1, 2.1, 3.1, 4.1, 5.1, 6.1,
138  1.2, 2.2, 3.2, 4.2, 5.2, 6.2,
139  1.3, 2.3, 3.3, 4.3, 5.3, 6.3,
140  1.4, 2.4, 3.4, 4.4, 5.4, 6.4,
141  1.5, 2.5, 3.5, 4.5, 5.5, 6.5,
142  1.6, 2.6, 3.6, 4.6, 5.6, 6.6,
143  1.7, 2.7, 3.7, 4.7, 5.7, 6.7,
144  1.8, 2.8, 3.8, 4.8, 5.8, 6.8 };
145 
146  // /////////////////////////////
147  // Test Constructors
148 
149  if(out)
150  *out << "\n***\n*** Testing constructors\n***\n";
151 
152  // DMatrixSlice
153  if(out) *out << "\nGenMatrixSlice gms1;\n";
154  DMatrixSlice gms1;
155  if(out) *out << "gms1 =\n" << gms1;
156  update_success( result = (gms1.rows() == 0 && gms1.cols() == 0 ), &success );
157  if(out)
158  *out << "((gms1.rows() -> "<<gms1.rows()
159  << ") == 0 && (gms1.cols() -> "<<gms1.cols()<<") == 0 ) : "
160  << result << std::endl;
161 
162 
163  if(out) *out << "\nGenMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );\n";
164  const DMatrixSlice gms2( const_cast<value_type*>(ptr), m*n, m, m, n );
165  if(out) *out << "gms2 =\n" << gms2;
166 
167  if(out) *out << "\nGenMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );\n";
168  const DMatrixSlice gms3( const_cast<DMatrixSlice&>(gms2), Range1D(1,m), Range1D(1,n) );
169  if(out) *out << "gms3 =\n" << gms3;
170 
171  // DMatrix
172 
173  if(out) *out << "\nGenMatrix gm1;\n";
174  DMatrix gm1;
175  if(out) *out << "gm1 =\n" << gm1();
176  update_success( result = (gm1.rows() == 0 && gm1.cols() == 0 ), &success );
177  if(out)
178  *out << "((gm1.rows() -> "<<gm1.rows()
179  << ") == 0 && (gm1.cols() -> "<<gm1.cols()<<") == 0 ) : "
180  << result << std::endl;
181 
182  if(out) *out << "\nGenMatrix gm2(m,n);\n";
183  DMatrix gm2(m,n);
184  if(out) *out << "gm2 =\n" << gm2();
185 
186  if(out) *out << "\nGenMatrix gm3(1.0,m,n);\n";
187  DMatrix gm3(1.0,m,n);
188  if(out) *out << "gm3 =\n" << gm3();
189  update_success( result = comp( gm3(), 1.0 ), &success );
190  if(out) *out << "gm3 == 1.0 : " << result << std::endl;
191 
192  if(out) *out << "\nGenMatrix gm4(ptr,m,n);\n";
193  DMatrix gm4(ptr,m,n);
194  if(out) *out << "gm4 =\n" << gm4();
195 
196  if(out) *out << "\nGenMatrix gm5(gms2);\n";
197  DMatrix gm5(gms2);
198  if(out) *out << "gm5 =\n" << gm5();
199 
200  // ////////////////////////////
201  // Test DMatrixSlice binding
202 
203  if(out)
204  *out << "\n***\n*** Testing DMatrixSlice binding\n***\n";
205 
206  if(out) *out << "\ngms1.bind(gm4());\n";
207  gms1.bind(gm4());
208  if(out) *out << "gms1 =\n" << gms1();
209 
210  // ////////////////////////////
211  // Test DMatrix resizing
212 
213  if(out)
214  *out << "\n***\n*** Testing DMatrix resizing\n***\n";
215 
216  if(out) *out << "\ngm1.resize(m,n,1.0);\n";
217  gm1.resize(m,n,1.0);
218  if(out) *out << "gm1 =\n" << gm1();
219  update_success( result = comp( gm1(), 1.0 ), &success );
220  if(out) *out << "gm1 == 1.0 : " << result << std::endl;
221 
222  // ///////////////////////////////////////////////
223  // Test row, col, diag access and element access
224 
225  // DMatrixSlice
226 
227  if(out)
228  *out << "\n***\n*** Testing row, col, diag access and element access\n***\n";
229 
230  if(out) *out << "\nLet M = gms1\n";
231  check_access( gms1, 0, 0, out, &success );
232 
233  if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gms1)\n";
234  check_access( const_cast<const DMatrixSlice&>(gms1), 0, 0, out, &success );
235 
236  // DMatrix
237 
238  if(out) *out << "\nLet M = gm4\n";
239  check_access( gm4, 0, 0, out, &success );
240 
241  if(out) *out << "\nLet M = const_cast<const DMatrix&>(gm4)\n";
242  check_access( const_cast<const DMatrix&>(gm4), 0, 0, out, &success );
243 
244  // ////////////////////////////
245  // Test submatrix access
246 
247  if(out)
248  *out << "\n***\n*** Testing submatrix access\n***\n";
249 
250  if(out) *out << "\nRange1D r_rng(2,m-1), c_rng(2,n-1);\n";
251  Range1D r_rng(2,m-1), c_rng(2,n-1);
252 
253  // DMatrixSlice
254 
255  if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng)\n";
256  gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
257  if(out) *out << "M =\n" << gms1;
258  check_access( gms1, 1, 1, out, &success );
259 
260  if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(r_rng,c_rng))\n";
261  gms1.bind( const_cast<DMatrixSlice&>(gms2)(r_rng,c_rng) );
262  if(out) *out << "M =\n" << gms1;
263  check_access( gms1, 1, 1, out, &success );
264 
265  if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1)\n";
266  gms1.bind(const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
267  if(out) *out << "M =\n" << gms1;
268  check_access( gms1, 1, 1, out, &success );
269 
270  if(out) *out << "\nLet M = const_cast<DMatrixSlice&>(gms2(2,m-1,2,n-1))\n";
271  gms1.bind( const_cast<DMatrixSlice&>(gms2)(2,m-1,2,n-1) );
272  if(out) *out << "M =\n" << gms1;
273  check_access( gms1, 1, 1, out, &success );
274 
275  // DMatrix
276 
277  if(out) *out << "\nLet M = gm4(r_rng,c_rng)\n";
278  gms1.bind( gm4(r_rng,c_rng) );
279  if(out) *out << "M =\n" << gms1;
280  check_access( gms1, 1, 1, out, &success );
281 
282  if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(r_rng,c_rng)\n";
283  gms1.bind( const_cast<const DMatrix&>(gm4)(r_rng,c_rng) );
284  if(out) *out << "M =\n" << gms1;
285  check_access( gms1, 1, 1, out, &success );
286 
287  if(out) *out << "\nLet M = gm4(2,m-1,2,n-1)\n";
288  gms1.bind( gm4(2,m-1,2,n-1) );
289  if(out) *out << "M =\n" << gms1;
290  check_access( gms1, 1, 1, out, &success );
291 
292  if(out) *out << "\nLet M = const_cast<const DMatrixSlice&>(gm4)(2,m-1,2,n-1)\n";
293  gms1.bind( const_cast<const DMatrix&>(gm4)(2,m-1,2,n-1) );
294  if(out) *out << "M =\n" << gms1;
295  check_access( gms1, 1, 1, out, &success );
296 
297  // ////////////////////
298  // Test matrix overlap
299 
300  if(out)
301  *out << "\n***\n*** matrix overlap\n***\n";
302 
303  EOverLap ovlap;
304 
305  // DMatrixSlice
306 
307  if(out) *out << "(gms2.overlap(gms2) -> ";
308  ovlap = gms2.overlap(gms2);
309  result = update_success( ovlap == SAME_MEM, &success );
310  if(out) *out << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
311 
312  if(out) *out << "(gms2.overlap(gms2(r_rng,c_rng)) -> ";
313  ovlap = gms2.overlap(gms2(r_rng,c_rng));
314  result = update_success( ovlap == SOME_OVERLAP, &success );
315  if(out) *out << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
316 
317  if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n)) -> ";
318  ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2,m,n/2,n));
319  result = update_success( ovlap == SOME_OVERLAP, &success );
320  if(out) *out << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
321 
322  if(out) *out << "(gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n)) -> ";
323  ovlap = gms2(1,m/2,1,n/2).overlap(gms2(m/2+1,m,n/2+1,n));
324  result = update_success( ovlap == NO_OVERLAP, &success );
325  if(out) *out << overlap_str(ovlap) << ") == NO_OVERLAP : " << result << std::endl;
326 
327  // DMatrix
328 
329  if(out) *out << "(gm4.overlap(gm4) -> ";
330  ovlap = gm4.overlap(gm4());
331  result = update_success( ovlap == SAME_MEM, &success );
332  if(out) *out << overlap_str(ovlap) << ") == SAME_MEM : " << result << std::endl;
333 
334  if(out) *out << "(gm4.overlap(gm4(r_rng,c_rng)) -> ";
335  ovlap = gm4.overlap(gm4(r_rng,c_rng));
336  result = update_success( ovlap == SOME_OVERLAP, &success );
337  if(out) *out << overlap_str(ovlap) << ") == SOME_OVERLAP : " << result << std::endl;
338 
339  // //////////////////////////////////////////////////////
340  // Test vector overlap (continuation of vector testing)
341  //
342  // ToDo: Finish this someday once you get it figured out.
343 
344  // ///////////////////////////
345  // Test assignment operators
346 
347  if(out)
348  *out << "\n***\n*** assignment operators\n***\n";
349 
350  // DMatrixSlice
351 
352  if(out) *out << "\ngms1.bind(gm1());\n";
353  gms1.bind(gm1());
354 
355  if(out) *out << "\ngms1 = 2.0;\n";
356  gms1 = 2.0;
357  if(out) *out << "gms1 =\n" << gms1;
358  update_success( result = comp(gms1,2.0), &success );
359  if(out) *out << "gms1 == 2.0 : " << result << std::endl;
360 
361  if(out) *out << "\ngms1 = gms2;\n";
362  gms1 = gms2;
363  if(out) *out << "gms1 =\n" << gms1;
364  update_success( result = comp(gms1,gms2), &success );
365  if(out) *out << "gms1 == gms2 : " << result << std::endl;
366 
367  // DMatrix
368 
369  if(out) *out << "\ngm1 = 3.0;\n";
370  gm1 = 3.0;
371  if(out) *out << "gm1 =\n" << gm1();
372  update_success( result = comp(gm1(),3.0), &success );
373  if(out) *out << "gm1 == 3.0 : " << result << std::endl;
374 
375  if(out) *out << "\ngm1.resize(0,0); gm1 = gms2;\n";
376  gm1.resize(0,0);
377  gm1 = gms2;
378  if(out) *out << "gm1 =\n" << gm1();
379  update_success( result = comp(gm1(),gms2()), &success );
380  if(out) *out << "gm1 == gms2 : " << result << std::endl;
381 
382  if(out) *out << "\ngm1.resize(0,0); gm1 = gm4;\n";
383  gm1.resize(0,0);
384  gm1 = gm4;
385  if(out) *out << "gm1 =\n" << gm1();
386  update_success( result = comp(gm1(),gm4()), &success );
387  if(out) *out << "gm1 == gm4 : " << result << std::endl;
388 
389  } // end try
390  catch( const std::exception& excpt ) {
391  success = false;
392  if(out)
393  (*out) << "\nError, a standard exception was thrown: "
394  << typeName(excpt) << ": "
395  << excpt.what() << std::endl;
396  }
397  catch(...) {
398  success = false;
399  if(out)
400  (*out) << "\nError, an unknown exception was thrown\n";
401  }
402 
403  if(out) {
404  if(success)
405  (*out)
406  << "\n*** Congradulations, DMatrix and DMatrixSlice seem to check out. ***\n";
407  else
408  (*out)
409  << "\n*** Oops, all of the tests for DMatrix and DMatrixSlice "
410  "where not successful. ***\n";
411  }
412 
413  return success;
414 }
415 
DenseLinAlgPack::size_type
Teuchos::Ordinal size_type
Typedef for the size type of elements that are used by the library.
Definition: DenseLinAlgPack_Options.hpp:104
LinAlgOpPack::size_type
AbstractLinAlgPack::size_type size_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:54
typeName
std::string typeName(const T &t)
DenseLinAlgPack_MatVecCompare.hpp
DenseLinAlgPack_DMatrixOut.hpp
DenseLinAlgPack_TestDenseLinAlgPack.hpp
RangePack::Range1D
. One-based subregion index range class.
Definition: RangePack_Range1D.hpp:76
QPKWIK_Output::out
std::ostream * out
M
const LAPACK_C_Decl::f_int & M
Definition: DenseLinAlgPack_LAPACK_C_Decl.cpp:82
DenseLinAlgPack_DMatrixClass.hpp
DenseLinAlgPack::EOverLap
EOverLap
Enumeration for returning the amount of overlap between two objects.
Definition: DenseLinAlgPack_Types.hpp:81
LinAlgOpPack::value_type
AbstractLinAlgPack::value_type value_type
Definition: AbstractLinAlgPack_LinAlgOpPack.hpp:56
DenseLinAlgPack::value_type
FortranTypes::f_dbl_prec value_type
Typedef for the value type of elements that is used for the library.
Definition: DenseLinAlgPack_Options.hpp:100
DenseLinAlgPack_DVectorOut.hpp
TestingHelperPack::update_success
bool update_success(bool result_check, bool *success)
Helper function for updating a flag for if an operation returned false.
Definition: TestingHelperPack_update_success.cpp:48
DenseLinAlgPack::comp
bool comp(const DVectorSlice &vs1, const DVectorSlice &vs2)
Definition: DenseLinAlgPack_MatVecCompare.cpp:72
n
int n
Generated on Wed Dec 2 2015 08:50:25 for MOOCHO (Single Doxygen Collection) by   doxygen 1.8.6