Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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
69template<typename TYPE>
70int PrintTestResults(std::string, TYPE, TYPE, bool);
71
72int ReturnCodeCheck(std::string, int, int, bool);
73
76
77// Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
78template<typename TYPE>
79TYPE GetRandom(TYPE, TYPE);
80
81// Returns a random integer between the two input parameters, inclusive
82template<>
83int GetRandom(int, int);
84
85// Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
86template<>
87double GetRandom(double, double);
88
89template<typename T>
90std::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
103int 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 ) );
174 solver1.solveWithTransposeFlag( Teuchos::TRANS );
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 ) );
184 solver1.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
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
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 ) );
252 solver2.solveWithTransposeFlag( Teuchos::TRANS );
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 ) );
262 solver2.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
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 ) );
310 solver3.solveWithTransposeFlag( Teuchos::TRANS );
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 ) );
320 solver3.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
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 );
332 solver3.solveWithTransposeFlag( Teuchos::NO_TRANS );
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
355template<typename TYPE>
356int 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
372int 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
404template<typename TYPE>
405TYPE GetRandom(TYPE Low, TYPE High)
406{
407 return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
408}
409
410template<typename T>
411std::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
420template<>
421int GetRandom(int Low, int High)
422{
423 return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
424}
425
426template<>
427double 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}
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 class for solving dense linear problems.
Templated serial dense vector class.
Concrete serial communicator subclass.
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...
TYPE GetRandom(TYPE, TYPE)
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
SerialDenseVector< OTYPE, STYPE > DVector
int PrintTestResults(std::string, TYPE, TYPE, bool)
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance, bool verbose)
#define SCALARMAX
SerialDenseMatrix< OTYPE, STYPE > DMatrix
int ReturnCodeCheck(std::string, int, int, bool)
Teuchos::RCP< DVector > GetRandomVector(int n)
int main()
Definition evilMain.cpp:75
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 T one()
Returns representation of one for this scalar type.