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