From 06042f32e1b0bdf96f741ac061db2af263e8fabd Mon Sep 17 00:00:00 2001 From: "Heidi K. Thornquist" Date: Fri, 21 May 2021 15:21:17 -0600 Subject: [PATCH] Adds more un-preconditioned Tpetra tests to Belos This commit translates more of the Epetra-based tests to Tpetra linear algebra. 1) The Tpetra solver factory test was extended to include all of the linear solvers that are accessible through the Belos linear solver factory. This resulted in the detection of a bug in the test that was causing solver failures, which has been addressed. 2) Adds diagonal test to MINRES and TFQMR solvers for Epetra. The extension of the Tpetra solver factory test illustrated a problem in the MINRES solver for diagonal/identity matrices, which has been addressed in this commit. 3) Adds BiCGStab, RCG and TFQMR test to Tpetra-based linear algebra. Testing for these solvers has been absent before the extention to the Tpetra solver factory and the addition of these tests. --- .../belos/epetra/test/MINRES/CMakeLists.txt | 7 + .../epetra/test/MINRES/test_minres_diag.cpp | 422 ++++++++++++++++++ .../epetra/test/MINRES/test_minres_hb.cpp | 20 +- .../belos/epetra/test/TFQMR/CMakeLists.txt | 6 + .../epetra/test/TFQMR/test_tfqmr_diag.cpp | 422 ++++++++++++++++++ packages/belos/src/BelosCGIter.hpp | 10 +- packages/belos/src/BelosMinresIter.hpp | 24 +- .../belos/tpetra/test/BiCGStab/CMakeLists.txt | 21 + .../tpetra/test/BiCGStab/test_bicgstab_hb.cpp | 219 +++++++++ packages/belos/tpetra/test/CMakeLists.txt | 5 +- .../LinearSolverFactory/SolverFactory.cpp | 27 +- packages/belos/tpetra/test/RCG/CMakeLists.txt | 21 + .../belos/tpetra/test/RCG/test_rcg_hb.cpp | 232 ++++++++++ .../belos/tpetra/test/TFQMR/CMakeLists.txt | 22 + .../belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp | 231 ++++++++++ 15 files changed, 1659 insertions(+), 30 deletions(-) create mode 100644 packages/belos/epetra/test/MINRES/test_minres_diag.cpp create mode 100644 packages/belos/epetra/test/TFQMR/test_tfqmr_diag.cpp create mode 100644 packages/belos/tpetra/test/BiCGStab/CMakeLists.txt create mode 100644 packages/belos/tpetra/test/BiCGStab/test_bicgstab_hb.cpp create mode 100644 packages/belos/tpetra/test/RCG/CMakeLists.txt create mode 100644 packages/belos/tpetra/test/RCG/test_rcg_hb.cpp create mode 100644 packages/belos/tpetra/test/TFQMR/CMakeLists.txt create mode 100644 packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp 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