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_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  Teuchos::RCP<SDMatrix> A3bak = Teuchos::rcp( new SDMatrix( *A3 ) );
248 
249  // Create a serial dense solver.
251  solver3.factorWithEquilibration( true );
252 
253  // Pass in matrix and vectors
254  solver3.setMatrix( A3 );
255  solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
256 
257  // Factor and solve with matrix equilibration.
258  returnCode = solver3.factor();
259  testName = "Solve with matrix equilibration: factor() random A:";
260  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
261 
262  // Non-transpose solve
263  returnCode = solver3.solve();
264  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
265  failedComparison = CompareVectors( *x3, xhat, tol );
266  if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
267  numberFailedTests += failedComparison;
268  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
269 
270  // Factor and solve with matrix equilibration, where only solve is called.
271  solver3.setMatrix( A3bak );
272  solver3.setVectors( Teuchos::rcp( &xhat, false ), b3bak );
273  returnCode = solver3.solve();
274  testName = "Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
275  failedComparison = CompareVectors( *x3, xhat, tol );
276  if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
277  numberFailedTests += failedComparison;
278  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
279 
280  //
281  // If a test failed output the number of failed tests.
282  //
283  if(numberFailedTests > 0)
284  {
285  if (verbose) {
286  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
287  std::cout << "End Result: TEST FAILED" << std::endl;
288  return -1;
289  }
290  }
291  if(numberFailedTests == 0)
292  std::cout << "End Result: TEST PASSED" << std::endl;
293 
294  return 0;
295 
296 }
297 
298 template<typename TYPE>
299 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
300 {
301  int result;
302  if(calculatedResult == expectedResult)
303  {
304  if(verbose) std::cout << testName << " successful." << std::endl;
305  result = 0;
306  }
307  else
308  {
309  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
310  result = 1;
311  }
312  return result;
313 }
314 
315 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
316 {
317  int result;
318  if(expectedResult == 0)
319  {
320  if(returnCode == 0)
321  {
322  if(verbose) std::cout << testName << " test successful." << std::endl;
323  result = 0;
324  }
325  else
326  {
327  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
328  result = 1;
329  }
330  }
331  else
332  {
333  if(returnCode != 0)
334  {
335  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
336  result = 0;
337  }
338  else
339  {
340  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
341  result = 1;
342  }
343  }
344  return result;
345 }
346 
347 template<typename TYPE>
348 TYPE GetRandom(TYPE Low, TYPE High)
349 {
350  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
351 }
352 
353 template<typename T>
354 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
355 {
356  T lowMag = Low.real();
357  T highMag = High.real();
358  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
359  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
360  return std::complex<T>( real, imag );
361 }
362 
363 template<>
364 int GetRandom(int Low, int High)
365 {
366  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
367 }
368 
369 template<>
370 double GetRandom(double Low, double High)
371 {
372  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
373 }
374 
376 {
377 
379  Teuchos::RCP<DMatrix> tmp = Teuchos::rcp( new DMatrix(n,n) );
380  Teuchos::RCP<DMatrix> tmp2 = Teuchos::rcp( new DMatrix(n,n) );
383 
384  // QR factor a random matrix and get the orthogonal Q
387  solver.setMatrix( A );
388  solver.factor();
389  solver.formQ();
390  Q = solver.getQ();
391 
392  // Get n random positive eigenvalues and put them in a diagonal matrix D
394  for (int i=0; i<n; ++i) {
397  }
398 
399  // Form the spd matrix Q' D Q
402 
403  mat->setUpper();
404  for (int j=0; j<n; ++j) {
405  for (int i=0; i<=j; ++i) {
406  (*mat)(i,j) = (*tmp2)(i,j);
407  }
408  }
409 
410  return mat;
411 }
412 
414 {
415  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
416 
417  // Fill dense matrix with random entries.
418  for (int i=0; i<m; i++)
419  for (int j=0; j<n; j++)
420  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
421 
422  return newmat;
423 }
424 
426 {
427  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
428 
429  // Fill dense vector with random entries.
430  for (int i=0; i<n; i++)
431  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
432 
433  return newvec;
434 }
435 
436 /* Function: CompareVectors
437  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
438 */
440  const SerialDenseVector<OTYPE,STYPE>& Vector2,
442 {
443  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
444 
445  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
446  diff -= Vector2;
447 
448  MagnitudeType norm_diff = diff.normFrobenius();
449  MagnitudeType norm_v2 = Vector2.normFrobenius();
450  MagnitudeType temp = norm_diff;
451  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
452  temp /= norm_v2;
453 
454  if (temp > Tolerance)
455  return 1;
456  else
457  return 0;
458 }
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.
int PrintTestResults(std::string, TYPE, TYPE, bool)
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.
Teuchos::SerialDenseVector< int, std::complex< Real > > DVector
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...
Templated class for solving dense linear problems.
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.
#define SCALARMAX
A class for solving dense linear problems.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
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.
TYPE GetRandom(TYPE, TYPE)
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).
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)
Templated serial, dense, symmetric matrix class.
int main(int argc, char *argv[])
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.
Teuchos::RCP< DVector > GetRandomVector(int n)
Smart reference counting pointer class for automatic garbage collection.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
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...