Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  bool verbose);
102 
103 int main(int argc, char* argv[])
104 {
105  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
106 
107  int n=10, m=8;
108  (void) m; // forestall "unused variable" compiler warning
109  MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
110 
111  bool verbose = 0;
112  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
113 
114  if (verbose)
115  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
116 
117  int numberFailedTests = 0;
118  int returnCode = 0;
119  std::string testName = "", testType = "";
120 
121 #ifdef HAVE_TEUCHOS_COMPLEX
122  testType = "COMPLEX";
123 #else
124  testType = "REAL";
125 #endif
126 
127  if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
128 
129  // Create dense matrix and vector.
132  DVector xhat(n), b(n), bt(n);
133 
134  // Compute the right-hand side vector using multiplication.
136  testName = "Generating right-hand side vector using A*x, where x is a random vector:";
137  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
138 
140  testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
141  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
142 
143 #ifdef HAVE_TEUCHOS_COMPLEX
144  DVector bct(n);
146  testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
147  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
148 #endif
149 
150  // Fill the solution vector with zeros.
151  xhat.putScalar( ScalarTraits<STYPE>::zero() );
152 
153  // Create a serial dense solver.
155 
156  // Pass in matrix and vectors
157  solver1.setMatrix( A1 );
158  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
159 
160  // Test1: Simple factor and solve
161  returnCode = solver1.factor();
162  testName = "Simple solve: factor() random A:";
163  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
164 
165  // Non-transpose solve
166  returnCode = solver1.solve();
167  testName = "Simple solve: solve() random A (NO_TRANS):";
168  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
169  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
170 
171  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
172  xhat.putScalar( ScalarTraits<STYPE>::zero() );
173  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
175  returnCode = solver1.solve();
176  testName = "Simple solve: solve() random A (TRANS):";
177  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
178  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
179 
180 #ifdef HAVE_TEUCHOS_COMPLEX
181  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
182  xhat.putScalar( ScalarTraits<STYPE>::zero() );
183  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
185  returnCode = solver1.solve();
186  testName = "Simple solve: solve() random A (CONJ_TRANS):";
187  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
188  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
189 #endif
190 
191  // Test2: Invert the matrix, inverse should be in A.
192  returnCode = solver1.invert();
193  testName = "Simple solve: invert() random A:";
194  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
195 
196  // Compute the solution vector using multiplication and the inverse.
198  testName = "Computing solution using inverted random A (NO_TRANS):";
199  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
200  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
201 
202  returnCode = xhat.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bt, ScalarTraits<STYPE>::zero());
203  testName = "Computing solution using inverted random A (TRANS):";
204  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
205  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
206 
207 #ifdef HAVE_TEUCHOS_COMPLEX
208  returnCode = xhat.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bct, ScalarTraits<STYPE>::zero());
209  testName = "Computing solution using inverted random A (CONJ_TRANS):";
210  numberFailedTests += CompareVectors( *x1, xhat, tol, verbose );
211  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
212 #endif
213 
214  // Test3: Solve with iterative refinement.
215 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
216  // Iterative refinement not implemented in Eigen
217 #else
218  // Create random linear system
221 
222  // Create LHS through multiplication with A2
223  xhat.putScalar( ScalarTraits<STYPE>::zero() );
226 #ifdef HAVE_TEUCHOS_COMPLEX
228 #endif
229 
230  // Create a serial dense solver.
232  solver2.solveToRefinedSolution( true );
233 
234  // Pass in matrix and vectors
235  solver2.setMatrix( A2 );
236  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
237 
238  // Factor and solve with iterative refinement.
239  returnCode = solver2.factor();
240  testName = "Solve with iterative refinement: factor() random A:";
241  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
242 
243  // Non-transpose solve
244  returnCode = solver2.solve();
245  testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
246  numberFailedTests += CompareVectors( *x2, xhat, tol, verbose );
247  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
248 
249  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
250  xhat.putScalar( ScalarTraits<STYPE>::zero() );
251  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
253  returnCode = solver2.solve();
254  testName = "Solve with iterative refinement: solve() random A (TRANS):";
255  numberFailedTests += CompareVectors( *x2, xhat, tol, verbose );
256  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
257 
258 #ifdef HAVE_TEUCHOS_COMPLEX
259  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
260  xhat.putScalar( ScalarTraits<STYPE>::zero() );
261  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
263  returnCode = solver2.solve();
264  testName = "Solve with iterative refinement: solve() random A (CONJ_TRANS):";
265  numberFailedTests += CompareVectors( *x2, xhat, tol, verbose );
266  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
267 #endif
268 #endif
269 
270  // Test4: Solve with matrix equilibration.
271 
272  // Create random linear system
275 
276  // Create LHS through multiplication with A3
277  xhat.putScalar( ScalarTraits<STYPE>::zero() );
280 #ifdef HAVE_TEUCHOS_COMPLEX
282 #endif
283 
284  // Save backups for multiple solves.
287 
288  // Create a serial dense solver.
290  solver3.factorWithEquilibration( true );
291 
292  // Pass in matrix and vectors
293  solver3.setMatrix( A3 );
294  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
295 
296  // Factor and solve with matrix equilibration.
297  returnCode = solver3.factor();
298  testName = "Solve with matrix equilibration: factor() random A:";
299  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
300 
301  // Non-transpose solve
302  returnCode = solver3.solve();
303  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
304  numberFailedTests += CompareVectors( *x3, xhat, tol, verbose );
305  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
306 
307  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
308  xhat.putScalar( ScalarTraits<STYPE>::zero() );
309  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
311  returnCode = solver3.solve();
312  testName = "Solve with matrix equilibration: solve() random A (TRANS):";
313  numberFailedTests += CompareVectors( *x3, xhat, tol, verbose );
314  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
315 
316 #ifdef HAVE_TEUCHOS_COMPLEX
317  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
318  xhat.putScalar( ScalarTraits<STYPE>::zero() );
319  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
321  returnCode = solver3.solve();
322  testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
323  numberFailedTests += CompareVectors( *x3, xhat, tol, verbose );
324  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
325 #endif
326 
327  // Factor and solve with matrix equilibration, only call solve not factor.
328  // Use copy of A3 and b, they were overwritten in last factor() call.
329  xhat.putScalar( ScalarTraits<STYPE>::zero() );
330  solver3.setMatrix( A3bak );
331  solver3.setVectors( Teuchos::rcp( &xhat, false ), b3bak );
333  returnCode = solver3.solve();
334  testName = "Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
335  numberFailedTests += CompareVectors( *x3, xhat, tol, verbose );
336  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
337 
338  //
339  // If a test failed output the number of failed tests.
340  //
341  if(numberFailedTests > 0)
342  {
343  if (verbose) {
344  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
345  std::cout << "End Result: TEST FAILED" << std::endl;
346  return -1;
347  }
348  }
349  if(numberFailedTests == 0)
350  std::cout << "End Result: TEST PASSED" << std::endl;
351 
352  return 0;
353 }
354 
355 template<typename TYPE>
356 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
357 {
358  int result;
359  if(calculatedResult == expectedResult)
360  {
361  if(verbose) std::cout << testName << " successful." << std::endl;
362  result = 0;
363  }
364  else
365  {
366  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
367  result = 1;
368  }
369  return result;
370 }
371 
372 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
373 {
374  int result;
375  if(expectedResult == 0)
376  {
377  if(returnCode == 0)
378  {
379  if(verbose) std::cout << testName << " test successful." << std::endl;
380  result = 0;
381  }
382  else
383  {
384  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
385  result = 1;
386  }
387  }
388  else
389  {
390  if(returnCode != 0)
391  {
392  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
393  result = 0;
394  }
395  else
396  {
397  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
398  result = 1;
399  }
400  }
401  return result;
402 }
403 
404 template<typename TYPE>
405 TYPE GetRandom(TYPE Low, TYPE High)
406 {
407  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
408 }
409 
410 template<typename T>
411 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
412 {
413  T lowMag = Low.real();
414  T highMag = High.real();
415  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
416  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
417  return std::complex<T>( real, imag );
418 }
419 
420 template<>
421 int GetRandom(int Low, int High)
422 {
423  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
424 }
425 
426 template<>
427 double GetRandom(double Low, double High)
428 {
429  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
430 }
431 
433 {
434  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
435 
436  // Fill dense matrix with random entries.
437  for (int i=0; i<m; i++)
438  for (int j=0; j<n; j++)
439  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
440 
441  return newmat;
442 }
443 
445 {
446  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
447 
448  // Fill dense vector with random entries.
449  for (int i=0; i<n; i++)
450  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
451 
452  return newvec;
453 }
454 
455 /* Function: CompareVectors
456  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
457 */
459  const SerialDenseVector<OTYPE,STYPE>& Vector2,
461  bool verbose)
462 {
463  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
464 
465  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
466  diff -= Vector2;
467 
468  MagnitudeType norm_diff = diff.normFrobenius();
469  MagnitudeType norm_v2 = Vector2.normFrobenius();
470  MagnitudeType temp = norm_diff;
471  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
472  temp /= norm_v2;
473 
474  if (temp > Tolerance)
475  {
476  if (verbose)
477  std::cout << "COMPARISON FAILED : ";
478  return 1;
479  }
480  else
481  return 0;
482 }
int PrintTestResults(std::string, TYPE, TYPE, bool)
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
Templated serial dense matrix class.
Teuchos::SerialDenseVector< int, std::complex< Real > > DVector
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.
This class creates and provides basic support for dense vectors of templated type as a specialization...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
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.
TYPE GetRandom(TYPE, TYPE)
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.
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
bool CompareVectors(TYPE1 *Vector1, TYPE2 *Vector2, OType Size, TYPE2 Tolerance)
Teuchos::SerialDenseMatrix< int, std::complex< Real > > DMatrix
std::string Teuchos_Version()
int ReturnCodeCheck(std::string, int, int, bool)
int main(int argc, char *argv[])
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).
Templated serial dense vector class.
Defines basic traits for the scalar field type.
Teuchos::RCP< DVector > GetRandomVector(int n)
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...