diff --git a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt index ada8a971ea85..ae4eaf2e73dd 100644 --- a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt +++ b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt @@ -1,15 +1,31 @@ + +TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_tfqmr_diag + SOURCES test_tfqmr_diag.cpp + COMM serial mpi + STANDARD_PASS_OUTPUT +) + 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 --not-pseudo --filename=orsirr1_scaled.hb" + "--verbose --not-pseudo --explicit --filename=orsirr1_scaled.hb" + "--verbose --not-pseudo --recursive --filename=orsirr1_scaled.hb" + STANDARD_PASS_OUTPUT + ) + + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_pseudo_tfqmr_hb + SOURCES test_tfqmr_hb.cpp + COMM serial mpi + ARGS "--verbose --pseudo --filename=orsirr1_scaled.hb" - "--verbose --explicit --filename=orsirr1_scaled.hb" - "--verbose --recursive --filename=orsirr1_scaled.hb" STANDARD_PASS_OUTPUT ) diff --git a/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp b/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp new file mode 100644 index 000000000000..de0f6bb502be --- /dev/null +++ b/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp @@ -0,0 +1,428 @@ +//@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. +// + +// Teuchos +#include +#include + +// Tpetra +#include +#include +#include +#include +#include + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraAdapter.hpp" +#include "BelosTpetraOperator.hpp" +#include "BelosTFQMRSolMgr.hpp" + +using namespace Belos; + +//************************************************************************************************ + +template +class VectorOperator +{ + public: + + VectorOperator(int m_in, int n_in) : m(m_in), n(n_in) {}; + + virtual ~VectorOperator() {}; + + virtual void operator () (const MV &x, MV &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. + VectorOperator( const VectorOperator& ): m(0), n(0) {}; + VectorOperator* operator=( const VectorOperator& ) { return NULL; }; + +}; + +//************************************************************************************************ + +template +class DiagonalOperator : public VectorOperator +{ + public: + + DiagonalOperator(int n_in, ST v_in) : VectorOperator(n_in, n_in), v(v_in) { }; + + ~DiagonalOperator() { }; + + void operator () (const MV &x, MV &y) + { + y.scale(v, x); + }; + + private: + + ST v; +}; + +//************************************************************************************************ + +template +class DiagonalOperator2 : public VectorOperator +{ + public: + + DiagonalOperator2(int n_in, int min_gid_in, ST v_in) + : VectorOperator(n_in, n_in), min_gid(min_gid_in), v(v_in) {} + + ~DiagonalOperator2() { }; + + void operator () (const MV &x, MV &y) + { + auto yLocalData = y.getLocalViewHost(Tpetra::Access::ReadWrite); + auto xLocalData = x.getLocalViewHost(Tpetra::Access::ReadOnly); + + for (size_t j = 0; j < x.getNumVectors(); ++j) { + for (size_t i = 0; i < x.getLocalLength(); ++i) { + yLocalData(i, j) = (min_gid + i + 1) * v * xLocalData(i, j); // NOTE: square operator! + } + } + } + + private: + + int min_gid; + ST v; +}; + +//************************************************************************************************ + +template +class ComposedOperator : public VectorOperator +{ + public: + + ComposedOperator(int n, + const Teuchos::RCP>& pA_in, + const Teuchos::RCP>& pB_in); + + virtual ~ComposedOperator() {}; + + virtual void operator () (const MV &x, MV &y); + + private: + + Teuchos::RCP> pA; + Teuchos::RCP> pB; +}; + +template +ComposedOperator::ComposedOperator(int n_in, + const Teuchos::RCP>& pA_in, + const Teuchos::RCP>& pB_in) +: VectorOperator(n_in, n_in), pA(pA_in), pB(pB_in) +{ +} + +template +void ComposedOperator::operator () (const MV &x, MV &y) +{ + MV ytemp(y.getMap(), y.getNumVectors(), false); + (*pB)( x, ytemp ); + (*pA)( ytemp, y ); +} + +//************************************************************************************************ + +template +class TrilinosInterface : public OP +{ + public: + + TrilinosInterface(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) + {} + + void apply (const MV &X, + MV &Y, + Teuchos::ETransp mode = Teuchos::NO_TRANS, + ST alpha = Teuchos::ScalarTraits::one(), + ST beta = Teuchos::ScalarTraits::zero()) const override; + + virtual ~TrilinosInterface() {}; + + bool hasTransposeApply() const {return(use_transpose);}; + + Teuchos::RCP getDomainMap() const override {return pMap; } + + Teuchos::RCP getRangeMap() const override {return pMap; } + + private: + + Teuchos::RCP> pA; + Teuchos::RCP> pComm; + Teuchos::RCP pMap; + + bool use_transpose; +}; + +template +void TrilinosInterface::apply (const MV &X, MV &Y, + Teuchos::ETransp mode, ST alpha, ST beta) const +{ + (*pA)(X,Y); +} + +//************************************************************************************************ + +template +class IterativeInverseOperator : public VectorOperator +{ + public: + + IterativeInverseOperator(int n_in, int blocksize, + const Teuchos::RCP>& pA_in, + std::string opString="Iterative Solver", bool print_in=false); + + virtual ~IterativeInverseOperator() {} + + virtual void operator () (const MV &b, MV &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; +}; + +template +IterativeInverseOperator::IterativeInverseOperator(int n_in, int blocksize, + const Teuchos::RCP>& pA_in, + std::string opString, bool print_in) +: VectorOperator(n_in, n_in), // square operator + pA(pA_in), + print(print_in), + timer(opString) +{ + int nGlobal; + + pComm = Tpetra::getDefaultComm(); + Teuchos::reduceAll(*pComm, Teuchos::REDUCE_SUM, 1, &n_in, &nGlobal); + + pMap = Teuchos::rcp( new MP(nGlobal, n_in, 0, pComm) ); + pPE = Teuchos::rcp( new TrilinosInterface(pA, pComm, pMap ) ); + + pProb = Teuchos::rcp( new LinearProblem() ); + pProb->setOperator( pPE ); + + int max_iter = 100; + ST tol = sqrt(std::numeric_limits::epsilon()); + 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) ); +} + +template +void IterativeInverseOperator::operator () (const MV &b, MV &x) +{ + int pid = pComm->getRank(); + + // 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; + } +} + +//************************************************************************************************ +//************************************************************************************************ + +template +int run(int argc, char *argv[]) +{ + // Get default Tpetra template types + using ST = typename Tpetra::MultiVector::scalar_type; + using LO = typename Tpetra::MultiVector<>::local_ordinal_type; + using GO = typename Tpetra::MultiVector<>::global_ordinal_type; + using NT = typename Tpetra::MultiVector<>::node_type; + + // Init Tpetra types + using OP = typename Tpetra::Operator; + using MV = typename Tpetra::MultiVector; + using MP = typename Tpetra::Map; + + using Teuchos::RCP; + using Teuchos::rcp; + + Teuchos::GlobalMPISession session(&argc, &argv, NULL); + RCP > comm = Tpetra::getDefaultComm(); + int pid = rank(*comm); + + bool verbose = false; + bool success = true; + + ST tol = sqrt(std::numeric_limits::epsilon()); + + try { + int n(10); + int numRHS=1; + + RCP map = RCP(new MP(n, 0, comm)); + MV X(map, numRHS), Y(map, numRHS); + + X.putScalar(1.0); + + // Inner computes inv(D2)*y + RCP> D2 = rcp(new DiagonalOperator2(n, map->getMinGlobalIndex(), 1.0)); + IterativeInverseOperator 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 + RCP> D = rcp(new DiagonalOperator(n, 4.0)); + RCP> Inner = + rcp(new IterativeInverseOperator(n, 1, D, "Belos (inv(D))", false)); + + // Composed_Operator computed inv(D)*B*x + RCP> B = rcp(new DiagonalOperator(n, 4.0)); + RCP> C = rcp(new ComposedOperator(n, Inner, B)); + + // Outer computes inv(C) = inv(inv(D)*B)*x = inv(B)*D*x = x + RCP> Outer = + rcp(new IterativeInverseOperator(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.getNumVectors()); + Y.update(-1.0, X, 1.0); + Y.norm2(norm_Y); + + if (pid==0) + std::cout << "Two-norm of std::vector (Y-1.0) : "<< norm_Y[0] << std::endl; + + success = (norm_Y[0] < tol && !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; +} + +int main(int argc, char *argv[]) { + return run(argc, argv); + + // wrapped with a check: CMake option Trilinos_ENABLE_FLOAT=ON + // run(argc, argv); +} diff --git a/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp b/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp index 1bf0add4730f..48d67bc767d1 100644 --- a/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp +++ b/packages/belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp @@ -40,85 +40,97 @@ //@HEADER // // This driver reads a problem from a Harwell-Boeing (HB) file. +// The right-hand-side from the HB file is used instead of random vectors. // 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 +// Tpetra #include #include +#include +// Teuchos +#include +#include +#include -int main(int argc, char *argv[]) { - // - typedef Tpetra::MultiVector<>::scalar_type 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; +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTFQMRSolMgr.hpp" +#include "BelosPseudoBlockTFQMRSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" +template +int run(int argc, char *argv[]) { + // Teuchos + using SCT = typename Teuchos::ScalarTraits; + using MT = typename SCT::magnitudeType; + using Teuchos::ParameterList; using Teuchos::RCP; using Teuchos::rcp; - const ST one = SCT::one(); + // Tpetra + using ST = typename Tpetra::MultiVector::scalar_type; + using LO = typename Tpetra::MultiVector<>::local_ordinal_type; + using GO = typename Tpetra::MultiVector<>::global_ordinal_type; + using NT = typename Tpetra::MultiVector<>::node_type; - Teuchos::GlobalMPISession session(&argc, &argv, NULL); + using OP = Tpetra::Operator; + using MV = Tpetra::MultiVector; + using tcrsmatrix_t = Tpetra::CrsMatrix; + + // Belos + using OPT = typename Belos::OperatorTraits; + using MVT = typename Belos::MultiVecTraits; + + // MPI session + Teuchos::GlobalMPISession session(&argc, &argv, NULL); RCP > comm = Tpetra::getDefaultComm(); int MyPID = rank(*comm); - - // + + const ST one = SCT::one(); bool success = false; bool verbose = false; - try { - // + + try + { // Get test parameters from command-line processor - // - bool proc_verbose = false; - bool explicit_test = false; - bool comp_recursive = false; + bool procVerbose = false; + bool explicitTest = false; + bool compRecursive = 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 - + MT 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("explicit","implicit-only",&explicitTest,"Compute explicit residuals."); + cmdp.setOption("recursive","native",&compRecursive,"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("tol",&tol,"Relative residual tolerance used by TFQMR or pseudo-block 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)."); + 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 - // + // Get the problem - // - Belos::Tpetra::HarwellBoeingReader > reader( comm ); - RCP > A = reader.readFromFile( filename ); + 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) ); @@ -126,24 +138,23 @@ int main(int argc, char *argv[]) { B = rcp( new MV(map,numrhs) ); OPT::Apply( *A, *X, *B ); MVT::MvInit( *X, 0.0 ); - - proc_verbose = ( verbose && (MyPID==0) ); - // + procVerbose = ( 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) + if (explicitTest) { 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) + if (compRecursive) { belosList.set( "Compute Recursive Residuals", true ); } @@ -155,29 +166,25 @@ int main(int argc, char *argv[]) { } 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) + if (procVerbose) 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)) ); - // + solver = rcp( new Belos::TFQMRSolMgr(rcp(&problem, false), rcp(&belosList, false))); + // **********Print out information about problem******************* - // - if (proc_verbose) { + if (procVerbose) { std::cout << std::endl << std::endl; std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl; std::cout << "Number of right-hand sides: " << numrhs << std::endl; @@ -185,48 +192,49 @@ int main(int argc, char *argv[]) { 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 ); + std::vector actualResids( numrhs ); + std::vector rhsNorm( 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) { + MVT::MvNorm( resid, actualResids ); + MVT::MvNorm( *B, rhsNorm ); + + if (procVerbose) { std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; - } - - if ( ret!=Belos::Converged || badRes) { - if (proc_verbose) { - std::cout << "\nEnd Result: TEST FAILED" << std::endl; + for ( int i=0; i tol) badRes = true; } - return -1; } - // - // Default return value - // - if (proc_verbose) { - std::cout << "\nEnd Result: TEST PASSED" << std::endl; + + success = ret==Belos::Converged && !badRes; + + if (success) { + if (procVerbose) + std::cout << std::endl << "End Result: TEST PASSED" << std::endl; + } else { + if (procVerbose) + std::cout << std::endl << "End Result: TEST FAILED" << std::endl; } - success = true; } TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); - + return ( success ? EXIT_SUCCESS : EXIT_FAILURE ); +} + +int main(int argc, char *argv[]) { + return run(argc,argv); -} // end test_tfqmr_hb.cpp + // wrapped with a check: CMake option Trilinos_ENABLE_FLOAT=ON + // run(argc,argv); +}