Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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
72template<typename TYPE>
73int PrintTestResults(std::string, TYPE, TYPE, bool);
74
75int ReturnCodeCheck(std::string, int, int, bool);
76
80
81// Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
82template<typename TYPE>
83TYPE GetRandom(TYPE, TYPE);
84
85// Returns a random integer between the two input parameters, inclusive
86template<>
87int GetRandom(int, int);
88
89// Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
90template<>
91double GetRandom(double, double);
92
93template<typename T>
94std::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
107int 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
298template<typename TYPE>
299int 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
315int 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
347template<typename TYPE>
348TYPE GetRandom(TYPE Low, TYPE High)
349{
350 return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
351}
352
353template<typename T>
354std::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
363template<>
364int GetRandom(int Low, int High)
365{
366 return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
367}
368
369template<>
370double 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
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}
Reference-counted pointer class and non-member templated function implementations.
Defines basic traits for the scalar field type.
Non-member helper functions on the templated serial, dense matrix/vector classes.
Templated serial dense matrix class.
Templated serial dense vector class.
Templated class for solving dense linear problems.
Templated class for constructing and using Hermitian positive definite dense matrices.
Templated serial, dense, symmetric matrix class.
Smart reference counting pointer class for automatic garbage collection.
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.
This class creates and provides basic support for dense vectors of templated type as a specialization...
A class for solving dense linear problems.
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...
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int formQ()
Explicitly forms the unitary matrix Q.
A class for constructing and using Hermitian positive definite dense matrices.
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).
int factor()
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
void solveToRefinedSolution(bool flag)
Causes all solves to compute solution to best ability using iterative refinement.
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
Teuchos::RCP< SDMatrix > GetRandomSpdMatrix(int n)
TYPE GetRandom(TYPE, TYPE)
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
SerialDenseVector< OTYPE, STYPE > DVector
int PrintTestResults(std::string, TYPE, TYPE, bool)
SerialSymDenseMatrix< OTYPE, STYPE > SDMatrix
#define SCALARMAX
SerialDenseMatrix< OTYPE, STYPE > DMatrix
int ReturnCodeCheck(std::string, int, int, bool)
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Teuchos::RCP< DVector > GetRandomVector(int n)
int main()
Definition: evilMain.cpp:75
Definition: PackageA.cpp:3
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
std::string Teuchos_Version()
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
static T one()
Returns representation of one for this scalar type.