Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
test/DenseMatrix/cxx_main_sym.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
46 #include "Teuchos_Version.hpp"
47 
48 #define OTYPE int
49 #define STYPE std::complex<double>
50 
51 template<typename TYPE>
52 int PrintTestResults(std::string, TYPE, TYPE, bool);
53 
54 int ReturnCodeCheck(std::string, int, int, bool);
55 
59 
60 int main(int argc, char* argv[])
61 {
62 
63  int i;
64  bool verbose = 0;
65  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
66 
67  if (verbose)
68  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
69 
70  int numberFailedTests = 0;
71  int returnCode = 0;
72  std::string testName = "";
73 
74  if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL SYMMETRIC DENSE MATRIX **********"<<std::endl<<std::endl;
75 
76  // default constructor test
77  SDMatrix DefConTest;
78  if (verbose) std::cout <<"default constructor -- construct empty matrix ";
79  if ( DefConTest.values()!=NULL || DefConTest.numCols()!=0 || DefConTest.numRows()!=0 ||DefConTest.stride()!=0 ||DefConTest.empty()!=true ) {
80  if (verbose) std::cout << "unsuccessful."<<std::endl;
81  numberFailedTests++;
82  } else {
83  if (verbose) std::cout << "successful."<<std::endl;
84  }
85 
86  // constructor 1 (matrix w/ dimension but empty)
87 
88  SDMatrix Con1Test( 4 );
89  if (verbose) std::cout <<"constructor 1 -- empty matrix with given dimensions ";
90  if ( Con1Test.numRows()!=4 || Con1Test.numCols()!=4 || Con1Test( 1, 2 )!=0.0 ) {
91  if (verbose) std::cout << "unsuccessful."<<std::endl;
92  numberFailedTests++;
93  } else {
94  if (verbose) std::cout << "successful."<<std::endl;
95  }
96 
97  // constructor 2 (from array) tests
98 
99  STYPE a[9];
100  for(i = 0; i < 9; i++)
101  {
102  a[i] = i;
103  }
104  SDMatrix Con2Test1ExpRes;
105  Con2Test1ExpRes.shape(3);
106  Con2Test1ExpRes(0, 0) = 0;
107  Con2Test1ExpRes(1, 0) = 1; Con2Test1ExpRes(1, 1) = 4;
108  Con2Test1ExpRes(2, 0) = 2; Con2Test1ExpRes(2, 1) = 5; Con2Test1ExpRes(2, 2) = 8;
109 
110  // Create another lower triangular matrix with a view of 'a'.
111  SDMatrix Con2Test1(Teuchos::View, false, a, 3, 3);
112  numberFailedTests += PrintTestResults("constructor 2 -- construct matrix from array subrange", Con2Test1, Con2Test1ExpRes, verbose);
113 
114 
115  // constructor 3 (copy constructor)
116 
117  SDMatrix Con3TestCopy( Con2Test1ExpRes );
118  if(verbose) std::cout <<"constructor 3 -- copy constructor ";
119  if ( Con3TestCopy != Con2Test1ExpRes ) {
120  if (verbose) std::cout << "unsuccessful."<<std::endl;
121  numberFailedTests++;
122  } else {
123  if (verbose) std::cout << "successful."<<std::endl;
124  }
125 
126  SDMatrix Con3TestCopyTrans( Con2Test1ExpRes );
127  Con3TestCopyTrans.setUpper();
128  if(verbose) std::cout <<"constructor 3 -- copy constructor (upper active storage) ";
129  if ( Con3TestCopyTrans(2, 0) != Con2Test1ExpRes(2, 0) ) {
130  if (verbose) std::cout << "unsuccessful."<<std::endl;
131  numberFailedTests++;
132  } else {
133  if (verbose) std::cout << "successful."<<std::endl;
134  }
135 
136  // constructor 4 (submatrix)
137 
138  SDMatrix Con4TestOrig(Teuchos::Copy, false, a, 3, 3);
139  SDMatrix Con4TestSubmatrix;
140  Con4TestSubmatrix.shape( 2 );
141  Con4TestSubmatrix(0, 0) = 4;
142  Con4TestSubmatrix(1, 0) = 5; Con4TestSubmatrix(1, 1) = 8;
143  SDMatrix Con4TestCopy1(Teuchos::Copy, Con4TestOrig, 2, 1);
144  numberFailedTests += PrintTestResults("constructor 4 -- submatrix copy", Con4TestCopy1, Con4TestSubmatrix, verbose);
145  SDMatrix Con4TestCopy2(Teuchos::Copy, Con4TestOrig, 3, 0);
146  numberFailedTests += PrintTestResults("constructor 4 -- full matrix copy", Con4TestCopy2, Con4TestOrig, verbose);
147  SDMatrix Con4TestView1(Teuchos::View, Con4TestOrig, 2, 1);
148  numberFailedTests += PrintTestResults("constructor 4 -- full matrix view", Con4TestView1, Con4TestSubmatrix, verbose);
149  SDMatrix Con4TestView2(Teuchos::View, Con4TestOrig, 3, 0);
150  numberFailedTests += PrintTestResults("constructor 4 -- submatrix view", Con4TestView2, Con4TestOrig, verbose);
151 
152  // Norm Tests
153 
154  SDMatrix AAA;
155  AAA.shape( 3 );
156  AAA(0, 0) = 8;
157  AAA(1, 0) = 1; AAA(1, 1) = 8;
158  AAA(2, 0) = 2; AAA(2, 1) = 3; AAA(2, 2) = 8;
159  SDMatrix BBB;
160  numberFailedTests += PrintTestResults("normOne of a 3x3", AAA.normOne(), 13.0, verbose);
161  numberFailedTests += PrintTestResults("normInf of a 3x3", AAA.normInf(), 13.0, verbose);
163  numberFailedTests += PrintTestResults("normFrobenius of a 3x3", AAA.normFrobenius(), 3.0, verbose);
164  numberFailedTests += PrintTestResults("normOne of a 0x0", BBB.normOne(), 0.0, verbose);
165  numberFailedTests += PrintTestResults("normInf of a 0x0", BBB.normInf(), 0.0, verbose);
166  numberFailedTests += PrintTestResults("normFrobenius of a 0x0", BBB.normFrobenius(), 0.0, verbose);
167 
168  // Multiplication Tests
169 
170  // Reset values of AAA.
171  AAA(0, 0) = 8;
172  AAA(1, 0) = 1; AAA(1, 1) = 8;
173  AAA(2, 0) = 2; AAA(2, 1) = 3; AAA(2, 2) = 8;
174 
175  DMatrix My_Prod( 4, 3 ), My_GenMatrix( 4, 3 );
176  My_GenMatrix = Teuchos::ScalarTraits<STYPE>::one();
177 
178  // Matrix multiplication ( My_Prod = 1.0*My_GenMatrix*My_Matrix )
179  My_Prod.multiply( Teuchos::RIGHT_SIDE, 1.0, AAA, My_GenMatrix, 0.0 );
180  numberFailedTests += PrintTestResults("multiply() -- general times symmetric matrix (storage = lower tri)", My_Prod.normOne(), 52.0, verbose);
181  AAA.setUpper();
182  AAA(2, 1) = 1.0;
183  My_Prod.multiply( Teuchos::RIGHT_SIDE, 1.0, AAA, My_GenMatrix, 0.0 );
184  numberFailedTests += PrintTestResults("multiply() -- general times symmetric matrix (storage = upper tri)", My_Prod.normOne(), 44.0, verbose);
185 
186  // Set Method Tests.
187 
188  SDMatrix CCC( 5 );
189  // Randomize the entries in CCC.
190  testName = "random() -- enter random entries into matrix";
191  returnCode = CCC.random();
192  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
193  // Set the entries of CCC to 1.0.
194  testName = "putScalar() -- set every entry of this matrix to 1.0";
195  returnCode = CCC.putScalar(Teuchos::ScalarTraits<STYPE>::one());
196  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
197  // Check assignment operator.
198  SDMatrix CCC2( 5 );
199  CCC2.assign( CCC );
200  if (verbose) std::cout << "assign() -- copy the values of an input matrix ";
201  if ( CCC( 3, 4 ) == Teuchos::ScalarTraits<STYPE>::one() ) {
202  if (verbose) std::cout<< "successful" <<std::endl;
203  } else {
204  if (verbose) std::cout<< "unsuccessful" <<std::endl;
205  numberFailedTests++;
206  }
207  // Create a view into a submatrix of CCC
208  SDMatrix CCCview( Teuchos::View, CCC, 3 );
209  SDMatrix CCCtest1( 2 );
210  CCCtest1 = CCCview;
211  if (verbose) std::cout << "operator= -- small(empty) = large(view) ";
212  if (CCCtest1.numRows()==3 && CCCtest1.values()==CCC.values()) {
213  if (verbose) std::cout<< "successful" <<std::endl;
214  } else {
215  if (verbose) std::cout<< "unsuccessful" <<std::endl;
216  numberFailedTests++;
217  }
218  CCCtest1 = CCC;
219  if (verbose) std::cout << "operator= -- small(view) = large(copy) ";
220  if (CCCtest1.numRows()==5 && CCCtest1.values()!=CCC.values()) {
221  if (verbose) std::cout<< "successful"<<std::endl;
222  } else {
223  if (verbose) std::cout<< "unsuccessful"<<std::endl;
224  numberFailedTests++;
225  }
226  SDMatrix CCCtest2( 2 );
228  CCCtest1 = CCCtest2;
229  if (verbose) std::cout << "operator= -- large(copy) = small(copy) ";
230  if (CCCtest1.numRows()==2 ) {
231  if (verbose) std::cout<< "successful"<<std::endl;
232  } else {
233  if (verbose) std::cout<< "unsuccessful"<<std::endl;
234  numberFailedTests++;
235  }
236  CCCtest1 = CCCview;
237  if (verbose) std::cout << "operator= -- large(copy) = small(view) ";
238  if (CCCtest1.numRows()==3 && CCCtest1.stride()==5) {
239  if(verbose) std::cout<<"successful" <<std::endl;
240  } else {
241  if (verbose) std::cout<<"unsuccessful"<<std::endl;
242  numberFailedTests++;
243  }
244 
245  SDMatrix CCCtest3( CCCview );
246  CCCtest1 += CCCtest3;
247  if (verbose) std::cout << "operator+= -- add two matrices of the same size, but different leading dimension ";
248  if (CCCtest1(1,1)==2.0) {
249  if(verbose) std::cout<<"successful" <<std::endl;
250  } else {
251  if (verbose) std::cout<<"unsuccessful"<<std::endl;
252  numberFailedTests++;
253  }
254  if (verbose) std::cout << "operator+= -- add two matrices of different size (nothing should change) ";
255  CCCtest1 += CCC;
256  if (CCCtest1(1,1)==2.0) {
257  if(verbose) std::cout<<"successful" <<std::endl;
258  } else {
259  if (verbose) std::cout<<"unsuccessful"<<std::endl;
260  numberFailedTests++;
261  }
262 
263  // Scale Tests.
264 
265  SDMatrix ScalTest( 8 );
267  // Scale the entries by 8, it should be 8.
268  // The matrix is lower triangular, by default, so check a lower triangular entry.
269  if (verbose) std::cout << "operator*= -- scale matrix by some number ";
270  ScalTest *= 8.0;
271  if (ScalTest(7, 1) == 8.0) {
272  if (verbose) std::cout<< "successful." <<std::endl;
273  } else {
274  if (verbose) std::cout<< "unsuccessful." <<std::endl;
275  numberFailedTests++;
276  }
277 
278 
279  // Matrix Triple-Product Test
281  DMatrix W(3,2);
282  SDMatrix A1(2), A2(3);
283  A1(0,0) = 1.0, A1(1,1) = 2.0;
284  A2(0,0) = 1.0, A2(1,1) = 2.0, A2(2,2) = 3.00;
286 
287  SDMatrix C1upper(3), C1lower(3), C2upper(2), C2lower(2);
288  C1upper.setUpper(); C2upper.setUpper();
289  C1lower.setLower(); C2lower.setLower();
290 
291  // Test all combinations of triple products.
292 
293  // These should return a matrix with 1.5 in all entries
294  STYPE C1result = 1.5*Teuchos::ScalarTraits<STYPE>::one();
295  Teuchos::symMatTripleProduct<OTYPE,STYPE>( Teuchos::NO_TRANS, alpha, A1, W, C1upper );
296  Teuchos::symMatTripleProduct<OTYPE,STYPE>( Teuchos::NO_TRANS, alpha, A1, W, C1lower );
297 
298  // These should return a matrix with 3 in all entries
299  STYPE C2result = 3.0*Teuchos::ScalarTraits<STYPE>::one();
300  Teuchos::symMatTripleProduct<OTYPE,STYPE>( Teuchos::TRANS, alpha, A2, W, C2upper );
301  Teuchos::symMatTripleProduct<OTYPE,STYPE>( Teuchos::TRANS, alpha, A2, W, C2lower );
302 
303  if (verbose) std::cout << "triple product -- compute C = W'*A*W or C = W*A*W' ";
304  if (C1upper(2,1)==C1result && C1lower(1,2)==C1result && C2upper(1,0)==C2result && C2lower(0,1)==C2result) {
305  if (verbose) std::cout<< "successful." <<std::endl;
306  } else {
307  if (verbose) std::cout<< "unsuccessful." <<std::endl;
308  numberFailedTests++;
309  }
310 
311 
312 
313  //
314  // If a test failed output the number of failed tests.
315  //
316  if(numberFailedTests > 0)
317  {
318  if (verbose) {
319  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
320  std::cout << "End Result: TEST FAILED" << std::endl;
321  return -1;
322  }
323  }
324  if(numberFailedTests == 0)
325  std::cout << "End Result: TEST PASSED" << std::endl;
326 
327  return 0;
328 }
329 
330 template<typename TYPE>
331 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
332 {
333  int result;
334  if(calculatedResult == expectedResult)
335  {
336  if(verbose) std::cout << testName << " successful." << std::endl;
337  result = 0;
338  }
339  else
340  {
341  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
342  result = 1;
343  }
344  return result;
345 }
346 
347 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
348 {
349  int result;
350  if(expectedResult == 0)
351  {
352  if(returnCode == 0)
353  {
354  if(verbose) std::cout << testName << " test successful." << std::endl;
355  result = 0;
356  }
357  else
358  {
359  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
360  result = 1;
361  }
362  }
363  else
364  {
365  if(returnCode != 0)
366  {
367  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
368  result = 0;
369  }
370  else
371  {
372  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
373  result = 1;
374  }
375  }
376  return result;
377 }
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
Templated serial dense matrix class.
bool empty() const
Returns whether this matrix is empty.
Teuchos::SerialDenseMatrix< OTYPE, STYPE > DMatrix
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
Teuchos::SerialDenseVector< OTYPE, STYPE > DVector
This class creates and provides basic support for dense vectors of templated type as a specialization...
This structure defines some basic traits for a scalar field type.
void setLower()
Specify that the lower triangle of the this matrix should be used.
Teuchos::SerialSymDenseMatrix< OTYPE, STYPE > SDMatrix
void setUpper()
Specify that the upper triangle of the this matrix should be used.
OrdinalType numRows() const
Returns the row dimension of this matrix.
int main(int argc, char *argv[])
int PrintTestResults(std::string, TYPE, TYPE, bool)
std::string Teuchos_Version()
int ReturnCodeCheck(std::string, int, int, bool)
Templated serial, dense, symmetric matrix class.
Templated serial dense vector class.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
static T one()
Returns representation of one for this scalar type.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
This class creates and provides basic support for dense rectangular matrix of templated type...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.