diff --git a/packages/belos/epetra/test/MINRES/CMakeLists.txt b/packages/belos/epetra/test/MINRES/CMakeLists.txt index 3bf807fbb3ed..09a6a9a46fae 100644 --- a/packages/belos/epetra/test/MINRES/CMakeLists.txt +++ b/packages/belos/epetra/test/MINRES/CMakeLists.txt @@ -1,4 +1,10 @@ +TRIBITS_ADD_EXECUTABLE_AND_TEST( + minres_diag + SOURCES test_minres_diag.cpp + COMM serial mpi + STANDARD_PASS_OUTPUT +) TRIBITS_ADD_EXECUTABLE_AND_TEST( minres_indefinite @@ -18,6 +24,7 @@ IF (${PACKAGE_NAME}_ENABLE_Triutils) COMM serial mpi ARGS "--verbose --filename=bcsstk14.hb --tol=1e-5" + "--verbose --filename=bcsstk14.hb --num-rhs=2 --tol=1e-5" STANDARD_PASS_OUTPUT ) diff --git a/packages/belos/epetra/test/MINRES/test_minres_diag.cpp b/packages/belos/epetra/test/MINRES/test_minres_diag.cpp new file mode 100644 index 000000000000..06acaef3b2ac --- /dev/null +++ b/packages/belos/epetra/test/MINRES/test_minres_diag.cpp @@ -0,0 +1,422 @@ +//@HEADER +// ************************************************************************ +// +// Belos: Block Linear Solvers Package +// Copyright 2004 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// ************************************************************************ +//@HEADER +// +// This test generates diagonal matrices for MINRES to solve. +// +// NOTE: No preconditioner is used in this case. +// + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include "Teuchos_StandardCatchMacros.hpp" + +#ifdef EPETRA_MPI +#include +#include +#else +#include +#endif + +using std::vector; +using namespace Belos; + +//************************************************************************************************ + +class Vector_Operator +{ + public: + + Vector_Operator(int m_in, int n_in) : m(m_in), n(n_in) {}; + + virtual ~Vector_Operator() {}; + + virtual void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) = 0; + + int size (int dim) const { return (dim == 1) ? m : n; }; + + protected: + + int m, n; // an (m x n) operator + + private: + + // Not allowing copy construction. + Vector_Operator( const Vector_Operator& ): m(0), n(0) {}; + Vector_Operator* operator=( const Vector_Operator& ) { return NULL; }; + +}; + +//************************************************************************************************ + +class Diagonal_Operator : public Vector_Operator +{ + public: + + Diagonal_Operator(int n_in, double v_in) : Vector_Operator(n_in, n_in), v(v_in) { }; + + ~Diagonal_Operator() { }; + + void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) + { + y.Scale( v, x ); + }; + + private: + + double v; +}; + +//************************************************************************************************ + +class Diagonal_Operator_2 : public Vector_Operator +{ + public: + + Diagonal_Operator_2(int n_in, int min_gid_in, double v_in) + : Vector_Operator(n_in, n_in), min_gid(min_gid_in), v(v_in) {} + + ~Diagonal_Operator_2() { }; + + void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) + { + int myCols = y.MyLength(); + for (int j=0; j < x.NumVectors(); ++j) { + for (int i=0; i < myCols; ++i) (*y(j))[i] = (min_gid+i+1)*v*(*x(j))[i]; // NOTE: square operator! + } + }; + + private: + + int min_gid; + double v; +}; + +//************************************************************************************************ + +class Composed_Operator : public Vector_Operator +{ + public: + + Composed_Operator(int n, + const Teuchos::RCP& pA_in, + const Teuchos::RCP& pB_in); + + virtual ~Composed_Operator() {}; + + virtual void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y); + + private: + + Teuchos::RCP pA; + Teuchos::RCP pB; +}; + +Composed_Operator::Composed_Operator(int n_in, + const Teuchos::RCP& pA_in, + const Teuchos::RCP& pB_in) +: Vector_Operator(n_in, n_in), pA(pA_in), pB(pB_in) +{ +} + +void Composed_Operator::operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) +{ + Epetra_MultiVector ytemp(y.Map(), y.NumVectors(), false); + (*pB)( x, ytemp ); + (*pA)( ytemp, y ); +} + +//************************************************************************************************ + +class Trilinos_Interface : public Epetra_Operator +{ + public: + + Trilinos_Interface(const Teuchos::RCP pA_in, + const Teuchos::RCP pComm_in, + const Teuchos::RCP pMap_in) + : pA (pA_in), + pComm (pComm_in), + pMap (pMap_in), + use_transpose (false) + {} + + int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const; + + int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const + { + return(Apply(X,Y)); // No inverse + }; + + virtual ~Trilinos_Interface() {}; + + const char * Label() const {return("Trilinos_Interface, an Epetra_Operator implementation");}; + + bool UseTranspose() const {return(use_transpose);}; // always set to false + + int SetUseTranspose(bool UseTranspose_in) { use_transpose = false; return(-1); }; + + bool HasNormInf() const {return(false);}; // cannot return inf-norm + + double NormInf() const {return(0.0);}; + + virtual const Epetra_Comm & Comm() const {return *pComm; } + + virtual const Epetra_Map & OperatorDomainMap() const {return *pMap; } + + virtual const Epetra_Map & OperatorRangeMap() const {return *pMap; } + + private: + + Teuchos::RCP pA; + Teuchos::RCP pComm; + Teuchos::RCP pMap; + + bool use_transpose; +}; + +int Trilinos_Interface::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const +{ + (*pA)(X,Y); + + return(0); +} + +//************************************************************************************************ + +class Iterative_Inverse_Operator : public Vector_Operator +{ + + typedef Epetra_MultiVector MV; + typedef Epetra_Operator OP; + + public: + + Iterative_Inverse_Operator(int n_in, int blocksize, + const Teuchos::RCP& pA_in, + std::string opString="Iterative Solver", bool print_in=false); + + virtual ~Iterative_Inverse_Operator() {} + + virtual void operator () (const Epetra_MultiVector &b, Epetra_MultiVector &x); + + private: + + Teuchos::RCP pA; // operator which will be inverted + // supplies a matrix std::vector multiply + const bool print; + + Teuchos::Time timer; + Teuchos::RCP pComm; + Teuchos::RCP pMap; + + Teuchos::RCP pPE; + Teuchos::RCP pList; + Teuchos::RCP > pProb; + Teuchos::RCP > pBelos; +}; + +Iterative_Inverse_Operator::Iterative_Inverse_Operator(int n_in, int blocksize, + const Teuchos::RCP& pA_in, + std::string opString, bool print_in) +: Vector_Operator(n_in, n_in), // square operator + pA(pA_in), + print(print_in), + timer(opString) +{ + + int n_global; + +#ifdef EPETRA_MPI + MPI_Allreduce(&n, &n_global, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + pComm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) ); +#else + pComm = Teuchos::rcp( new Epetra_SerialComm() ); + n_global = n; +#endif + pMap = Teuchos::rcp( new Epetra_Map(n_global, n, 0, *pComm) ); + + pPE = Teuchos::rcp( new Trilinos_Interface(pA, pComm, pMap ) ); + + pProb = Teuchos::rcp( new LinearProblem() ); + pProb->setOperator( pPE ); + + int max_iter = 100; + double tol = 1.0e-10; + int verbosity = Belos::Errors + Belos::Warnings; + if (print) + verbosity += Belos::TimingDetails + Belos::StatusTestDetails; + + pList = Teuchos::rcp( new Teuchos::ParameterList ); + pList->set( "Maximum Iterations", max_iter ); + pList->set( "Convergence Tolerance", tol ); + pList->set( "Verbosity", verbosity ); + + pBelos = Teuchos::rcp( new MinresSolMgr(pProb, pList) ); +} + +void Iterative_Inverse_Operator::operator () (const Epetra_MultiVector &b, Epetra_MultiVector &x) +{ + int pid = pComm->MyPID(); + + // Initialize the solution to zero + x.PutScalar( 0.0 ); + + // Reset the solver, problem, and status test for next solve (HKT) + pProb->setProblem( Teuchos::rcp(&x, false), Teuchos::rcp(&b, false) ); + + timer.start(); + Belos::ReturnType ret = pBelos->solve(); + timer.stop(); + + if (pid == 0 && print) { + if (ret == Belos::Converged) + { + std::cout << std::endl << "pid[" << pid << "] Minres converged" << std::endl; + std::cout << "Solution time: " << timer.totalElapsedTime() << std::endl; + + } + else + std::cout << std::endl << "pid[" << pid << "] Minres did not converge" << std::endl; + } +} + +//************************************************************************************************ +//************************************************************************************************ + +int main(int argc, char *argv[]) +{ + int pid = -1; + + Teuchos::GlobalMPISession session(&argc, &argv, NULL); + +#ifdef EPETRA_MPI + Epetra_MpiComm Comm(MPI_COMM_WORLD); +#else + Epetra_SerialComm Comm; +#endif + + bool verbose = false; + bool success = true; + try { + + pid = Comm.MyPID(); + + int n(10); + int numRHS=1; + + Epetra_Map Map = Epetra_Map(n, 0, Comm); + Epetra_MultiVector X(Map, numRHS, false), Y(Map, numRHS, false); + X.PutScalar( 1.0 ); + + // Inner computes inv(D2)*y + Teuchos::RCP D2 = Teuchos::rcp(new Diagonal_Operator_2(n, Map.MinMyGID(), 1.0)); + Iterative_Inverse_Operator A2(n, 1, D2, "Belos (inv(D2))", true); + + // should return x=(1, 1/2, 1/3, ..., 1/10) + A2(X,Y); + + if (pid==0) { + std::cout << "Vector Y should have all entries [1, 1/2, 1/3, ..., 1/10]" << std::endl; + } + Y.Print(std::cout); + + // Inner computes inv(D)*x + Teuchos::RCP D = Teuchos::rcp(new Diagonal_Operator(n, 4.0)); + Teuchos::RCP Inner = + Teuchos::rcp(new Iterative_Inverse_Operator(n, 1, D, "Belos (inv(D))", false)); + + // Composed_Operator computed inv(D)*B*x + Teuchos::RCP B = Teuchos::rcp(new Diagonal_Operator(n, 4.0)); + Teuchos::RCP C = Teuchos::rcp(new Composed_Operator(n, Inner, B)); + + // Outer computes inv(C) = inv(inv(D)*B)*x = inv(B)*D*x = x + Teuchos::RCP Outer = + Teuchos::rcp(new Iterative_Inverse_Operator(n, 1, C, "Belos (inv(C)=inv(inv(D)*B))", true)); + + // should return x=1/4 + (*Inner)(X,Y); + + if (pid==0) { + std::cout << std::endl << "Vector Y should have all entries [1/4, 1/4, 1/4, ..., 1/4]" << std::endl; + } + Y.Print(std::cout); + + // should return x=1 + (*Outer)(X,Y); + + if (pid==0) { + std::cout << "Vector Y should have all entries [1, 1, 1, ..., 1]" << std::endl; + } + Y.Print(std::cout); + + // Compute the norm of Y - 1.0 + std::vector norm_Y(Y.NumVectors()); + Y.Update(-1.0, X, 1.0); + Y.Norm2(&norm_Y[0]); + + if (pid==0) + std::cout << "Two-norm of std::vector (Y-1.0) : "<< norm_Y[0] << std::endl; + + success = (norm_Y[0] < 1e-10 && !Teuchos::ScalarTraits::isnaninf( norm_Y[0] ) ); + + if (success) { + if (pid==0) + std::cout << "End Result: TEST PASSED" << std::endl; + } else { + if (pid==0) + std::cout << "End Result: TEST FAILED" << std::endl; + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/packages/belos/epetra/test/MINRES/test_minres_hb.cpp b/packages/belos/epetra/test/MINRES/test_minres_hb.cpp index 265c502596fd..e864bcf04437 100644 --- a/packages/belos/epetra/test/MINRES/test_minres_hb.cpp +++ b/packages/belos/epetra/test/MINRES/test_minres_hb.cpp @@ -119,10 +119,11 @@ main (int argc, char *argv[]) RCP A; RCP B, X; RCP rowMap; + int tmp_numRHS = numRHS; try { // This might change the number of right-hand sides, if we read in // a right-hand side from the Harwell-Boeing file. - Belos::Util::createEpetraProblem (filename, &rowMap, &A, &B, &X, &MyPID, numRHS); + Belos::Util::createEpetraProblem (filename, &rowMap, &A, &B, &X, &MyPID, tmp_numRHS); } catch (std::exception& e) { TEUCHOS_TEST_FOR_EXCEPTION (true, std::runtime_error, "Failed to create Epetra problem for matrix " @@ -130,6 +131,18 @@ main (int argc, char *argv[]) "createEpetraProblem() reports the following " "error: " << e.what()); } + + // + // *****Construct initial guess and random right-hand-sides ***** + // + if (tmp_numRHS != numRHS) + { + X = rcp( new Epetra_MultiVector( A->Map(), numRHS ) ); + MVT::MvRandom( *X ); + B = rcp( new Epetra_MultiVector( A->Map(), numRHS ) ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + } // // Compute the initial residual norm of the problem, so we can see // by how much it improved after the solve. @@ -310,11 +323,6 @@ main (int argc, char *argv[]) } } } - -# ifdef BELOS_TEUCHOS_TIME_MONITOR - Teuchos::TimeMonitor::summarize (verbOut); -# endif // BELOS_TEUCHOS_TIME_MONITOR - success = (ret == Belos::Converged && !badRes); if (success) { diff --git a/packages/belos/epetra/test/TFQMR/CMakeLists.txt b/packages/belos/epetra/test/TFQMR/CMakeLists.txt index f268e91a756d..a64972845d0b 100644 --- a/packages/belos/epetra/test/TFQMR/CMakeLists.txt +++ b/packages/belos/epetra/test/TFQMR/CMakeLists.txt @@ -1,4 +1,10 @@ +TRIBITS_ADD_EXECUTABLE_AND_TEST( + tfqmr_diag + SOURCES test_tfqmr_diag.cpp + COMM serial mpi + STANDARD_PASS_OUTPUT +) ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) IF (${PACKAGE_NAME}_ENABLE_Triutils) diff --git a/packages/belos/epetra/test/TFQMR/test_tfqmr_diag.cpp b/packages/belos/epetra/test/TFQMR/test_tfqmr_diag.cpp new file mode 100644 index 000000000000..60f3bede6fed --- /dev/null +++ b/packages/belos/epetra/test/TFQMR/test_tfqmr_diag.cpp @@ -0,0 +1,422 @@ +//@HEADER +// ************************************************************************ +// +// Belos: Block Linear Solvers Package +// Copyright 2004 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// ************************************************************************ +//@HEADER +// +// This test generates diagonal matrices for TFQMR to solve. +// +// NOTE: No preconditioner is used in this case. +// + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include "Teuchos_StandardCatchMacros.hpp" + +#ifdef EPETRA_MPI +#include +#include +#else +#include +#endif + +using std::vector; +using namespace Belos; + +//************************************************************************************************ + +class Vector_Operator +{ + public: + + Vector_Operator(int m_in, int n_in) : m(m_in), n(n_in) {}; + + virtual ~Vector_Operator() {}; + + virtual void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) = 0; + + int size (int dim) const { return (dim == 1) ? m : n; }; + + protected: + + int m, n; // an (m x n) operator + + private: + + // Not allowing copy construction. + Vector_Operator( const Vector_Operator& ): m(0), n(0) {}; + Vector_Operator* operator=( const Vector_Operator& ) { return NULL; }; + +}; + +//************************************************************************************************ + +class Diagonal_Operator : public Vector_Operator +{ + public: + + Diagonal_Operator(int n_in, double v_in) : Vector_Operator(n_in, n_in), v(v_in) { }; + + ~Diagonal_Operator() { }; + + void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) + { + y.Scale( v, x ); + }; + + private: + + double v; +}; + +//************************************************************************************************ + +class Diagonal_Operator_2 : public Vector_Operator +{ + public: + + Diagonal_Operator_2(int n_in, int min_gid_in, double v_in) + : Vector_Operator(n_in, n_in), min_gid(min_gid_in), v(v_in) {} + + ~Diagonal_Operator_2() { }; + + void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) + { + int myCols = y.MyLength(); + for (int j=0; j < x.NumVectors(); ++j) { + for (int i=0; i < myCols; ++i) (*y(j))[i] = (min_gid+i+1)*v*(*x(j))[i]; // NOTE: square operator! + } + }; + + private: + + int min_gid; + double v; +}; + +//************************************************************************************************ + +class Composed_Operator : public Vector_Operator +{ + public: + + Composed_Operator(int n, + const Teuchos::RCP& pA_in, + const Teuchos::RCP& pB_in); + + virtual ~Composed_Operator() {}; + + virtual void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y); + + private: + + Teuchos::RCP pA; + Teuchos::RCP pB; +}; + +Composed_Operator::Composed_Operator(int n_in, + const Teuchos::RCP& pA_in, + const Teuchos::RCP& pB_in) +: Vector_Operator(n_in, n_in), pA(pA_in), pB(pB_in) +{ +} + +void Composed_Operator::operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y) +{ + Epetra_MultiVector ytemp(y.Map(), y.NumVectors(), false); + (*pB)( x, ytemp ); + (*pA)( ytemp, y ); +} + +//************************************************************************************************ + +class Trilinos_Interface : public Epetra_Operator +{ + public: + + Trilinos_Interface(const Teuchos::RCP pA_in, + const Teuchos::RCP pComm_in, + const Teuchos::RCP pMap_in) + : pA (pA_in), + pComm (pComm_in), + pMap (pMap_in), + use_transpose (false) + {} + + int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const; + + int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const + { + return(Apply(X,Y)); // No inverse + }; + + virtual ~Trilinos_Interface() {}; + + const char * Label() const {return("Trilinos_Interface, an Epetra_Operator implementation");}; + + bool UseTranspose() const {return(use_transpose);}; // always set to false + + int SetUseTranspose(bool UseTranspose_in) { use_transpose = false; return(-1); }; + + bool HasNormInf() const {return(false);}; // cannot return inf-norm + + double NormInf() const {return(0.0);}; + + virtual const Epetra_Comm & Comm() const {return *pComm; } + + virtual const Epetra_Map & OperatorDomainMap() const {return *pMap; } + + virtual const Epetra_Map & OperatorRangeMap() const {return *pMap; } + + private: + + Teuchos::RCP pA; + Teuchos::RCP pComm; + Teuchos::RCP pMap; + + bool use_transpose; +}; + +int Trilinos_Interface::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const +{ + (*pA)(X,Y); + + return(0); +} + +//************************************************************************************************ + +class Iterative_Inverse_Operator : public Vector_Operator +{ + + typedef Epetra_MultiVector MV; + typedef Epetra_Operator OP; + + public: + + Iterative_Inverse_Operator(int n_in, int blocksize, + const Teuchos::RCP& pA_in, + std::string opString="Iterative Solver", bool print_in=false); + + virtual ~Iterative_Inverse_Operator() {} + + virtual void operator () (const Epetra_MultiVector &b, Epetra_MultiVector &x); + + private: + + Teuchos::RCP pA; // operator which will be inverted + // supplies a matrix std::vector multiply + const bool print; + + Teuchos::Time timer; + Teuchos::RCP pComm; + Teuchos::RCP pMap; + + Teuchos::RCP pPE; + Teuchos::RCP pList; + Teuchos::RCP > pProb; + Teuchos::RCP > pBelos; +}; + +Iterative_Inverse_Operator::Iterative_Inverse_Operator(int n_in, int blocksize, + const Teuchos::RCP& pA_in, + std::string opString, bool print_in) +: Vector_Operator(n_in, n_in), // square operator + pA(pA_in), + print(print_in), + timer(opString) +{ + + int n_global; + +#ifdef EPETRA_MPI + MPI_Allreduce(&n, &n_global, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + pComm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) ); +#else + pComm = Teuchos::rcp( new Epetra_SerialComm() ); + n_global = n; +#endif + pMap = Teuchos::rcp( new Epetra_Map(n_global, n, 0, *pComm) ); + + pPE = Teuchos::rcp( new Trilinos_Interface(pA, pComm, pMap ) ); + + pProb = Teuchos::rcp( new LinearProblem() ); + pProb->setOperator( pPE ); + + int max_iter = 100; + double tol = 1.0e-10; + int verbosity = Belos::Errors + Belos::Warnings; + if (print) + verbosity += Belos::TimingDetails + Belos::StatusTestDetails; + + pList = Teuchos::rcp( new Teuchos::ParameterList ); + pList->set( "Maximum Iterations", max_iter ); + pList->set( "Convergence Tolerance", tol ); + pList->set( "Verbosity", verbosity ); + + pBelos = Teuchos::rcp( new TFQMRSolMgr(pProb, pList) ); +} + +void Iterative_Inverse_Operator::operator () (const Epetra_MultiVector &b, Epetra_MultiVector &x) +{ + int pid = pComm->MyPID(); + + // Initialize the solution to zero + x.PutScalar( 0.0 ); + + // Reset the solver, problem, and status test for next solve (HKT) + pProb->setProblem( Teuchos::rcp(&x, false), Teuchos::rcp(&b, false) ); + + timer.start(); + Belos::ReturnType ret = pBelos->solve(); + timer.stop(); + + if (pid == 0 && print) { + if (ret == Belos::Converged) + { + std::cout << std::endl << "pid[" << pid << "] TFQMR converged" << std::endl; + std::cout << "Solution time: " << timer.totalElapsedTime() << std::endl; + + } + else + std::cout << std::endl << "pid[" << pid << "] TFQMR did not converge" << std::endl; + } +} + +//************************************************************************************************ +//************************************************************************************************ + +int main(int argc, char *argv[]) +{ + int pid = -1; + + Teuchos::GlobalMPISession session(&argc, &argv, NULL); + +#ifdef EPETRA_MPI + Epetra_MpiComm Comm(MPI_COMM_WORLD); +#else + Epetra_SerialComm Comm; +#endif + + bool verbose = false; + bool success = true; + try { + + pid = Comm.MyPID(); + + int n(10); + int numRHS=1; + + Epetra_Map Map = Epetra_Map(n, 0, Comm); + Epetra_MultiVector X(Map, numRHS, false), Y(Map, numRHS, false); + X.PutScalar( 1.0 ); + + // Inner computes inv(D2)*y + Teuchos::RCP D2 = Teuchos::rcp(new Diagonal_Operator_2(n, Map.MinMyGID(), 1.0)); + Iterative_Inverse_Operator A2(n, 1, D2, "Belos (inv(D2))", true); + + // should return x=(1, 1/2, 1/3, ..., 1/10) + A2(X,Y); + + if (pid==0) { + std::cout << "Vector Y should have all entries [1, 1/2, 1/3, ..., 1/10]" << std::endl; + } + Y.Print(std::cout); + + // Inner computes inv(D)*x + Teuchos::RCP D = Teuchos::rcp(new Diagonal_Operator(n, 4.0)); + Teuchos::RCP Inner = + Teuchos::rcp(new Iterative_Inverse_Operator(n, 1, D, "Belos (inv(D))", false)); + + // Composed_Operator computed inv(D)*B*x + Teuchos::RCP B = Teuchos::rcp(new Diagonal_Operator(n, 4.0)); + Teuchos::RCP C = Teuchos::rcp(new Composed_Operator(n, Inner, B)); + + // Outer computes inv(C) = inv(inv(D)*B)*x = inv(B)*D*x = x + Teuchos::RCP Outer = + Teuchos::rcp(new Iterative_Inverse_Operator(n, 1, C, "Belos (inv(C)=inv(inv(D)*B))", true)); + + // should return x=1/4 + (*Inner)(X,Y); + + if (pid==0) { + std::cout << std::endl << "Vector Y should have all entries [1/4, 1/4, 1/4, ..., 1/4]" << std::endl; + } + Y.Print(std::cout); + + // should return x=1 + (*Outer)(X,Y); + + if (pid==0) { + std::cout << "Vector Y should have all entries [1, 1, 1, ..., 1]" << std::endl; + } + Y.Print(std::cout); + + // Compute the norm of Y - 1.0 + std::vector norm_Y(Y.NumVectors()); + Y.Update(-1.0, X, 1.0); + Y.Norm2(&norm_Y[0]); + + if (pid==0) + std::cout << "Two-norm of std::vector (Y-1.0) : "<< norm_Y[0] << std::endl; + + success = (norm_Y[0] < 1e-10 && !Teuchos::ScalarTraits::isnaninf( norm_Y[0] ) ); + + if (success) { + if (pid==0) + std::cout << "End Result: TEST PASSED" << std::endl; + } else { + if (pid==0) + std::cout << "End Result: TEST FAILED" << std::endl; + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/packages/belos/src/BelosCGIter.hpp b/packages/belos/src/BelosCGIter.hpp index 4686b64cee18..d58f994314b2 100644 --- a/packages/belos/src/BelosCGIter.hpp +++ b/packages/belos/src/BelosCGIter.hpp @@ -399,10 +399,6 @@ class CGIter : virtual public CGIteration { // std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width."); - // Create convenience variables for zero and one. - const ScalarType one = Teuchos::ScalarTraits::one(); - const ScalarType zero = Teuchos::ScalarTraits::zero(); - if (newstate.R != Teuchos::null) { TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_), @@ -413,7 +409,7 @@ class CGIter : virtual public CGIteration { // Copy basis vectors from newstate into V if (newstate.R != R_) { // copy over the initial residual (unpreconditioned). - MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ ); + MVT::Assign( *newstate.R, *R_ ); } // Compute initial direction vectors @@ -432,7 +428,7 @@ class CGIter : virtual public CGIteration { else { MVT::Assign( *R_, *Z_ ); } - MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ ); + MVT::Assign( *Z_, *P_ ); } else { @@ -484,7 +480,7 @@ class CGIter : virtual public CGIteration { rHz[0] = rHs(0,1); } else MVT::MvDot( *R_, *Z_, rHz ); - + //////////////////////////////////////////////////////////////// // Iterate until the status test tells us to stop. // diff --git a/packages/belos/src/BelosMinresIter.hpp b/packages/belos/src/BelosMinresIter.hpp index 0153b7fab3d0..5cc4a452c6ce 100644 --- a/packages/belos/src/BelosMinresIter.hpp +++ b/packages/belos/src/BelosMinresIter.hpp @@ -403,12 +403,13 @@ class MinresIter : virtual public MinresIteration { // Create convenience variables for zero, one. const ScalarType one = SCT::one(); const MagnitudeType zero = SMT::zero(); + const MagnitudeType m_zero = SMT::zero(); // Set up y and v for the first Lanczos vector v_1. // y = beta1_ P' v1, where P = C**(-1). // v is really P' v1. - MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R2_ ); - MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R1_ ); + MVT::Assign( *newstate.Y, *R2_ ); + MVT::Assign( *newstate.Y, *R1_ ); // Initialize the W's to 0. MVT::MvInit ( *W_ ); @@ -420,7 +421,7 @@ class MinresIter : virtual public MinresIteration { else { if (newstate.Y != Y_) { // copy over the initial residual (unpreconditioned). - MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *Y_ ); + MVT::Assign( *newstate.Y, *Y_ ); } } @@ -428,7 +429,7 @@ class MinresIter : virtual public MinresIteration { beta1_ = Teuchos::SerialDenseMatrix( 1, 1 ); MVT::MvTransMv( one, *newstate.Y, *Y_, beta1_ ); - TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < zero, + TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < m_zero, std::invalid_argument, "The preconditioner is not positive definite." ); @@ -462,7 +463,8 @@ class MinresIter : virtual public MinresIteration { // Create convenience variables for zero and one. const ScalarType one = SCT::one(); - const MagnitudeType zero = SMT::zero(); + const ScalarType zero = SCT::zero(); + const MagnitudeType m_zero = SMT::zero(); // Allocate memory for scalars. Teuchos::SerialDenseMatrix alpha( 1, 1 ); @@ -553,10 +555,10 @@ class MinresIter : virtual public MinresIteration { // algebra library to compute a posteriori rounding error bounds // for the inner product, and then changing // Belos::MultiVecTraits to make this information available). - TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) <= zero, + TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) < m_zero, MinresIterateFailure, - "Belos::MinresIter::iterate(): Encountered nonpositi" - "ve value " << beta(0,0) << " for r2^H*M*r2 at itera" + "Belos::MinresIter::iterate(): Encountered negative " + "value " << beta(0,0) << " for r2^H*M*r2 at itera" "tion " << iter_ << ": MINRES cannot continue." ); beta(0,0) = SCT::squareroot( beta(0,0) ); @@ -579,7 +581,7 @@ class MinresIter : virtual public MinresIteration { // w1 = w2; // w2 = w; - MVT::MvAddMv( one, *W_, zero, *W_, *W1_ ); + MVT::Assign( *W_, *W1_ ); tmpW = W1_; W1_ = W2_; W2_ = W_; @@ -592,8 +594,8 @@ class MinresIter : virtual public MinresIteration { // Update x: // x = x + phi*w; - MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec ); - lp_->updateSolution(); + //MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec ); + lp_->updateSolution( W_, true, phi ); } // end while (sTest_->checkStatus(this) != Passed) } diff --git a/packages/belos/tpetra/test/BiCGStab/CMakeLists.txt b/packages/belos/tpetra/test/BiCGStab/CMakeLists.txt new file mode 100644 index 000000000000..776ace7b12f5 --- /dev/null +++ b/packages/belos/tpetra/test/BiCGStab/CMakeLists.txt @@ -0,0 +1,21 @@ +ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) + +IF (${PACKAGE_NAME}_ENABLE_Triutils) + + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_BiCGStab_hb + SOURCES test_bicgstab_hb.cpp + COMM serial mpi + ARGS "--verbose --max-iters=10" + "--verbose --max-iters=10 --tol=1e-10" + "--verbose --max-iters=10 --num-rhs=2" + STANDARD_PASS_OUTPUT + ) + + TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestBiCGSTABFiles + SOURCE_DIR ${${PACKAGE_NAME}_SOURCE_DIR}/epetra/test/BiCGStab + SOURCE_FILES cage4.hb + EXEDEPS Tpetra_BiCGStab_hb + ) + +ENDIF() diff --git a/packages/belos/tpetra/test/BiCGStab/test_bicgstab_hb.cpp b/packages/belos/tpetra/test/BiCGStab/test_bicgstab_hb.cpp new file mode 100644 index 000000000000..4479a36f327a --- /dev/null +++ b/packages/belos/tpetra/test/BiCGStab/test_bicgstab_hb.cpp @@ -0,0 +1,219 @@ +//@HEADER +// ************************************************************************ +// +// Belos: Block Linear Solvers Package +// Copyright 2004 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// ************************************************************************ +//@HEADER +// +// This driver reads a problem from a file, which can only be a Harwell-Boeing (*.hb) +// matrix. The right-hand side is generated using a random solution vector that has +// the matrix applied to it. The initial guesses are all set to zero. +// +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraAdapter.hpp" +#include "BelosBiCGStabSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" + +#include +#include +#include +#include "Teuchos_StandardCatchMacros.hpp" +#include +#include + +int main (int argc, char *argv[]) +{ + using Teuchos::ParameterList; + using Teuchos::RCP; + using Teuchos::rcp; + using std::cout; + using std::endl; + + typedef double ST; + typedef Teuchos::ScalarTraits SCT; + typedef SCT::magnitudeType MT; + typedef Tpetra::Operator OP; + typedef Tpetra::MultiVector MV; + typedef Belos::OperatorTraits OPT; + typedef Belos::MultiVecTraits MVT; + + Teuchos::GlobalMPISession mpisess(&argc,&argv,&cout); + + const ST one = SCT::one(); + + RCP > comm = Tpetra::getDefaultComm(); + int MyPID = rank(*comm); + + bool verbose = false; + bool success = true; + + // This "try" relates to TEUCHOS_STANDARD_CATCH_STATEMENTS near the + // bottom of main(). That macro has the corresponding "catch". + try { + bool proc_verbose = false; + int frequency = -1; // frequency of status test output. + int numrhs = 1; // number of right-hand sides to solve for + int maxiters = -1; // maximum number of iterations allowed per linear system + std::string filename ("cage4.hb"); + MT tol = 1.0e-5; // relative residual tolerance + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption ("verbose", "quiet", &verbose, "Whether to print messages " + "and results."); + cmdp.setOption ("frequency", &frequency, "Frequency of solver output " + "(1 means every iteration; -1 means never)."); + cmdp.setOption ("filename", &filename, "Test matrix filename. " + "Allowed file extensions: *.hb, *.mtx, *.triU, *.triS"); + cmdp.setOption ("tol", &tol, "Relative residual tolerance for solver."); + cmdp.setOption ("num-rhs", &numrhs, "Number of right-hand sides to solve."); + cmdp.setOption ("max-iters", &maxiters, "Maximum number of iterations per " + "linear system (-1 = adapted to problem/block size)."); + if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (! verbose) { + frequency = -1; // reset frequency if test is not verbose + } + + proc_verbose = ( verbose && (MyPID==0) ); + + // + // *************Get the problem********************* + // + Belos::Tpetra::HarwellBoeingReader > reader( comm ); + RCP > A = reader.readFromFile( filename ); + RCP > map = A->getDomainMap(); + + // Create initial vectors + RCP B, X; + X = rcp( new MV(map,numrhs) ); + MVT::MvRandom( *X ); + B = rcp( new MV(map,numrhs) ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + + // + // Create parameter list for the Belos solver + // + const int NumGlobalElements = B->getGlobalLength (); + if (maxiters == -1) { + maxiters = NumGlobalElements - 1; // maximum number of iterations to run + } + RCP belosList (new ParameterList ("Belos")); + belosList->set ("Maximum Iterations", maxiters); + belosList->set ("Convergence Tolerance", tol); + if (verbose) { + belosList->set ("Verbosity", Belos::Errors + Belos::Warnings + + Belos::TimingDetails + Belos::StatusTestDetails); + if (frequency > 0) { + belosList->set ("Output Frequency", frequency); + } + } + else { + belosList->set ("Verbosity", Belos::Errors + Belos::Warnings + + Belos::FinalSummary); + } + // + // Construct a preconditioned linear problem + // + RCP > problem + = rcp (new Belos::LinearProblem (A, X, B)); + bool set = problem->setProblem (); + if (! set) { + if (proc_verbose) { + cout << endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << endl; + } + return -1; + } + + // Create a Belos solver. + RCP > solver + = rcp (new Belos::BiCGStabSolMgr (problem, belosList)); + + if (proc_verbose) { + cout << endl << endl; + cout << "Dimension of matrix: " << NumGlobalElements << endl; + cout << "Number of right-hand sides: " << numrhs << endl; + cout << "Max number of BiCGStab iterations: " << maxiters << endl; + cout << "Relative residual tolerance: " << tol << endl; + cout << endl; + } + + // Ask Belos to solve the linear system. + Belos::ReturnType ret = solver->solve(); + + // + // Compute actual residuals. + // + bool badRes = false; + std::vector actual_resids( numrhs ); + std::vector rhs_norm( numrhs ); + MV resid(map, numrhs); + OPT::Apply( *A, *X, resid ); + MVT::MvAddMv( -one, resid, one, *B, resid ); + MVT::MvNorm( resid, actual_resids ); + MVT::MvNorm( *B, rhs_norm ); + if (proc_verbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; + } + + if ( ret!=Belos::Converged || badRes) { + if (proc_verbose) { + std::cout << "\nEnd Result: TEST FAILED" << std::endl; + } + return -1; + } + // + // Default return value + // + if (proc_verbose) { + std::cout << "\nEnd Result: TEST PASSED" << std::endl; + } + + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/packages/belos/tpetra/test/CMakeLists.txt b/packages/belos/tpetra/test/CMakeLists.txt index 42b36e7d553e..71bb66a37f74 100644 --- a/packages/belos/tpetra/test/CMakeLists.txt +++ b/packages/belos/tpetra/test/CMakeLists.txt @@ -1,7 +1,10 @@ ADD_SUBDIRECTORY(GCRODR) +ADD_SUBDIRECTORY(RCG) +ADD_SUBDIRECTORY(BiCGStab) ADD_SUBDIRECTORY(BlockCG) -ADD_SUBDIRECTORY(FixedPoint) ADD_SUBDIRECTORY(BlockGmres) +ADD_SUBDIRECTORY(TFQMR) +ADD_SUBDIRECTORY(FixedPoint) ADD_SUBDIRECTORY(LinearSolverFactory) ADD_SUBDIRECTORY(MultipleSolves) ADD_SUBDIRECTORY(MVOPTester) diff --git a/packages/belos/tpetra/test/LinearSolverFactory/SolverFactory.cpp b/packages/belos/tpetra/test/LinearSolverFactory/SolverFactory.cpp index 69c2350ed5bc..f7483a2c084b 100644 --- a/packages/belos/tpetra/test/LinearSolverFactory/SolverFactory.cpp +++ b/packages/belos/tpetra/test/LinearSolverFactory/SolverFactory.cpp @@ -159,8 +159,8 @@ testSolver (Teuchos::FancyOStream& out, Teuchos::OSTab tab1 (out); RCP > solver; + Belos::SolverFactory factory; try { - Belos::SolverFactory factory; solver = factory.create (solverName, Teuchos::null); } catch (std::exception& e) { out << "*** FAILED: Belos::SolverFactory::create threw an exception: " @@ -178,6 +178,7 @@ testSolver (Teuchos::FancyOStream& out, out << "Create the Belos::LinearProblem to solve" << endl; typedef Belos::LinearProblem linear_problem_type; + X->putScalar (STS::zero ()); RCP problem (new linear_problem_type (A, X, B)); problem->setProblem (); @@ -186,8 +187,14 @@ testSolver (Teuchos::FancyOStream& out, out << "Apply solver to \"solve\" AX=B for X. Belos already has solver " "tests; the point is to check that solve() doesn't throw." << endl; - X->putScalar (STS::zero ()); - solver->solve (); + try { + solver->solve (); + } catch (std::exception& e) { + out << "*** FAILED: Belos::SolverFactory::solve threw an exception: " + << e.what () << endl; + success = false; + return; + } } template @@ -290,8 +297,18 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL( SolverFactory, CreateAndSolve, SC, LO, GO, NT // // Teuchos::Array solverNames = factory.supportedSolverNames (); // const int numSolvers = static_cast (solverNames.size ()); - const char* solverNames[2] = {"CG", "GMRES"}; - const int numSolvers = 2; + const char* solverNames[10] = { + "BICGSTAB", + "BLOCK CG", + "BLOCK GMRES", + "FIXED POINT", + "GCRODR", + "MINRES", + "PSEUDOBLOCK CG", + "PSEUDOBLOCK GMRES", + "PSEUDOBLOCK TFQMR", + "TFQMR"}; + const int numSolvers = 10; int numSolversTested = 0; for (int k = 0; k < numSolvers; ++k) { diff --git a/packages/belos/tpetra/test/RCG/CMakeLists.txt b/packages/belos/tpetra/test/RCG/CMakeLists.txt new file mode 100644 index 000000000000..6623e64cbb0b --- /dev/null +++ b/packages/belos/tpetra/test/RCG/CMakeLists.txt @@ -0,0 +1,21 @@ +ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) + +IF (${PACKAGE_NAME}_ENABLE_Triutils) + +TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_rcg_hb_test + SOURCES test_rcg_hb.cpp + ARGS + "--verbose --tol=1e-6 --filename=bcsstk14.hb --num-rhs=3 --max-subspace=100 --recycle=10 --max-iters=4000" + COMM serial mpi + ) + +ASSERT_DEFINED(Anasazi_SOURCE_DIR) +TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestRCGFiles + SOURCE_DIR ${Anasazi_SOURCE_DIR}/testmatrices + SOURCE_FILES bcsstk14.hb + EXEDEPS Tpetra_rcg_hb_test + ) + +ENDIF() + diff --git a/packages/belos/tpetra/test/RCG/test_rcg_hb.cpp b/packages/belos/tpetra/test/RCG/test_rcg_hb.cpp new file mode 100644 index 000000000000..e66044df42df --- /dev/null +++ b/packages/belos/tpetra/test/RCG/test_rcg_hb.cpp @@ -0,0 +1,232 @@ +//@HEADER +// ************************************************************************ +// +// Belos: Block Linear Solvers Package +// Copyright 2004 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// ************************************************************************ +//@HEADER +// +// This driver reads a problem from a Harwell-Boeing (HB) file. +// The right-hand-side corresponds to a randomly generated solution. +// The initial guesses are all set to zero. +// +// NOTE: No preconditioner is used in this case. +// +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraAdapter.hpp" +#include "BelosRCGSolMgr.hpp" + +#include +#include +#include +#include +#include + +// I/O for Harwell-Boeing files +#define HIDE_TPETRA_INOUT_IMPLEMENTATIONS +#include + +using namespace Teuchos; +using Tpetra::Operator; +using Tpetra::CrsMatrix; +using Tpetra::MultiVector; +using std::endl; +using std::cout; +using std::vector; +using Teuchos::tuple; + +int main(int argc, char *argv[]) { + + typedef double ST; + typedef ScalarTraits SCT; + typedef SCT::magnitudeType MT; + typedef Tpetra::Operator OP; + typedef Tpetra::MultiVector MV; + typedef Belos::OperatorTraits OPT; + typedef Belos::MultiVecTraits MVT; + + GlobalMPISession mpisess(&argc,&argv,&cout); + + const ST one = SCT::one(); + + RCP > comm = Tpetra::getDefaultComm(); + int MyPID = rank(*comm); + + // + // Get test parameters from command-line processor + // + bool verbose = false, proc_verbose = false, debug = false; + int frequency = -1; // how often residuals are printed by solver + std::string filename("bcsstk14.hb"); + int numBlocks = 30; // maximum number of blocks the solver can use for the Krylov subspace + int recycleBlocks = 3; // maximum number of blocks the solver can use for the recycle space + int numrhs = 1; // number of right-hand sides to solve for + int maxiters = -1; // maximum number of iterations allowed per linear system + MT tol = 1.0e-5; // relative residual tolerance + + CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + cmdp.setOption("debug","nodebug",&debug,"Run debugging checks."); + cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by the RCG solver."); + cmdp.setOption("max-subspace",&numBlocks,"Maximum number of vectors in search space (not including recycle space)."); + cmdp.setOption("recycle",&recycleBlocks,"Number of vectors in recycle space."); + cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for."); + cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size)."); + cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); + if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (debug) { + verbose = true; + } + if (!verbose) { + frequency = -1; // reset frequency if test is not verbose + } + + proc_verbose = ( verbose && (MyPID==0) ); + + if (proc_verbose) { + std::cout << Belos::Belos_Version() << std::endl << std::endl; + } + + RCP > A; + Tpetra::Utils::readHBMatrix(filename,comm,A); + RCP > map = A->getDomainMap(); + + // Create initial vectors + RCP B, X; + X = rcp( new MV(map,numrhs) ); + MVT::MvRandom( *X ); + B = rcp( new MV(map,numrhs) ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + + // + // ********Other information used by block solver*********** + // *****************(can be user specified)****************** + // + int numIters1; + const int NumGlobalElements = B->getGlobalLength(); + if (maxiters == -1) + maxiters = NumGlobalElements - 1; // maximum number of iterations to run + // + ParameterList belosList; + belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed + belosList.set( "Num Blocks", numBlocks); // Maximum number of blocks in Krylov space + belosList.set( "Num Recycled Blocks", recycleBlocks ); // Number of vectors in recycle space + belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + if (verbose) { + belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + + Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails ); + if (frequency > 0) + belosList.set( "Output Frequency", frequency ); + } + else + belosList.set( "Verbosity", Belos::Errors + Belos::Warnings ); + // + // Construct an unpreconditioned linear problem instance. + // + Belos::LinearProblem problem( A, X, B ); + bool set = problem.setProblem(); + if (set == false) { + if (proc_verbose) + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + // + // ******************************************************************* + // *********************Start the RCG iteration*********************** + // ******************************************************************* + // + Belos::RCGSolMgr solver( rcpFromRef(problem), rcpFromRef(belosList) ); + + // + // **********Print out information about problem******************* + // + if (proc_verbose) { + std::cout << std::endl << std::endl; + std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl; + std::cout << "Number of right-hand sides: " << numrhs << std::endl; + std::cout << "Max number of RCG iterations: " << maxiters << std::endl; + std::cout << "Relative residual tolerance: " << tol << std::endl; + std::cout << std::endl; + } + // + // Perform solve + // + Belos::ReturnType ret = solver.solve(); + numIters1=solver.getNumIters(); + // + // Compute actual residuals. + // + bool badRes = false; + std::vector actual_resids( numrhs ); + std::vector rhs_norm( numrhs ); + MV resid(map, numrhs); + OPT::Apply( *A, *X, resid ); + MVT::MvAddMv( -one, resid, one, *B, resid ); + MVT::MvNorm( resid, actual_resids ); + MVT::MvNorm( *B, rhs_norm ); + if (proc_verbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; + } + + if (proc_verbose) { std::cout << "Solve took " << numIters1 << " iterations." << std::endl; } + + if ( ret!=Belos::Converged || badRes ) { + if (proc_verbose) { + std::cout << "\nEnd Result: TEST FAILED" << std::endl; + } + return -1; + } + // + // Default return value + // + if (proc_verbose) { + std::cout << "\nEnd Result: TEST PASSED" << std::endl; + } + return 0; + // +} // end test_rcg_hb.cpp diff --git a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt new file mode 100644 index 000000000000..ada8a971ea85 --- /dev/null +++ b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt @@ -0,0 +1,22 @@ +ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) +IF (${PACKAGE_NAME}_ENABLE_Triutils) + + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_tfqmr_hb + SOURCES test_tfqmr_hb.cpp + COMM serial mpi + ARGS + "--verbose --filename=orsirr1_scaled.hb" + "--verbose --pseudo --filename=orsirr1_scaled.hb" + "--verbose --explicit --filename=orsirr1_scaled.hb" + "--verbose --recursive --filename=orsirr1_scaled.hb" + STANDARD_PASS_OUTPUT + ) + + TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestTFQMRFiles + SOURCE_DIR ${Belos_SOURCE_DIR}/epetra/example/BlockGmres + SOURCE_FILES orsirr1.hb orsirr1_scaled.hb + EXEDEPS Tpetra_tfqmr_hb + ) + +ENDIF(${PACKAGE_NAME}_ENABLE_Triutils) diff --git a/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp b/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp new file mode 100644 index 000000000000..bb4363d680ac --- /dev/null +++ b/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp @@ -0,0 +1,231 @@ +//@HEADER +// ************************************************************************ +// +// Belos: Block Linear Solvers Package +// Copyright 2004 Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Michael A. Heroux (maherou@sandia.gov) +// +// ************************************************************************ +//@HEADER +// +// This driver reads a problem from a Harwell-Boeing (HB) file. +// The initial guesses are all set to zero. +// +// NOTE: No preconditioner is used in this case. +// +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraAdapter.hpp" +#include "BelosTFQMRSolMgr.hpp" +#include "BelosPseudoBlockTFQMRSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" + +#include +#include +#include +#include +#include + + +int main(int argc, char *argv[]) { + // + typedef double ST; + typedef Teuchos::ScalarTraits SCT; + typedef SCT::magnitudeType MT; + typedef Tpetra::Operator OP; + typedef Tpetra::MultiVector MV; + typedef Belos::OperatorTraits OPT; + typedef Belos::MultiVecTraits MVT; + + using Teuchos::ParameterList; + using Teuchos::RCP; + using Teuchos::rcp; + + const ST one = SCT::one(); + + Teuchos::GlobalMPISession session(&argc, &argv, NULL); + + RCP > comm = Tpetra::getDefaultComm(); + int MyPID = rank(*comm); + + // + bool success = false; + bool verbose = false; + try { + // + // Get test parameters from command-line processor + // + bool proc_verbose = false; + bool explicit_test = false; + bool comp_recursive = false; + bool pseudo = false; + int frequency = -1; // how often residuals are printed by solver + int numrhs = 1; // total number of right-hand sides to solve for + int maxiters = -1; // maximum number of iterations for solver to use + std::string filename("orsirr1.hb"); + double tol = 1.0e-5; // relative residual tolerance + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + + cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); + cmdp.setOption("explicit","implicit-only",&explicit_test,"Compute explicit residuals."); + cmdp.setOption("recursive","native",&comp_recursive,"Compute recursive residuals."); + cmdp.setOption("pseudo","not-pseudo",&pseudo,"Use pseudo-block TFQMR solver."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by TFQMR solver."); + cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); + cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for."); + cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 := adapted to problem size)."); + if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (!verbose) + frequency = -1; // reset frequency if test is not verbose + // + // Get the problem + // + Belos::Tpetra::HarwellBoeingReader > reader( comm ); + RCP > A = reader.readFromFile( filename ); + RCP > map = A->getDomainMap(); + + // Create initial vectors + RCP B, X; + X = rcp( new MV(map,numrhs) ); + MVT::MvRandom( *X ); + B = rcp( new MV(map,numrhs) ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + + proc_verbose = ( verbose && (MyPID==0) ); + // + // Solve using Belos + // + const int NumGlobalElements = B->getGlobalLength(); + if (maxiters == -1) + maxiters = NumGlobalElements; // maximum number of iterations to run + // + ParameterList belosList; + belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed + belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + if (explicit_test) + { + belosList.set( "Explicit Residual Test", true ); // Scale by norm of right-hand side vector." + belosList.set( "Explicit Residual Scaling", "Norm of RHS" ); // Scale by norm of right-hand side vector." + } + if (comp_recursive) + { + belosList.set( "Compute Recursive Residuals", true ); + } + if (verbose) { + belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + + Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails ); + if (frequency > 0) + belosList.set( "Output Frequency", frequency ); + } + else + belosList.set( "Verbosity", Belos::Errors + Belos::Warnings ); + // + // Construct an unpreconditioned linear problem instance. + // + Belos::LinearProblem problem( A, X, B ); + bool set = problem.setProblem(); + if (set == false) { + if (proc_verbose) + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + // + // Create an iterative solver manager. + // + RCP< Belos::SolverManager > solver; + + if (pseudo) + solver = rcp( new Belos::PseudoBlockTFQMRSolMgr(rcp(&problem,false), rcp(&belosList,false)) ); + else + solver = rcp( new Belos::TFQMRSolMgr(rcp(&problem,false), rcp(&belosList,false)) ); + // + // **********Print out information about problem******************* + // + if (proc_verbose) { + std::cout << std::endl << std::endl; + std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl; + std::cout << "Number of right-hand sides: " << numrhs << std::endl; + std::cout << "Max number of TFQMR iterations: " << maxiters << std::endl; + std::cout << "Relative residual tolerance: " << tol << std::endl; + std::cout << std::endl; + } + // + // Perform solve + // + Belos::ReturnType ret = solver->solve(); + // + // Compute actual residuals. + // + bool badRes = false; + std::vector actual_resids( numrhs ); + std::vector rhs_norm( numrhs ); + MV resid(map, numrhs); + OPT::Apply( *A, *X, resid ); + MVT::MvAddMv( -one, resid, one, *B, resid ); + MVT::MvNorm( resid, actual_resids ); + MVT::MvNorm( *B, rhs_norm ); + if (proc_verbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; + } + + if ( ret!=Belos::Converged || badRes) { + if (proc_verbose) { + std::cout << "\nEnd Result: TEST FAILED" << std::endl; + } + return -1; + } + // + // Default return value + // + if (proc_verbose) { + std::cout << "\nEnd Result: TEST PASSED" << std::endl; + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); + + return ( success ? EXIT_SUCCESS : EXIT_FAILURE ); + +} // end test_tfqmr_hb.cpp