Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cxx_main_band_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_RCP.hpp"
49#include "Teuchos_Version.hpp"
50
55
56#define OTYPE int
57#ifdef HAVE_TEUCHOS_COMPLEX
58#define STYPE std::complex<double>
59#else
60#define STYPE double
61#endif
62
63// SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
64// random numbers in [-SCALARMAX, SCALARMAX] will be generated.
65#ifdef HAVE_TEUCHOS_COMPLEX
66#define SCALARMAX STYPE(10,0)
67#else
68#define SCALARMAX STYPE(10)
69#endif
70
71template<typename TYPE>
72int PrintTestResults(std::string, TYPE, TYPE, bool);
73
74int ReturnCodeCheck(std::string, int, int, bool);
75
79
80// Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
81template<typename TYPE>
82TYPE GetRandom(TYPE, TYPE);
83
84// Returns a random integer between the two input parameters, inclusive
85template<>
86int GetRandom(int, int);
87
88// Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
89template<>
90double GetRandom(double, double);
91
92template<typename T>
93std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
94
95// Generates random matrices and vectors using GetRandom()
96Teuchos::RCP<DMatrix> GetRandomBandedMatrix(int m, int n, int kl, int ku);
98
99// Compares the difference between two vectors using relative euclidean norms
100// Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
102 const SerialDenseVector<OTYPE,STYPE>& Vector2,
104
105// Compares the difference between two matrices using relative euclidean norms
106// Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
108 const SerialDenseMatrix<OTYPE,STYPE>& Matrix2,
110
111int main(int argc, char* argv[])
112
113{
114
115 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
116
117 int n=10, kl=2, ku=3;
118 MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
119
120 bool verbose = 0;
121 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
122
123 if (verbose)
124 std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
125
126 int numberFailedTests = 0;
127 int returnCode = 0;
128 std::string testName = "", testType = "";
129
130#ifdef HAVE_TEUCHOS_COMPLEX
131 testType = "COMPLEX";
132#else
133 testType = "REAL";
134#endif
135
136 if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL BAND DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
137
138 // Create dense matrix and vector.
141 DVector xhat(n), b(n), bt(n);
142
143 // Create a serial dense solver.
145
148
149 // Convert the dense matrix to a matrix in LAPACK banded storage
150 AB1 = Teuchos::generalToBanded( A1, kl, ku, true );
151
152 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
153 C1 = Teuchos::bandedToGeneral( AB1 );
154 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
155 numberFailedTests += CompareMatrices( *A1, *C1, tol );
156 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
157
158 // Compute the right-hand side vector using multiplication.
160 testName = "Generating right-hand side vector using A*x, where x is a random vector:";
161 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
163 testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
164 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
165
166#ifdef HAVE_TEUCHOS_COMPLEX
167 DVector bct(n);
169 testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
170 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
171#endif
172
173 // Fill the solution vector with zeros.
174 xhat.putScalar( ScalarTraits<STYPE>::zero() );
175
176 // Pass in matrix and vectors
177 solver1.setMatrix( AB1 );
178 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
179
180 // Test1: Simple factor and solve
181 returnCode = solver1.factor();
182 testName = "Simple solve: factor() random A:";
183 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
184
185 // Non-transpose solve
186 returnCode = solver1.solve();
187 testName = "Simple solve: solve() random A (NO_TRANS):";
188 numberFailedTests += CompareVectors( *x1, xhat, tol );
189 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
190
191 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
192 xhat.putScalar( ScalarTraits<STYPE>::zero() );
193 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
195 returnCode = solver1.solve();
196 testName = "Simple solve: solve() random A (TRANS):";
197 numberFailedTests += CompareVectors( *x1, xhat, tol );
198 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
199
200#ifdef HAVE_TEUCHOS_COMPLEX
201 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
202 xhat.putScalar( ScalarTraits<STYPE>::zero() );
203 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
205 returnCode = solver1.solve();
206 testName = "Simple solve: solve() random A (CONJ_TRANS):";
207 numberFailedTests += CompareVectors( *x1, xhat, tol );
208 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
209#endif
210
211 // Test2: Solve with iterative refinement.
212#ifdef HAVE_TEUCHOSNUMERICS_EIGEN
213 // Iterative refinement not implemented in Eigen
214#else
215 // Create random linear system
218
219 // Create LHS through multiplication with A2
220 xhat.putScalar( ScalarTraits<STYPE>::zero() );
223#ifdef HAVE_TEUCHOS_COMPLEX
225#endif
226
227 // Create a serial dense solver.
229 solver2.solveToRefinedSolution( true );
230
233 // Convert the dense matrix to a matrix in LAPACK banded storage
234 AB2 = Teuchos::generalToBanded( A2, kl, ku, true );
235
236 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
237 C2 = Teuchos::bandedToGeneral( AB2 );
238 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
239 numberFailedTests += CompareMatrices( *A2, *C2, tol );
240 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
241
242 // Pass in matrix and vectors
243 solver2.setMatrix( AB2 );
244 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
245
246 // Factor and solve with iterative refinement.
247 returnCode = solver2.factor();
248 testName = "Solve with iterative refinement: factor() random A:";
249 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
250
251 // Non-transpose solve
252 returnCode = solver2.solve();
253 testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
254 numberFailedTests += CompareVectors( *x2, xhat, tol );
255 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
256
257 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
258 xhat.putScalar( ScalarTraits<STYPE>::zero() );
259 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
261 returnCode = solver2.solve();
262 testName = "Solve with iterative refinement: solve() random A (TRANS):";
263 numberFailedTests += CompareVectors( *x2, xhat, tol );
264 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
265
266#ifdef HAVE_TEUCHOS_COMPLEX
267 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
268 xhat.putScalar( ScalarTraits<STYPE>::zero() );
269 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
271 returnCode = solver2.solve();
272 testName = "Solve with iterative refinement: solve() random A (CONJ_TRANS):";
273 numberFailedTests += CompareVectors( *x2, xhat, tol );
274 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
275#endif
276#endif
277
278 // Test3: Solve with matrix equilibration.
279
280 // Create random linear system
283
284 // Create LHS through multiplication with A3
285 xhat.putScalar( ScalarTraits<STYPE>::zero() );
288#ifdef HAVE_TEUCHOS_COMPLEX
290#endif
291
292 // Create a serial dense solver.
294 solver3.factorWithEquilibration( true );
295
298 // Convert the dense matrix to a matrix in LAPACK banded storage
299 AB3 = Teuchos::generalToBanded( A3, kl, ku, true );
300
301 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
302 C3 = Teuchos::bandedToGeneral( AB3 );
303 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
304 numberFailedTests += CompareMatrices( *A3, *C3, tol );
305 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
306
307 // Pass in matrix and vectors
308 solver3.setMatrix( AB3 );
309 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
310
311 Teuchos::RCP<BDMatrix> AB3bak = Teuchos::rcp( new BDMatrix( *AB3 ) );
313
314 // Factor and solve with matrix equilibration.
315 returnCode = solver3.factor();
316 testName = "Solve with matrix equilibration: factor() random A:";
317 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
318
319 // Non-transpose solve
320 returnCode = solver3.solve();
321 testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
322 numberFailedTests += CompareVectors( *x3, xhat, tol );
323 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
324
325 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
326 xhat.putScalar( ScalarTraits<STYPE>::zero() );
327 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
329 returnCode = solver3.solve();
330 testName = "Solve with matrix equilibration: solve() random A (TRANS):";
331 numberFailedTests += CompareVectors( *x3, xhat, tol );
332 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
333
334#ifdef HAVE_TEUCHOS_COMPLEX
335 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
336 xhat.putScalar( ScalarTraits<STYPE>::zero() );
337 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
339 returnCode = solver3.solve();
340 testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
341 numberFailedTests += CompareVectors( *x3, xhat, tol );
342 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
343#endif
344
345 // Non-tranpose solve where factor is not called
346 xhat.putScalar( ScalarTraits<STYPE>::zero() );
347 solver3.setMatrix( AB3bak );
348 solver3.setVectors( Teuchos::rcp( &xhat, false ), b3bak );
350 returnCode = solver3.solve();
351 testName = "Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
352 numberFailedTests += CompareVectors( *x3, xhat, tol );
353 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
354
355 //
356 // If a test failed output the number of failed tests.
357 //
358 if(numberFailedTests > 0)
359 {
360 if (verbose) {
361 std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
362 std::cout << "End Result: TEST FAILED" << std::endl;
363 return -1;
364 }
365 }
366 if(numberFailedTests == 0)
367 std::cout << "End Result: TEST PASSED!" << std::endl;
368
369 return 0;
370}
371
372template<typename TYPE>
373int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
374{
375 int result;
376 if(calculatedResult == expectedResult)
377 {
378 if(verbose) std::cout << testName << " successful." << std::endl;
379 result = 0;
380 }
381 else
382 {
383 if(verbose) std::cout << testName << " unsuccessful." << std::endl;
384 result = 1;
385 }
386 return result;
387}
388
389int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
390{
391 int result;
392 if(expectedResult == 0)
393 {
394 if(returnCode == 0)
395 {
396 if(verbose) std::cout << testName << " test successful." << std::endl;
397 result = 0;
398 }
399 else
400 {
401 if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
402 result = 1;
403 }
404 }
405 else
406 {
407 if(returnCode != 0)
408 {
409 if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
410 result = 0;
411 }
412 else
413 {
414 if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
415 result = 1;
416 }
417 }
418 return result;
419}
420
421template<typename TYPE>
422TYPE GetRandom(TYPE Low, TYPE High)
423{
424 return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
425}
426
427template<typename T>
428std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
429{
430 T lowMag = Low.real();
431 T highMag = High.real();
432 T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
433 T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
434 return std::complex<T>( real, imag );
435}
436
437template<>
438int GetRandom(int Low, int High)
439{
440 return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
441}
442
443template<>
444double GetRandom(double Low, double High)
445{
446 return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
447}
448
449Teuchos::RCP<DMatrix> GetRandomBandedMatrix(int m, int n, int kl, int ku)
450{
451 Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
452
453 // Fill dense matrix with random entries.
454 for (int i=0; i<m; i++)
455 for (int j=0; j<n; j++)
456 if (j>= i-kl && j<=i+ku)
457 (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
458
459 return newmat;
460}
461
463{
464 Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
465
466 // Fill dense vector with random entries.
467 for (int i=0; i<n; i++)
468 (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
469
470 return newvec;
471}
472
473/* Function: CompareVectors
474 Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
475*/
477 const SerialDenseVector<OTYPE,STYPE>& Vector2,
479{
480 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
481
482 SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
483 diff -= Vector2;
484
485 MagnitudeType norm_diff = diff.normFrobenius();
486 MagnitudeType norm_v2 = Vector2.normFrobenius();
487 MagnitudeType temp = norm_diff;
488 if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
489 temp /= norm_v2;
490
491 if (temp > Tolerance)
492 return 1;
493 else
494 return 0;
495}
496
497/* Function: CompareVectors
498 Purpose: Compares the difference between two matrices using relative euclidean-norms, i.e. ||m_1-m_2||_\inf/||m_2||_\inf
499*/
501 const SerialDenseMatrix<OTYPE,STYPE>& Matrix2,
503{
504 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
505
506 SerialDenseMatrix<OTYPE,STYPE> diff( Matrix1 );
507 diff -= Matrix2;
508
509 MagnitudeType norm_diff = diff.normInf();
510 MagnitudeType norm_m2 = Matrix2.normInf();
511 MagnitudeType temp = norm_diff;
512 if (norm_m2 != ScalarTraits<MagnitudeType>::zero())
513 temp /= norm_m2;
514
515 if (temp > Tolerance)
516 return 1;
517 else
518 return 0;
519}
Reference-counted pointer class and non-member templated function implementations.
Defines basic traits for the scalar field type.
Templated serial dense matrix class.
Non-member helper functions on the templated serial, dense matrix/vector classes.
Templated serial dense matrix class.
Templated serial dense vector class.
Smart reference counting pointer class for automatic garbage collection.
This class creates and provides basic support for banded dense matrices of templated type.
A class for representing and solving banded dense linear systems.
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 setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
void solveWithTransposeFlag(Teuchos::ETransp trans)
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
This class creates and provides basic support for dense rectangular matrix of templated type.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
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.
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...
Teuchos::RCP< DMatrix > GetRandomBandedMatrix(int m, int n, int kl, int ku)
TYPE GetRandom(TYPE, TYPE)
SerialDenseVector< OTYPE, STYPE > DVector
int PrintTestResults(std::string, TYPE, TYPE, bool)
#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 CompareMatrices(const SerialDenseMatrix< OTYPE, STYPE > &Matrix1, const SerialDenseMatrix< OTYPE, STYPE > &Matrix2, ScalarTraits< STYPE >::magnitudeType Tolerance)
SerialBandDenseMatrix< OTYPE, STYPE > BDMatrix
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.