Teuchos Package Browser (Single Doxygen Collection) Version of the Day
cxx_main_solver.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "Teuchos_SerialDenseMatrix.hpp"
00043 #include "Teuchos_SerialDenseVector.hpp"
00044 #include "Teuchos_SerialDenseHelpers.hpp"
00045 #include "Teuchos_SerialDenseSolver.hpp"
00046 #include "Teuchos_ScalarTraits.hpp"
00047 #include "Teuchos_RCP.hpp"
00048 #include "Teuchos_Version.hpp"
00049 
00050 using Teuchos::ScalarTraits;
00051 using Teuchos::SerialDenseMatrix;
00052 using Teuchos::SerialDenseVector;
00053 
00054 #define OTYPE int
00055 #ifdef HAVE_TEUCHOS_COMPLEX
00056 #define STYPE std::complex<double>
00057 #else
00058 #define STYPE double
00059 #endif
00060 
00061 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
00062 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
00063 #ifdef HAVE_TEUCHOS_COMPLEX
00064 #define SCALARMAX  STYPE(10,0)
00065 #else
00066 #define SCALARMAX  STYPE(10)
00067 #endif
00068 
00069 template<typename TYPE>
00070 int PrintTestResults(std::string, TYPE, TYPE, bool);
00071 
00072 int ReturnCodeCheck(std::string, int, int, bool);
00073 
00074 typedef SerialDenseVector<OTYPE, STYPE> DVector;
00075 typedef SerialDenseMatrix<OTYPE, STYPE> DMatrix;
00076 
00077 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
00078 template<typename TYPE>
00079 TYPE GetRandom(TYPE, TYPE);
00080 
00081 // Returns a random integer between the two input parameters, inclusive
00082 template<>
00083 int GetRandom(int, int);
00084 
00085 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
00086 template<>
00087 double GetRandom(double, double);
00088 
00089 template<typename T>
00090 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
00091 
00092 // Generates random matrices and vectors using GetRandom()
00093 Teuchos::RCP<DMatrix> GetRandomMatrix(int m, int n);
00094 Teuchos::RCP<DVector> GetRandomVector(int n);
00095 
00096 // Compares the difference between two vectors using relative euclidean norms
00097 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
00098 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
00099                    const SerialDenseVector<OTYPE,STYPE>& Vector2,
00100                    ScalarTraits<STYPE>::magnitudeType Tolerance );
00101 
00102 int main(int argc, char* argv[])
00103 {
00104   typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
00105 
00106   int n=10, m=8;
00107   (void) m; // forestall "unused variable" compiler warning
00108   MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
00109 
00110   bool verbose = 0;
00111   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00112 
00113   if (verbose)
00114     std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
00115 
00116   int numberFailedTests = 0;
00117   int returnCode = 0;
00118   std::string testName = "", testType = "";
00119 
00120 #ifdef HAVE_TEUCHOS_COMPLEX
00121   testType = "COMPLEX";
00122 #else
00123   testType = "REAL";
00124 #endif
00125 
00126   if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
00127 
00128   // Create dense matrix and vector.
00129   Teuchos::RCP<DMatrix> A1 = GetRandomMatrix(n,n);
00130   Teuchos::RCP<DVector> x1 = GetRandomVector(n);
00131   DVector xhat(n), b(n), bt(n);
00132 
00133   // Compute the right-hand side vector using multiplication.
00134   returnCode = b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
00135   testName = "Generating right-hand side vector using A*x, where x is a random vector:";
00136   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00137 
00138   returnCode = bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
00139   testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
00140   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00141 
00142 #ifdef HAVE_TEUCHOS_COMPLEX
00143   DVector bct(n);
00144   returnCode = bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
00145   testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
00146   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00147 #endif
00148 
00149   // Fill the solution vector with zeros.
00150   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00151 
00152   // Create a serial dense solver.
00153   Teuchos::SerialDenseSolver<OTYPE, STYPE> solver1;
00154 
00155   // Pass in matrix and vectors
00156   solver1.setMatrix( A1 );
00157   solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00158 
00159   // Test1:  Simple factor and solve
00160   returnCode = solver1.factor();
00161   testName = "Simple solve: factor() random A:";
00162   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00163 
00164   // Non-transpose solve
00165   returnCode = solver1.solve();
00166   testName = "Simple solve: solve() random A (NO_TRANS):";
00167   numberFailedTests += CompareVectors( *x1, xhat, tol );
00168   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00169 
00170   // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00171   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00172   solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
00173   solver1.solveWithTransposeFlag( Teuchos::TRANS );
00174   returnCode = solver1.solve();
00175   testName = "Simple solve: solve() random A (TRANS):";
00176   numberFailedTests += CompareVectors( *x1, xhat, tol );
00177   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00178 
00179 #ifdef HAVE_TEUCHOS_COMPLEX
00180   // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00181   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00182   solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
00183   solver1.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
00184   returnCode = solver1.solve();
00185   testName = "Simple solve: solve() random A (CONJ_TRANS):";
00186   numberFailedTests += CompareVectors( *x1, xhat, tol );
00187   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00188 #endif
00189 
00190   // Test2: Invert the matrix, inverse should be in A.
00191   returnCode = solver1.invert();
00192   testName = "Simple solve: invert() random A:";
00193   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00194 
00195   // Compute the solution vector using multiplication and the inverse.
00196   returnCode = xhat.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, b, ScalarTraits<STYPE>::zero());
00197   testName = "Computing solution using inverted random A (NO_TRANS):";
00198   numberFailedTests += CompareVectors( *x1, xhat, tol );
00199   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00200 
00201   returnCode = xhat.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bt, ScalarTraits<STYPE>::zero());
00202   testName = "Computing solution using inverted random A (TRANS):";
00203   numberFailedTests += CompareVectors( *x1, xhat, tol );
00204   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00205 
00206 #ifdef HAVE_TEUCHOS_COMPLEX
00207   returnCode = xhat.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, bct, ScalarTraits<STYPE>::zero());
00208   testName = "Computing solution using inverted random A (CONJ_TRANS):";
00209   numberFailedTests += CompareVectors( *x1, xhat, tol );
00210   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00211 #endif
00212 
00213   // Test3:  Solve with iterative refinement.
00214 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00215   // Iterative refinement not implemented in Eigen
00216 #else
00217   // Create random linear system
00218   Teuchos::RCP<DMatrix> A2 = GetRandomMatrix(n,n);
00219   Teuchos::RCP<DVector> x2 = GetRandomVector(n);
00220 
00221   // Create LHS through multiplication with A2
00222   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00223   b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
00224   bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
00225 #ifdef HAVE_TEUCHOS_COMPLEX
00226   bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
00227 #endif
00228 
00229   // Create a serial dense solver.
00230   Teuchos::SerialDenseSolver<OTYPE, STYPE> solver2;
00231   solver2.solveToRefinedSolution( true );
00232 
00233   // Pass in matrix and vectors
00234   solver2.setMatrix( A2 );
00235   solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00236 
00237   // Factor and solve with iterative refinement.
00238   returnCode = solver2.factor();
00239   testName = "Solve with iterative refinement: factor() random A:";
00240   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00241 
00242   // Non-transpose solve
00243   returnCode = solver2.solve();
00244   testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
00245   numberFailedTests += CompareVectors( *x2, xhat, tol );
00246   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00247 
00248   // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00249   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00250   solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
00251   solver2.solveWithTransposeFlag( Teuchos::TRANS );
00252   returnCode = solver2.solve();
00253   testName = "Solve with iterative refinement: solve() random A (TRANS):";
00254   numberFailedTests += CompareVectors( *x2, xhat, tol );
00255   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00256 
00257 #ifdef HAVE_TEUCHOS_COMPLEX
00258   // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00259   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00260   solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
00261   solver2.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
00262   returnCode = solver2.solve();
00263   testName = "Solve with iterative refinement: solve() random A (CONJ_TRANS):";
00264   numberFailedTests += CompareVectors( *x2, xhat, tol );
00265   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00266 #endif
00267 #endif
00268 
00269   // Test4:  Solve with matrix equilibration.
00270 
00271   // Create random linear system
00272   Teuchos::RCP<DMatrix> A3 = GetRandomMatrix(n,n);
00273   Teuchos::RCP<DVector> x3 = GetRandomVector(n);
00274 
00275   // Create LHS through multiplication with A3
00276   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00277   b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
00278   bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
00279 #ifdef HAVE_TEUCHOS_COMPLEX
00280   bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
00281 #endif
00282 
00283   // Create a serial dense solver.
00284   Teuchos::SerialDenseSolver<OTYPE, STYPE> solver3;
00285   solver3.factorWithEquilibration( true );
00286 
00287   // Pass in matrix and vectors
00288   solver3.setMatrix( A3 );
00289   solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00290 
00291   // Factor and solve with matrix equilibration.
00292   returnCode = solver3.factor();
00293   testName = "Solve with matrix equilibration: factor() random A:";
00294   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00295 
00296   // Non-transpose solve
00297   returnCode = solver3.solve();
00298   testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
00299   numberFailedTests += CompareVectors( *x3, xhat, tol );
00300   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00301 
00302   // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00303   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00304   solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
00305   solver3.solveWithTransposeFlag( Teuchos::TRANS );
00306   returnCode = solver3.solve();
00307   testName = "Solve with matrix equilibration: solve() random A (TRANS):";
00308   numberFailedTests += CompareVectors( *x3, xhat, tol );
00309   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00310 
00311 #ifdef HAVE_TEUCHOS_COMPLEX
00312   // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
00313   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00314   solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
00315   solver3.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
00316   returnCode = solver3.solve();
00317   testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
00318   numberFailedTests += CompareVectors( *x3, xhat, tol );
00319   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00320 #endif
00321 
00322   //
00323   // If a test failed output the number of failed tests.
00324   //
00325   if(numberFailedTests > 0)
00326   {
00327             if (verbose) {
00328                 std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
00329                 std::cout << "End Result: TEST FAILED" << std::endl;
00330                 return -1;
00331             }
00332         }
00333   if(numberFailedTests == 0)
00334     std::cout << "End Result: TEST PASSED" << std::endl;
00335 
00336   return 0;
00337 }
00338 
00339 template<typename TYPE>
00340 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
00341 {
00342   int result;
00343   if(calculatedResult == expectedResult)
00344     {
00345       if(verbose) std::cout << testName << " successful." << std::endl;
00346       result = 0;
00347     }
00348   else
00349     {
00350       if(verbose) std::cout << testName << " unsuccessful." << std::endl;
00351       result = 1;
00352     }
00353   return result;
00354 }
00355 
00356 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
00357 {
00358   int result;
00359   if(expectedResult == 0)
00360     {
00361       if(returnCode == 0)
00362         {
00363           if(verbose) std::cout << testName << " test successful." << std::endl;
00364           result = 0;
00365         }
00366       else
00367         {
00368           if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
00369           result = 1;
00370         }
00371     }
00372   else
00373     {
00374       if(returnCode != 0)
00375         {
00376           if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
00377           result = 0;
00378         }
00379       else
00380         {
00381           if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
00382           result = 1;
00383         }
00384     }
00385   return result;
00386 }
00387 
00388 template<typename TYPE>
00389 TYPE GetRandom(TYPE Low, TYPE High)
00390 {
00391   return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
00392 }
00393 
00394 template<typename T>
00395 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
00396 {
00397   T lowMag = Low.real();
00398   T highMag = High.real();
00399   T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
00400   T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
00401   return std::complex<T>( real, imag );
00402 }
00403 
00404 template<>
00405 int GetRandom(int Low, int High)
00406 {
00407   return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
00408 }
00409 
00410 template<>
00411 double GetRandom(double Low, double High)
00412 {
00413   return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
00414 }
00415 
00416 Teuchos::RCP<DMatrix> GetRandomMatrix(int m, int n)
00417 {
00418   Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
00419 
00420   // Fill dense matrix with random entries.
00421   for (int i=0; i<m; i++)
00422     for (int j=0; j<n; j++)
00423       (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
00424 
00425   return newmat;
00426 }
00427 
00428 Teuchos::RCP<DVector> GetRandomVector(int n)
00429 {
00430   Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
00431 
00432   // Fill dense vector with random entries.
00433   for (int i=0; i<n; i++)
00434     (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
00435 
00436   return newvec;
00437 }
00438 
00439 /*  Function:  CompareVectors
00440     Purpose:   Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
00441 */
00442 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
00443                    const SerialDenseVector<OTYPE,STYPE>& Vector2,
00444                    ScalarTraits<STYPE>::magnitudeType Tolerance )
00445 {
00446   typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
00447 
00448   SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
00449   diff -= Vector2;
00450 
00451   MagnitudeType norm_diff = diff.normFrobenius();
00452   MagnitudeType norm_v2 = Vector2.normFrobenius();
00453   MagnitudeType temp = norm_diff;
00454   if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
00455     temp /= norm_v2;
00456 
00457   if (temp > Tolerance)
00458     return 1;
00459   else
00460     return 0;
00461 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines