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_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 
47 #include "Teuchos_ScalarTraits.hpp"
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 
71 template<typename TYPE>
72 int PrintTestResults(std::string, TYPE, TYPE, bool);
73 
74 int ReturnCodeCheck(std::string, int, int, bool);
75 
79 
80 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
81 template<typename TYPE>
82 TYPE GetRandom(TYPE, TYPE);
83 
84 // Returns a random integer between the two input parameters, inclusive
85 template<>
86 int GetRandom(int, int);
87 
88 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
89 template<>
90 double GetRandom(double, double);
91 
92 template<typename T>
93 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
94 
95 // Generates random matrices and vectors using GetRandom()
96 Teuchos::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 
111 int 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 
372 template<typename TYPE>
373 int 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 
389 int 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 
421 template<typename TYPE>
422 TYPE GetRandom(TYPE Low, TYPE High)
423 {
424  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
425 }
426 
427 template<typename T>
428 std::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 
437 template<>
438 int GetRandom(int Low, int High)
439 {
440  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
441 }
442 
443 template<>
444 double 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 
449 Teuchos::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 }
int PrintTestResults(std::string, TYPE, TYPE, bool)
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
Templated serial dense matrix class.
void solveToRefinedSolution(bool flag)
Set whether or not to use iterative refinement to improve solutions to linear systems.
int setMatrix(const RCP< SerialBandDenseMatrix< OrdinalType, ScalarType > > &AB)
Sets the pointer for coefficient matrix.
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.
#define SCALARMAX
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).
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.
Templated serial dense matrix class.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
TYPE GetRandom(TYPE, TYPE)
This class creates and provides basic support for banded dense matrices of templated type...
A class for representing and solving banded dense linear systems.
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)
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
int main(int argc, char *argv[])
bool CompareMatrices(TYPE1 *Matrix1, TYPE2 *Matrix2, OType Rows, OType Columns, OType LDM, TYPE2 Tolerance)
Teuchos::RCP< DMatrix > GetRandomBandedMatrix(int m, int n, int kl, int ku)
Teuchos::SerialBandDenseMatrix< OTYPE, STYPE > BDMatrix
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GBTRF.
Templated serial dense vector class.
Defines basic traits for the scalar field type.
void solveWithTransposeFlag(Teuchos::ETransp trans)
Teuchos::RCP< DVector > GetRandomVector(int n)
Smart reference counting pointer class for automatic garbage collection.
Reference-counted pointer class and non-member templated function implementations.
This class creates and provides basic support for dense rectangular matrix of templated type...