Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
cxx_main_solver.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_ScalarTraits.hpp"
47 #include "Teuchos_RCP.hpp"
48 #include "Teuchos_Version.hpp"
49 
53 
54 #define OTYPE int
55 #ifdef HAVE_TEUCHOS_COMPLEX
56 #define STYPE std::complex<double>
57 #else
58 #define STYPE double
59 #endif
60 
61 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
62 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
63 #ifdef HAVE_TEUCHOS_COMPLEX
64 #define SCALARMAX STYPE(10,0)
65 #else
66 #define SCALARMAX STYPE(10)
67 #endif
68 
69 template<typename TYPE>
70 int PrintTestResults(std::string, TYPE, TYPE, bool);
71 
72 int ReturnCodeCheck(std::string, int, int, bool);
73 
76 
77 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
78 template<typename TYPE>
79 TYPE GetRandom(TYPE, TYPE);
80 
81 // Returns a random integer between the two input parameters, inclusive
82 template<>
83 int GetRandom(int, int);
84 
85 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
86 template<>
87 double GetRandom(double, double);
88 
89 template<typename T>
90 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
91 
92 // Generates random matrices and vectors using GetRandom()
95 
96 // Compares the difference between two vectors using relative euclidean norms
97 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
99  const SerialDenseVector<OTYPE,STYPE>& Vector2,
101 
102 int main(int argc, char* argv[])
103 {
104  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
105 
106  int n=10, m=8;
107  (void) m; // forestall "unused variable" compiler warning
108  MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
109 
110  bool verbose = 0;
111  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
112 
113  if (verbose)
114  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
115 
116  int numberFailedTests = 0;
117  int returnCode = 0;
118  std::string testName = "", testType = "";
119 
120 #ifdef HAVE_TEUCHOS_COMPLEX
121  testType = "COMPLEX";
122 #else
123  testType = "REAL";
124 #endif
125 
126  if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
127 
128  // Create dense matrix and vector.
131  DVector xhat(n), b(n), bt(n);
132 
133  // Compute the right-hand side vector using multiplication.
135  testName = "Generating right-hand side vector using A*x, where x is a random vector:";
136  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
137 
139  testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
140  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
141 
142 #ifdef HAVE_TEUCHOS_COMPLEX
143  DVector bct(n);
145  testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
146  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
147 #endif
148 
149  // Fill the solution vector with zeros.
150  xhat.putScalar( ScalarTraits<STYPE>::zero() );
151 
152  // Create a serial dense solver.
154 
155  // Pass in matrix and vectors
156  solver1.setMatrix( A1 );
157  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
158 
159  // Test1: Simple factor and solve
160  returnCode = solver1.factor();
161  testName = "Simple solve: factor() random A:";
162  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
163 
164  // Non-transpose solve
165  returnCode = solver1.solve();
166  testName = "Simple solve: solve() random A (NO_TRANS):";
167  numberFailedTests += CompareVectors( *x1, xhat, tol );
168  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
169 
170  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
171  xhat.putScalar( ScalarTraits<STYPE>::zero() );
172  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
174  returnCode = solver1.solve();
175  testName = "Simple solve: solve() random A (TRANS):";
176  numberFailedTests += CompareVectors( *x1, xhat, tol );
177  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
178 
179 #ifdef HAVE_TEUCHOS_COMPLEX
180  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
181  xhat.putScalar( ScalarTraits<STYPE>::zero() );
182  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
184  returnCode = solver1.solve();
185  testName = "Simple solve: solve() random A (CONJ_TRANS):";
186  numberFailedTests += CompareVectors( *x1, xhat, tol );
187  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
188 #endif
189 
190  // Test2: Invert the matrix, inverse should be in A.
191  returnCode = solver1.invert();
192  testName = "Simple solve: invert() random A:";
193  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
194 
195  // Compute the solution vector using multiplication and the inverse.
197  testName = "Computing solution using inverted random A (NO_TRANS):";
198  numberFailedTests += CompareVectors( *x1, xhat, tol );
199  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
200 
201  returnCode = xhat.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bt, ScalarTraits<STYPE>::zero());
202  testName = "Computing solution using inverted random A (TRANS):";
203  numberFailedTests += CompareVectors( *x1, xhat, tol );
204  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
205 
206 #ifdef HAVE_TEUCHOS_COMPLEX
207  returnCode = xhat.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bct, ScalarTraits<STYPE>::zero());
208  testName = "Computing solution using inverted random A (CONJ_TRANS):";
209  numberFailedTests += CompareVectors( *x1, xhat, tol );
210  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
211 #endif
212 
213  // Test3: Solve with iterative refinement.
214 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
215  // Iterative refinement not implemented in Eigen
216 #else
217  // Create random linear system
220 
221  // Create LHS through multiplication with A2
222  xhat.putScalar( ScalarTraits<STYPE>::zero() );
225 #ifdef HAVE_TEUCHOS_COMPLEX
227 #endif
228 
229  // Create a serial dense solver.
231  solver2.solveToRefinedSolution( true );
232 
233  // Pass in matrix and vectors
234  solver2.setMatrix( A2 );
235  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
236 
237  // Factor and solve with iterative refinement.
238  returnCode = solver2.factor();
239  testName = "Solve with iterative refinement: factor() random A:";
240  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
241 
242  // Non-transpose solve
243  returnCode = solver2.solve();
244  testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
245  numberFailedTests += CompareVectors( *x2, xhat, tol );
246  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
247 
248  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
249  xhat.putScalar( ScalarTraits<STYPE>::zero() );
250  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
252  returnCode = solver2.solve();
253  testName = "Solve with iterative refinement: solve() random A (TRANS):";
254  numberFailedTests += CompareVectors( *x2, xhat, tol );
255  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
256 
257 #ifdef HAVE_TEUCHOS_COMPLEX
258  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
259  xhat.putScalar( ScalarTraits<STYPE>::zero() );
260  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
262  returnCode = solver2.solve();
263  testName = "Solve with iterative refinement: solve() random A (CONJ_TRANS):";
264  numberFailedTests += CompareVectors( *x2, xhat, tol );
265  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
266 #endif
267 #endif
268 
269  // Test4: Solve with matrix equilibration.
270 
271  // Create random linear system
274 
275  // Create LHS through multiplication with A3
276  xhat.putScalar( ScalarTraits<STYPE>::zero() );
279 #ifdef HAVE_TEUCHOS_COMPLEX
281 #endif
282 
283  // Create a serial dense solver.
285  solver3.factorWithEquilibration( true );
286 
287  // Pass in matrix and vectors
288  solver3.setMatrix( A3 );
289  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
290 
291  // Factor and solve with matrix equilibration.
292  returnCode = solver3.factor();
293  testName = "Solve with matrix equilibration: factor() random A:";
294  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
295 
296  // Non-transpose solve
297  returnCode = solver3.solve();
298  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
299  numberFailedTests += CompareVectors( *x3, xhat, tol );
300  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
301 
302  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
303  xhat.putScalar( ScalarTraits<STYPE>::zero() );
304  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
306  returnCode = solver3.solve();
307  testName = "Solve with matrix equilibration: solve() random A (TRANS):";
308  numberFailedTests += CompareVectors( *x3, xhat, tol );
309  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
310 
311 #ifdef HAVE_TEUCHOS_COMPLEX
312  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
313  xhat.putScalar( ScalarTraits<STYPE>::zero() );
314  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
316  returnCode = solver3.solve();
317  testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
318  numberFailedTests += CompareVectors( *x3, xhat, tol );
319  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
320 #endif
321 
322  //
323  // If a test failed output the number of failed tests.
324  //
325  if(numberFailedTests > 0)
326  {
327  if (verbose) {
328  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
329  std::cout << "End Result: TEST FAILED" << std::endl;
330  return -1;
331  }
332  }
333  if(numberFailedTests == 0)
334  std::cout << "End Result: TEST PASSED" << std::endl;
335 
336  return 0;
337 }
338 
339 template<typename TYPE>
340 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
341 {
342  int result;
343  if(calculatedResult == expectedResult)
344  {
345  if(verbose) std::cout << testName << " successful." << std::endl;
346  result = 0;
347  }
348  else
349  {
350  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
351  result = 1;
352  }
353  return result;
354 }
355 
356 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
357 {
358  int result;
359  if(expectedResult == 0)
360  {
361  if(returnCode == 0)
362  {
363  if(verbose) std::cout << testName << " test successful." << std::endl;
364  result = 0;
365  }
366  else
367  {
368  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
369  result = 1;
370  }
371  }
372  else
373  {
374  if(returnCode != 0)
375  {
376  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
377  result = 0;
378  }
379  else
380  {
381  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
382  result = 1;
383  }
384  }
385  return result;
386 }
387 
388 template<typename TYPE>
389 TYPE GetRandom(TYPE Low, TYPE High)
390 {
391  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
392 }
393 
394 template<typename T>
395 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
396 {
397  T lowMag = Low.real();
398  T highMag = High.real();
399  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
400  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
401  return std::complex<T>( real, imag );
402 }
403 
404 template<>
405 int GetRandom(int Low, int High)
406 {
407  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
408 }
409 
410 template<>
411 double GetRandom(double Low, double High)
412 {
413  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
414 }
415 
417 {
418  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
419 
420  // Fill dense matrix with random entries.
421  for (int i=0; i<m; i++)
422  for (int j=0; j<n; j++)
423  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
424 
425  return newmat;
426 }
427 
429 {
430  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
431 
432  // Fill dense vector with random entries.
433  for (int i=0; i<n; i++)
434  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
435 
436  return newvec;
437 }
438 
439 /* Function: CompareVectors
440  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
441 */
443  const SerialDenseVector<OTYPE,STYPE>& Vector2,
445 {
446  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
447 
448  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
449  diff -= Vector2;
450 
451  MagnitudeType norm_diff = diff.normFrobenius();
452  MagnitudeType norm_v2 = Vector2.normFrobenius();
453  MagnitudeType temp = norm_diff;
454  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
455  temp /= norm_v2;
456 
457  if (temp > Tolerance)
458  return 1;
459  else
460  return 0;
461 }
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
Templated serial dense matrix class.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
Teuchos::RCP< DVector > GetRandomVector(int n)
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
int main(int argc, char *argv[])
TYPE GetRandom(TYPE, TYPE)
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 solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
Templated class for solving dense linear problems.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
int invert()
Inverts the this matrix.
std::string Teuchos_Version()
int PrintTestResults(std::string, TYPE, TYPE, bool)
int ReturnCodeCheck(std::string, int, int, bool)
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
SerialDenseMatrix< OTYPE, STYPE > DMatrix
Templated serial dense vector class.
Defines basic traits for the scalar field type.
SerialDenseVector< OTYPE, STYPE > DVector
Smart reference counting pointer class for automatic garbage collection.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
#define SCALARMAX
Reference-counted pointer class and non-member templated function implementations.
A class for solving dense linear problems.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...