Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
cxx_main_spd_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 
48 #include "Teuchos_ScalarTraits.hpp"
49 #include "Teuchos_RCP.hpp"
50 #include "Teuchos_Version.hpp"
51 
56 
57 #define OTYPE int
58 #ifdef HAVE_TEUCHOS_COMPLEX
59 #define STYPE std::complex<double>
60 #else
61 #define STYPE double
62 #endif
63 
64 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
65 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
66 #ifdef HAVE_TEUCHOS_COMPLEX
67 #define SCALARMAX STYPE(10,0)
68 #else
69 #define SCALARMAX STYPE(10)
70 #endif
71 
72 template<typename TYPE>
73 int PrintTestResults(std::string, TYPE, TYPE, bool);
74 
75 int ReturnCodeCheck(std::string, int, int, bool);
76 
80 
81 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
82 template<typename TYPE>
83 TYPE GetRandom(TYPE, TYPE);
84 
85 // Returns a random integer between the two input parameters, inclusive
86 template<>
87 int GetRandom(int, int);
88 
89 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
90 template<>
91 double GetRandom(double, double);
92 
93 template<typename T>
94 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
95 
96 // Generates random matrices and vectors using GetRandom()
100 
101 // Compares the difference between two vectors using relative euclidean norms
102 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
104  const SerialDenseVector<OTYPE,STYPE>& Vector2,
106 
107 int main(int argc, char* argv[])
108 {
109  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
110 
111  int n=10;
112  MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
113 
114  bool verbose = 0;
115  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
116 
117  if (verbose)
118  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
119 
120  int numberFailedTests = 0, failedComparison = 0;
121  int returnCode = 0;
122  std::string testName = "", testType = "";
123 
124 #ifdef HAVE_TEUCHOS_COMPLEX
125  testType = "COMPLEX";
126 #else
127  testType = "REAL";
128 #endif
129 
130  if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL SPD DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
131 
132  // Create dense matrix and vector.
135  DVector xhat(n), b(n), bt(n);
136 
137  Teuchos::RCP<DMatrix> A1full = Teuchos::rcp( new DMatrix(n,n) );
138  for (int j=0; j<n; ++j) {
139  for (int i=0; i<=j; ++i) {
140  (*A1full)(i,j) = (*A1)(i,j);
141  }
142  }
143  for (int j=0; j<n; ++j) {
144  for (int i=j+1; i<n; ++i) {
145  (*A1full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A1)(i,j) );
146  }
147  }
148 
149  // Compute the right-hand side vector using multiplication.
150  returnCode = b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1full, *x1, ScalarTraits<STYPE>::zero());
151  testName = "Generating right-hand side vector using A*x, where x is a random vector:";
152  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
153 
154  // Fill the solution vector with zeros.
155  xhat.putScalar( ScalarTraits<STYPE>::zero() );
156 
157  // Create a serial spd dense solver.
159 
160  // Pass in matrix and vectors
161  solver1.setMatrix( A1 );
162  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
163 
164  // Test1: Simple factor and solve
165  returnCode = solver1.factor();
166  testName = "Simple solve: factor() random A:";
167  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
168 
169  // Non-transpose solve
170  returnCode = solver1.solve();
171  testName = "Simple solve: solve() random A (NO_TRANS):";
172  failedComparison = CompareVectors( *x1, xhat, tol );
173  if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
174  numberFailedTests += failedComparison;
175  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
176 
177  // Test2: Solve with iterative refinement.
178 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
179  // Iterative refinement not implemented in Eigen
180 #else
181  // Create random linear system
184 
185  Teuchos::RCP<DMatrix> A2full = Teuchos::rcp( new DMatrix(n,n) );
186  for (int j=0; j<n; ++j) {
187  for (int i=0; i<=j; ++i) {
188  (*A2full)(i,j) = (*A2)(i,j);
189  }
190  }
191  for (int j=0; j<n; ++j) {
192  for (int i=j+1; i<n; ++i) {
193  (*A2full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A2)(i,j) );
194  }
195  }
196 
197  // Create LHS through multiplication with A2
198  xhat.putScalar( ScalarTraits<STYPE>::zero() );
200 
201  // Create a serial dense solver.
203  solver2.solveToRefinedSolution( true );
204 
205  // Pass in matrix and vectors
206  solver2.setMatrix( A2 );
207  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
208 
209  // Factor and solve with iterative refinement.
210  returnCode = solver2.factor();
211  testName = "Solve with iterative refinement: factor() random A:";
212  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
213 
214  // Non-transpose solve
215  returnCode = solver2.solve();
216  testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
217  failedComparison = CompareVectors( *x2, xhat, tol );
218  if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
219  numberFailedTests += failedComparison;
220  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
221 
222 #endif
223 
224  // Test3: Solve with matrix equilibration.
225 
226  // Create random linear system
229 
230  Teuchos::RCP<DMatrix> A3full = Teuchos::rcp( new DMatrix(n,n) );
231  for (int j=0; j<n; ++j) {
232  for (int i=0; i<=j; ++i) {
233  (*A3full)(i,j) = (*A3)(i,j);
234  }
235  }
236  for (int j=0; j<n; ++j) {
237  for (int i=j+1; i<n; ++i) {
238  (*A3full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A3)(i,j) );
239  }
240  }
241 
242  // Create LHS through multiplication with A3
243  xhat.putScalar( ScalarTraits<STYPE>::zero() );
245 
246  // Create a serial dense solver.
248  solver3.factorWithEquilibration( true );
249 
250  // Pass in matrix and vectors
251  solver3.setMatrix( A3 );
252  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
253 
254  // Factor and solve with matrix equilibration.
255  returnCode = solver3.factor();
256  testName = "Solve with matrix equilibration: factor() random A:";
257  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
258 
259  // Non-transpose solve
260  returnCode = solver3.solve();
261  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
262  failedComparison = CompareVectors( *x3, xhat, tol );
263  if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
264  numberFailedTests += failedComparison;
265  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
266 
267  //
268  // If a test failed output the number of failed tests.
269  //
270  if(numberFailedTests > 0)
271  {
272  if (verbose) {
273  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
274  std::cout << "End Result: TEST FAILED" << std::endl;
275  return -1;
276  }
277  }
278  if(numberFailedTests == 0)
279  std::cout << "End Result: TEST PASSED" << std::endl;
280 
281  return 0;
282 
283 }
284 
285 template<typename TYPE>
286 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
287 {
288  int result;
289  if(calculatedResult == expectedResult)
290  {
291  if(verbose) std::cout << testName << " successful." << std::endl;
292  result = 0;
293  }
294  else
295  {
296  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
297  result = 1;
298  }
299  return result;
300 }
301 
302 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
303 {
304  int result;
305  if(expectedResult == 0)
306  {
307  if(returnCode == 0)
308  {
309  if(verbose) std::cout << testName << " test successful." << std::endl;
310  result = 0;
311  }
312  else
313  {
314  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
315  result = 1;
316  }
317  }
318  else
319  {
320  if(returnCode != 0)
321  {
322  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
323  result = 0;
324  }
325  else
326  {
327  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
328  result = 1;
329  }
330  }
331  return result;
332 }
333 
334 template<typename TYPE>
335 TYPE GetRandom(TYPE Low, TYPE High)
336 {
337  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
338 }
339 
340 template<typename T>
341 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
342 {
343  T lowMag = Low.real();
344  T highMag = High.real();
345  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
346  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
347  return std::complex<T>( real, imag );
348 }
349 
350 template<>
351 int GetRandom(int Low, int High)
352 {
353  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
354 }
355 
356 template<>
357 double GetRandom(double Low, double High)
358 {
359  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
360 }
361 
363 {
364 
370 
371  // QR factor a random matrix and get the orthogonal Q
374  solver.setMatrix( A );
375  solver.factor();
376  solver.formQ();
377  Q = solver.getQ();
378 
379  // Get n random positive eigenvalues and put them in a diagonal matrix D
381  for (int i=0; i<n; ++i) {
384  }
385 
386  // Form the spd matrix Q' D Q
389 
390  mat->setUpper();
391  for (int j=0; j<n; ++j) {
392  for (int i=0; i<=j; ++i) {
393  (*mat)(i,j) = (*tmp2)(i,j);
394  }
395  }
396 
397  return mat;
398 }
399 
401 {
402  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
403 
404  // Fill dense matrix with random entries.
405  for (int i=0; i<m; i++)
406  for (int j=0; j<n; j++)
407  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
408 
409  return newmat;
410 }
411 
413 {
414  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
415 
416  // Fill dense vector with random entries.
417  for (int i=0; i<n; i++)
418  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
419 
420  return newvec;
421 }
422 
423 /* Function: CompareVectors
424  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
425 */
427  const SerialDenseVector<OTYPE,STYPE>& Vector2,
429 {
430  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
431 
432  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
433  diff -= Vector2;
434 
435  MagnitudeType norm_diff = diff.normFrobenius();
436  MagnitudeType norm_v2 = Vector2.normFrobenius();
437  MagnitudeType temp = norm_diff;
438  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
439  temp /= norm_v2;
440 
441  if (temp > Tolerance)
442  return 1;
443  else
444  return 0;
445 }
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
SerialSymDenseMatrix< OTYPE, STYPE > SDMatrix
Teuchos::RCP< SDMatrix > GetRandomSpdMatrix(int n)
A class for constructing and using Hermitian positive definite dense matrices.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Templated serial dense matrix class.
int formQ()
Explicitly forms the unitary matrix Q.
int ReturnCodeCheck(std::string, int, int, bool)
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int main(int argc, char *argv[])
Teuchos::RCP< DVector > GetRandomVector(int n)
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.
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Templated class for solving dense linear problems.
Definition: PackageA.cpp:3
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.
#define SCALARMAX
A class for solving dense linear problems.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
int PrintTestResults(std::string, TYPE, TYPE, bool)
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
TYPE GetRandom(TYPE, TYPE)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
int setMatrix(const RCP< SerialSymDenseMatrix< OrdinalType, ScalarType > > &A_in)
Sets the pointers for coefficient matrix.
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
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).
std::string Teuchos_Version()
SerialDenseVector< OTYPE, STYPE > DVector
Templated serial, dense, symmetric matrix class.
Templated class for constructing and using Hermitian positive definite dense matrices.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
Templated serial dense vector class.
Defines basic traits for the scalar field type.
Smart reference counting pointer class for automatic garbage collection.
SerialDenseMatrix< OTYPE, STYPE > DMatrix
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
Reference-counted pointer class and non-member templated function implementations.
This class creates and provides basic support for dense rectangular matrix of templated type...