Teuchos Package Browser (Single Doxygen Collection) Version of the Day
cxx_main_spd_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_SerialSymDenseMatrix.hpp"
00044 #include "Teuchos_SerialDenseVector.hpp"
00045 #include "Teuchos_SerialDenseHelpers.hpp"
00046 #include "Teuchos_SerialSpdDenseSolver.hpp"
00047 #include "Teuchos_SerialQRDenseSolver.hpp"
00048 #include "Teuchos_ScalarTraits.hpp"
00049 #include "Teuchos_RCP.hpp"
00050 #include "Teuchos_Version.hpp"
00051 
00052 using Teuchos::ScalarTraits;
00053 using Teuchos::SerialDenseMatrix;
00054 using Teuchos::SerialSymDenseMatrix;
00055 using Teuchos::SerialDenseVector;
00056 
00057 #define OTYPE int
00058 #ifdef HAVE_TEUCHOS_COMPLEX
00059 #define STYPE std::complex<double>
00060 #else
00061 #define STYPE double
00062 #endif
00063 
00064 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
00065 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
00066 #ifdef HAVE_TEUCHOS_COMPLEX
00067 #define SCALARMAX  STYPE(10,0)
00068 #else
00069 #define SCALARMAX  STYPE(10)
00070 #endif
00071 
00072 template<typename TYPE>
00073 int PrintTestResults(std::string, TYPE, TYPE, bool);
00074 
00075 int ReturnCodeCheck(std::string, int, int, bool);
00076 
00077 typedef SerialDenseVector<OTYPE, STYPE> DVector;
00078 typedef SerialDenseMatrix<OTYPE, STYPE> DMatrix;
00079 typedef SerialSymDenseMatrix<OTYPE, STYPE> SDMatrix;
00080 
00081 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
00082 template<typename TYPE>
00083 TYPE GetRandom(TYPE, TYPE);
00084 
00085 // Returns a random integer between the two input parameters, inclusive
00086 template<>
00087 int GetRandom(int, int);
00088 
00089 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
00090 template<>
00091 double GetRandom(double, double);
00092 
00093 template<typename T>
00094 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
00095 
00096 // Generates random matrices and vectors using GetRandom()
00097 Teuchos::RCP<SDMatrix> GetRandomSpdMatrix(int n);
00098 Teuchos::RCP<DMatrix> GetRandomMatrix(int m, int n);
00099 Teuchos::RCP<DVector> GetRandomVector(int n);
00100 
00101 // Compares the difference between two vectors using relative euclidean norms
00102 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
00103 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
00104        const SerialDenseVector<OTYPE,STYPE>& Vector2,
00105        ScalarTraits<STYPE>::magnitudeType Tolerance );
00106 
00107 int main(int argc, char* argv[])
00108 {
00109   typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
00110 
00111   int n=10;
00112   MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
00113 
00114   bool verbose = 0;
00115   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00116 
00117   if (verbose)
00118     std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
00119 
00120   int numberFailedTests = 0, failedComparison = 0;
00121   int returnCode = 0;
00122   std::string testName = "", testType = "";
00123 
00124 #ifdef HAVE_TEUCHOS_COMPLEX
00125   testType = "COMPLEX";
00126 #else
00127   testType = "REAL";
00128 #endif
00129 
00130   if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL SPD DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
00131 
00132   // Create dense matrix and vector.
00133   Teuchos::RCP<SDMatrix> A1 = GetRandomSpdMatrix(n);
00134   Teuchos::RCP<DVector> x1 = GetRandomVector(n);
00135   DVector xhat(n), b(n), bt(n);
00136 
00137   Teuchos::RCP<DMatrix> A1full = Teuchos::rcp( new DMatrix(n,n) );
00138   for (int j=0; j<n; ++j) {
00139     for (int i=0; i<=j; ++i) {
00140       (*A1full)(i,j) = (*A1)(i,j);
00141     }
00142   }
00143   for (int j=0; j<n; ++j) {
00144     for (int i=j+1; i<n; ++i) {
00145       (*A1full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A1)(i,j) );
00146     }
00147   }
00148 
00149   // Compute the right-hand side vector using multiplication.
00150   returnCode = b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1full, *x1, ScalarTraits<STYPE>::zero());
00151   testName = "Generating right-hand side vector using A*x, where x is a random vector:";
00152   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00153 
00154   // Fill the solution vector with zeros.
00155   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00156 
00157   // Create a serial spd dense solver.
00158   Teuchos::SerialSpdDenseSolver<OTYPE, STYPE> solver1;
00159 
00160   // Pass in matrix and vectors
00161   solver1.setMatrix( A1 );
00162   solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00163 
00164   // Test1:  Simple factor and solve
00165   returnCode = solver1.factor();
00166   testName = "Simple solve: factor() random A:";
00167   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00168 
00169   // Non-transpose solve
00170   returnCode = solver1.solve();
00171   testName = "Simple solve: solve() random A (NO_TRANS):";
00172   failedComparison = CompareVectors( *x1, xhat, tol );
00173   if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
00174   numberFailedTests += failedComparison;
00175   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00176 
00177   // Test2:  Solve with iterative refinement.
00178 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
00179   // Iterative refinement not implemented in Eigen
00180 #else
00181   // Create random linear system
00182   Teuchos::RCP<SDMatrix> A2 = GetRandomSpdMatrix(n);
00183   Teuchos::RCP<DVector> x2 = GetRandomVector(n);
00184 
00185   Teuchos::RCP<DMatrix> A2full = Teuchos::rcp( new DMatrix(n,n) );
00186   for (int j=0; j<n; ++j) {
00187     for (int i=0; i<=j; ++i) {
00188       (*A2full)(i,j) = (*A2)(i,j);
00189     }
00190   }
00191   for (int j=0; j<n; ++j) {
00192     for (int i=j+1; i<n; ++i) {
00193       (*A2full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A2)(i,j) );
00194     }
00195   }
00196 
00197   // Create LHS through multiplication with A2
00198   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00199   b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2full, *x2, ScalarTraits<STYPE>::zero());
00200 
00201   // Create a serial dense solver.
00202   Teuchos::SerialSpdDenseSolver<OTYPE, STYPE> solver2;
00203   solver2.solveToRefinedSolution( true );
00204 
00205   // Pass in matrix and vectors
00206   solver2.setMatrix( A2 );
00207   solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00208 
00209   // Factor and solve with iterative refinement.
00210   returnCode = solver2.factor();
00211   testName = "Solve with iterative refinement: factor() random A:";
00212   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00213 
00214   // Non-transpose solve
00215   returnCode = solver2.solve();
00216   testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
00217   failedComparison = CompareVectors( *x2, xhat, tol );
00218   if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
00219   numberFailedTests += failedComparison;
00220   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00221 
00222 #endif
00223 
00224   // Test3:  Solve with matrix equilibration.
00225 
00226   // Create random linear system
00227   Teuchos::RCP<SDMatrix> A3 = GetRandomSpdMatrix(n);
00228   Teuchos::RCP<DVector> x3 = GetRandomVector(n);
00229 
00230   Teuchos::RCP<DMatrix> A3full = Teuchos::rcp( new DMatrix(n,n) );
00231   for (int j=0; j<n; ++j) {
00232     for (int i=0; i<=j; ++i) {
00233       (*A3full)(i,j) = (*A3)(i,j);
00234     }
00235   }
00236   for (int j=0; j<n; ++j) {
00237     for (int i=j+1; i<n; ++i) {
00238       (*A3full)(i,j) = ScalarTraits<STYPE>::conjugate( (*A3)(i,j) );
00239     }
00240   }
00241 
00242   // Create LHS through multiplication with A3
00243   xhat.putScalar( ScalarTraits<STYPE>::zero() );
00244   b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3full, *x3, ScalarTraits<STYPE>::zero());
00245 
00246   // Create a serial dense solver.
00247   Teuchos::SerialSpdDenseSolver<OTYPE, STYPE> solver3;
00248   solver3.factorWithEquilibration( true );
00249 
00250   // Pass in matrix and vectors
00251   solver3.setMatrix( A3 );
00252   solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
00253 
00254   // Factor and solve with matrix equilibration.
00255   returnCode = solver3.factor();
00256   testName = "Solve with matrix equilibration: factor() random A:";
00257   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00258 
00259   // Non-transpose solve
00260   returnCode = solver3.solve();
00261   testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
00262   failedComparison = CompareVectors( *x3, xhat, tol );
00263   if(verbose && failedComparison>0) std::cout << "COMPARISON FAILED : ";
00264   numberFailedTests += failedComparison;
00265   numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
00266 
00267   //
00268   // If a test failed output the number of failed tests.
00269   //
00270   if(numberFailedTests > 0)
00271   {
00272       if (verbose) {
00273     std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
00274     std::cout << "End Result: TEST FAILED" << std::endl;
00275     return -1;
00276       }
00277   }
00278   if(numberFailedTests == 0)
00279     std::cout << "End Result: TEST PASSED" << std::endl;
00280 
00281   return 0;
00282 
00283 }
00284 
00285 template<typename TYPE>
00286 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
00287 {
00288   int result;
00289   if(calculatedResult == expectedResult)
00290     {
00291       if(verbose) std::cout << testName << " successful." << std::endl;
00292       result = 0;
00293     }
00294   else
00295     {
00296       if(verbose) std::cout << testName << " unsuccessful." << std::endl;
00297       result = 1;
00298     }
00299   return result;
00300 }
00301 
00302 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
00303 {
00304   int result;
00305   if(expectedResult == 0)
00306     {
00307       if(returnCode == 0)
00308   {
00309     if(verbose) std::cout << testName << " test successful." << std::endl;
00310     result = 0;
00311   }
00312       else
00313   {
00314     if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
00315     result = 1;
00316   }
00317     }
00318   else
00319     {
00320       if(returnCode != 0)
00321   {
00322     if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
00323     result = 0;
00324   }
00325       else
00326   {
00327     if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
00328     result = 1;
00329   }
00330     }
00331   return result;
00332 }
00333 
00334 template<typename TYPE>
00335 TYPE GetRandom(TYPE Low, TYPE High)
00336 {
00337   return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
00338 }
00339 
00340 template<typename T>
00341 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
00342 {
00343   T lowMag = Low.real();
00344   T highMag = High.real();
00345   T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
00346   T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
00347   return std::complex<T>( real, imag );
00348 }
00349 
00350 template<>
00351 int GetRandom(int Low, int High)
00352 {
00353   return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
00354 }
00355 
00356 template<>
00357 double GetRandom(double Low, double High)
00358 {
00359   return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
00360 }
00361 
00362 Teuchos::RCP<SDMatrix> GetRandomSpdMatrix(int n)
00363 {
00364 
00365   Teuchos::RCP<SDMatrix> mat = Teuchos::rcp( new SDMatrix(n) );
00366   Teuchos::RCP<DMatrix> tmp = Teuchos::rcp( new DMatrix(n,n) );
00367   Teuchos::RCP<DMatrix> tmp2 = Teuchos::rcp( new DMatrix(n,n) );
00368   Teuchos::RCP<DMatrix> Q = Teuchos::rcp( new DMatrix(n,n) );
00369   Teuchos::RCP<DMatrix> D = Teuchos::rcp( new DMatrix(n,n) );
00370 
00371   // QR factor a random matrix and get the orthogonal Q
00372   Teuchos::RCP<DMatrix> A = GetRandomMatrix(n,n);
00373   Teuchos::SerialQRDenseSolver<OTYPE, STYPE> solver;
00374   solver.setMatrix( A );
00375   solver.factor();
00376   solver.formQ();
00377   Q = solver.getQ();
00378 
00379   // Get n random positive eigenvalues and put them in a diagonal matrix D
00380   Teuchos::RCP<DVector> E = GetRandomVector(n);
00381   for (int i=0; i<n; ++i) {
00382     (*D)(i,i) = ScalarTraits<STYPE>::magnitude(ScalarTraits<STYPE>::one()) +
00383                 ScalarTraits<STYPE>::magnitude( (*E)(i) ) * ScalarTraits<STYPE>::magnitude( (*E)(i) );
00384   }
00385 
00386   // Form the spd matrix Q' D Q
00387   tmp->multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *D, *Q, ScalarTraits<STYPE>::zero());
00388   tmp2->multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *Q, *tmp, ScalarTraits<STYPE>::zero());
00389 
00390   mat->setUpper();
00391   for (int j=0; j<n; ++j) {
00392     for (int i=0; i<=j; ++i) {
00393       (*mat)(i,j) = (*tmp2)(i,j);
00394     }
00395   }
00396 
00397   return mat;
00398 }
00399 
00400 Teuchos::RCP<DMatrix> GetRandomMatrix(int m, int n)
00401 {
00402   Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
00403 
00404   // Fill dense matrix with random entries.
00405   for (int i=0; i<m; i++)
00406     for (int j=0; j<n; j++)
00407       (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
00408 
00409   return newmat;
00410 }
00411 
00412 Teuchos::RCP<DVector> GetRandomVector(int n)
00413 {
00414   Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
00415 
00416   // Fill dense vector with random entries.
00417   for (int i=0; i<n; i++)
00418     (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
00419 
00420   return newvec;
00421 }
00422 
00423 /*  Function:  CompareVectors
00424     Purpose:   Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
00425 */
00426 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
00427        const SerialDenseVector<OTYPE,STYPE>& Vector2,
00428        ScalarTraits<STYPE>::magnitudeType Tolerance )
00429 {
00430   typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
00431 
00432   SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
00433   diff -= Vector2;
00434 
00435   MagnitudeType norm_diff = diff.normFrobenius();
00436   MagnitudeType norm_v2 = Vector2.normFrobenius();
00437   MagnitudeType temp = norm_diff;
00438   if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
00439     temp /= norm_v2;
00440 
00441   if (temp > Tolerance)
00442     return 1;
00443   else
00444     return 0;
00445 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines