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