From 34b5b55b0980bbb8276da37ce979460d9031bc1d Mon Sep 17 00:00:00 2001 From: "antoine.meyer1" Date: Wed, 23 Aug 2023 05:54:39 -0400 Subject: [PATCH 1/3] Belos: add Tpetra TFQMR test --- .../belos/tpetra/test/TFQMR/CMakeLists.txt | 53 ++- .../tpetra/test/TFQMR/test_ptfqmr_hb.cpp | 282 ++++++++++++ .../tpetra/test/TFQMR/test_tfqmr_diag.cpp | 428 ++++++++++++++++++ .../belos/tpetra/test/TFQMR/test_tfqmr_hb.cpp | 180 ++++---- 4 files changed, 853 insertions(+), 90 deletions(-) create mode 100644 packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp create mode 100644 packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp diff --git a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt index ada8a971ea85..78eecf1e3d39 100644 --- a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt +++ b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt @@ -1,18 +1,63 @@ + +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 ) + #ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Ifpack2) + #IF (${PACKAGE_NAME}_ENABLE_Ifpack2) + + #TRIBITS_ADD_EXECUTABLE_AND_TEST( + # Tpetra_ptfqmr_hb + # SOURCES test_ptfqmr_hb.cpp + # COMM serial mpi + # ARGS + # "--verbose --not-pseudo --left-prec" + # "--verbose --not-pseudo --left-prec --num-rhs=10" + # "--verbose --not-pseudo --right-prec" + # "--verbose --not-pseudo --right-prec --num-rhs=10" + # STANDARD_PASS_OUTPUT + #) + + #TRIBITS_ADD_EXECUTABLE_AND_TEST( + # Tpetra_pseudo_ptfqmr_hb + # SOURCES test_ptfqmr_hb.cpp + # COMM serial mpi + # ARGS + # "--verbose --pseudo --left-prec" + # "--verbose --pseudo --left-prec --num-rhs=10" + # "--verbose --pseudo --right-prec" + # "--verbose --pseudo --right-prec --num-rhs=10" + # STANDARD_PASS_OUTPUT + #) + + #ENDIF(${PACKAGE_NAME}_ENABLE_Ifpack2) + TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestTFQMRFiles SOURCE_DIR ${Belos_SOURCE_DIR}/epetra/example/BlockGmres SOURCE_FILES orsirr1.hb orsirr1_scaled.hb diff --git a/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp b/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp new file mode 100644 index 000000000000..fd0ad6e39507 --- /dev/null +++ b/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp @@ -0,0 +1,282 @@ +//@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. +// Multiple right-hand-sides are created randomly. +// The initial guesses are all set to zero. +// + +// Tpetra +#include +#include +#include + +// Teuchos +#include +#include +#include +#include + +// Ifpack2 +#include +//#include // AM: todo, finds an equivalent + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraAdapter.hpp" +#include "BelosPseudoBlockTFQMRSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" + +template +int run(int argc, char *argv[]) { + // Teuchos + using Teuchos::ParameterList; + using Teuchos::RCP; + using Teuchos::rcp; + + // 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; + + using OP = Tpetra::Operator; + using MV = Tpetra::MultiVector; + + using tcrsmatrix_t = Tpetra::CrsMatrix; + using tmap_t = Tpetra::Map; + + // 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); + + bool success = false; + bool verbose = false; + + try { + bool procVerbose = false; + bool leftprec = true; // use left preconditioning to solve these linear systems + bool explicitTest = true; + bool pseudo = false; + int frequency = -1; // how often residuals are printed by solver + int numrhs = 1; + int maxiters = -1; // maximum iterations allowed + std::string filename("orsirr1.hb"); + ST tol = sqrt(std::numeric_limits::epsilon()); // relative residual tolerance + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + cmdp.setOption("left-prec","right-prec",&leftprec,"Left preconditioning or right."); + cmdp.setOption("explicit","implicit-only",&explicitTest,"Compute explicit residuals."); + cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); + cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); + cmdp.setOption("pseudo","not-pseudo",&pseudo,"Use pseudo-block TFQMR solver."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver."); + cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for."); + cmdp.setOption("maxiters",&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(); + procVerbose = verbose && (MyPID==0); /* Only print on zero processor */ + + // *****Construct the Preconditioner***** + if (procVerbose) std::cout << std::endl << std::endl; + if (procVerbose) std::cout << "Constructing ILU preconditioner" << std::endl; + int Lfill = 2; + // if (argc > 2) Lfill = atoi(argv[2]); + if (procVerbose) std::cout << "Using Lfill = " << Lfill << std::endl; + int Overlap = 2; + // if (argc > 3) Overlap = atoi(argv[3]); + if (procVerbose) std::cout << "Using Level Overlap = " << Overlap << std::endl; + double Athresh = 0.0; + // if (argc > 4) Athresh = atof(argv[4]); + if (procVerbose) std::cout << "Using Absolute Threshold Value of " << Athresh << std::endl; + double Rthresh = 1.0; + // if (argc >5) Rthresh = atof(argv[5]); + if (procVerbose) std::cout << "Using Relative Threshold Value of " << Rthresh << std::endl; + + // Init Ifpack2 classes + RCP ilukGraph; + // AM: TODO, IFPACK ISSUE + //RCP ilukFactors; + + // AM: TODO, IFPACK ISSUE + if (Lfill > -1) { + ilukGraph = rcp(new Ifpack2::IlukGraph(A->getGraph(), Lfill, Overlap)); + /*int info = ilukGraph->ConstructFilledGraph(); + TEUCHOS_ASSERT( info == 0 ); + ilukFactors = rcp(new Ifpack_CrsRiluk(*ilukGraph)); + int initerr = ilukFactors->InitValues(*A); + if (initerr != 0) std::cout << "InitValues error = " << initerr; + info = ilukFactors->Factor(); + TEUCHOS_ASSERT( info == 0 );*/ + } + + // AM: FROM HERE, CONVERSION TO DO BELOW + + // + bool transA = false; + double Cond_Est; + ilukFactors->Condest(transA, Cond_Est); + if (procVerbose) { + std::cout << "Condition number estimate for this preconditoner = " << Cond_Est << std::endl; + std::cout << std::endl; + } + + // + // Create the Belos preconditioned operator from the Ifpack preconditioner. + // NOTE: This is necessary because Belos expects an operator to apply the + // preconditioner with Apply() NOT ApplyInverse(). + RCP Prec = rcp( new Belos::EpetraPrecOp( ilukFactors ) ); + + // + // ********Other information used by block solver*********** + // *****************(can be user specified)****************** + // + const int NumGlobalElements = Map.NumGlobalElements(); + 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( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + belosList.set( "Explicit Residual Test", explicit_test ); // Need to check for the explicit residual before returning + 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 ); + + // *****Construct solution std::vector and random right-hand-sides ***** + RCP X = rcp( new MV(Map, numrhs) ); + RCP B = rcp( new MV(Map, numrhs) ); + MVT::MvRandom( *X ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + + Belos::LinearProblem problem( A, X, B ); + if (leftprec) + problem.setLeftPrec( Prec ); + else + problem.setRightPrec( Prec ); + + bool set = problem.setProblem(); + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + + // Start the TFQMR or Pseudo Block TFQMR iteration + 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 (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; + std::cout << "Max number of Pseudo Block 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 actualResids( numrhs ); + std::vector rhsNorm( numrhs ); + MV R(Map, numrhs); + OPT::Apply( *A, *X, R ); + MVT::MvAddMv( -1.0, R, 1.0, *B, R ); + MVT::MvNorm( R, actualResids ); + MVT::MvNorm( *B, rhsNorm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol ) badRes = true; + } + } + + success = ret==Belos::Converged && !badRes; + + if (success) { + if (procVerbose) + std::cout << "End Result: TEST PASSED" << std::endl; + } else { + if (procVerbose) + 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_diag.cpp b/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp new file mode 100644 index 000000000000..57534b8d054a --- /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; + + MPI_Allreduce(&n_in, &nGlobal, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + pComm = Tpetra::getDefaultComm(); + + 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); +} From ac0e8c2fc362a1b0a03c4eff905ef8b130684e47 Mon Sep 17 00:00:00 2001 From: antoinemeyer5 Date: Thu, 21 Sep 2023 14:44:31 +0200 Subject: [PATCH 2/3] Belos: Remove test_ptfqmr_hb because cannot be converted --- .../belos/tpetra/test/TFQMR/CMakeLists.txt | 29 -- .../tpetra/test/TFQMR/test_ptfqmr_hb.cpp | 282 ------------------ 2 files changed, 311 deletions(-) delete mode 100644 packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp diff --git a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt index 78eecf1e3d39..ae4eaf2e73dd 100644 --- a/packages/belos/tpetra/test/TFQMR/CMakeLists.txt +++ b/packages/belos/tpetra/test/TFQMR/CMakeLists.txt @@ -29,35 +29,6 @@ IF (${PACKAGE_NAME}_ENABLE_Triutils) STANDARD_PASS_OUTPUT ) - #ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Ifpack2) - #IF (${PACKAGE_NAME}_ENABLE_Ifpack2) - - #TRIBITS_ADD_EXECUTABLE_AND_TEST( - # Tpetra_ptfqmr_hb - # SOURCES test_ptfqmr_hb.cpp - # COMM serial mpi - # ARGS - # "--verbose --not-pseudo --left-prec" - # "--verbose --not-pseudo --left-prec --num-rhs=10" - # "--verbose --not-pseudo --right-prec" - # "--verbose --not-pseudo --right-prec --num-rhs=10" - # STANDARD_PASS_OUTPUT - #) - - #TRIBITS_ADD_EXECUTABLE_AND_TEST( - # Tpetra_pseudo_ptfqmr_hb - # SOURCES test_ptfqmr_hb.cpp - # COMM serial mpi - # ARGS - # "--verbose --pseudo --left-prec" - # "--verbose --pseudo --left-prec --num-rhs=10" - # "--verbose --pseudo --right-prec" - # "--verbose --pseudo --right-prec --num-rhs=10" - # STANDARD_PASS_OUTPUT - #) - - #ENDIF(${PACKAGE_NAME}_ENABLE_Ifpack2) - TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestTFQMRFiles SOURCE_DIR ${Belos_SOURCE_DIR}/epetra/example/BlockGmres SOURCE_FILES orsirr1.hb orsirr1_scaled.hb diff --git a/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp b/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp deleted file mode 100644 index fd0ad6e39507..000000000000 --- a/packages/belos/tpetra/test/TFQMR/test_ptfqmr_hb.cpp +++ /dev/null @@ -1,282 +0,0 @@ -//@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. -// Multiple right-hand-sides are created randomly. -// The initial guesses are all set to zero. -// - -// Tpetra -#include -#include -#include - -// Teuchos -#include -#include -#include -#include - -// Ifpack2 -#include -//#include // AM: todo, finds an equivalent - -// Belos -#include "BelosConfigDefs.hpp" -#include "BelosLinearProblem.hpp" -#include "BelosTpetraAdapter.hpp" -#include "BelosPseudoBlockTFQMRSolMgr.hpp" -#include "BelosTpetraTestFramework.hpp" - -template -int run(int argc, char *argv[]) { - // Teuchos - using Teuchos::ParameterList; - using Teuchos::RCP; - using Teuchos::rcp; - - // 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; - - using OP = Tpetra::Operator; - using MV = Tpetra::MultiVector; - - using tcrsmatrix_t = Tpetra::CrsMatrix; - using tmap_t = Tpetra::Map; - - // 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); - - bool success = false; - bool verbose = false; - - try { - bool procVerbose = false; - bool leftprec = true; // use left preconditioning to solve these linear systems - bool explicitTest = true; - bool pseudo = false; - int frequency = -1; // how often residuals are printed by solver - int numrhs = 1; - int maxiters = -1; // maximum iterations allowed - std::string filename("orsirr1.hb"); - ST tol = sqrt(std::numeric_limits::epsilon()); // relative residual tolerance - - Teuchos::CommandLineProcessor cmdp(false,true); - cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); - cmdp.setOption("left-prec","right-prec",&leftprec,"Left preconditioning or right."); - cmdp.setOption("explicit","implicit-only",&explicitTest,"Compute explicit residuals."); - cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); - cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); - cmdp.setOption("pseudo","not-pseudo",&pseudo,"Use pseudo-block TFQMR solver."); - cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver."); - cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for."); - cmdp.setOption("maxiters",&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(); - procVerbose = verbose && (MyPID==0); /* Only print on zero processor */ - - // *****Construct the Preconditioner***** - if (procVerbose) std::cout << std::endl << std::endl; - if (procVerbose) std::cout << "Constructing ILU preconditioner" << std::endl; - int Lfill = 2; - // if (argc > 2) Lfill = atoi(argv[2]); - if (procVerbose) std::cout << "Using Lfill = " << Lfill << std::endl; - int Overlap = 2; - // if (argc > 3) Overlap = atoi(argv[3]); - if (procVerbose) std::cout << "Using Level Overlap = " << Overlap << std::endl; - double Athresh = 0.0; - // if (argc > 4) Athresh = atof(argv[4]); - if (procVerbose) std::cout << "Using Absolute Threshold Value of " << Athresh << std::endl; - double Rthresh = 1.0; - // if (argc >5) Rthresh = atof(argv[5]); - if (procVerbose) std::cout << "Using Relative Threshold Value of " << Rthresh << std::endl; - - // Init Ifpack2 classes - RCP ilukGraph; - // AM: TODO, IFPACK ISSUE - //RCP ilukFactors; - - // AM: TODO, IFPACK ISSUE - if (Lfill > -1) { - ilukGraph = rcp(new Ifpack2::IlukGraph(A->getGraph(), Lfill, Overlap)); - /*int info = ilukGraph->ConstructFilledGraph(); - TEUCHOS_ASSERT( info == 0 ); - ilukFactors = rcp(new Ifpack_CrsRiluk(*ilukGraph)); - int initerr = ilukFactors->InitValues(*A); - if (initerr != 0) std::cout << "InitValues error = " << initerr; - info = ilukFactors->Factor(); - TEUCHOS_ASSERT( info == 0 );*/ - } - - // AM: FROM HERE, CONVERSION TO DO BELOW - - // - bool transA = false; - double Cond_Est; - ilukFactors->Condest(transA, Cond_Est); - if (procVerbose) { - std::cout << "Condition number estimate for this preconditoner = " << Cond_Est << std::endl; - std::cout << std::endl; - } - - // - // Create the Belos preconditioned operator from the Ifpack preconditioner. - // NOTE: This is necessary because Belos expects an operator to apply the - // preconditioner with Apply() NOT ApplyInverse(). - RCP Prec = rcp( new Belos::EpetraPrecOp( ilukFactors ) ); - - // - // ********Other information used by block solver*********** - // *****************(can be user specified)****************** - // - const int NumGlobalElements = Map.NumGlobalElements(); - 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( "Convergence Tolerance", tol ); // Relative convergence tolerance requested - belosList.set( "Explicit Residual Test", explicit_test ); // Need to check for the explicit residual before returning - 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 ); - - // *****Construct solution std::vector and random right-hand-sides ***** - RCP X = rcp( new MV(Map, numrhs) ); - RCP B = rcp( new MV(Map, numrhs) ); - MVT::MvRandom( *X ); - OPT::Apply( *A, *X, *B ); - MVT::MvInit( *X, 0.0 ); - - Belos::LinearProblem problem( A, X, B ); - if (leftprec) - problem.setLeftPrec( Prec ); - else - problem.setRightPrec( Prec ); - - bool set = problem.setProblem(); - if (set == false) { - if (procVerbose) - std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; - return -1; - } - - // Start the TFQMR or Pseudo Block TFQMR iteration - 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 (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; - std::cout << "Max number of Pseudo Block 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 actualResids( numrhs ); - std::vector rhsNorm( numrhs ); - MV R(Map, numrhs); - OPT::Apply( *A, *X, R ); - MVT::MvAddMv( -1.0, R, 1.0, *B, R ); - MVT::MvNorm( R, actualResids ); - MVT::MvNorm( *B, rhsNorm ); - if (procVerbose) { - std::cout<< "---------- Actual Residuals (normalized) ----------"< tol ) badRes = true; - } - } - - success = ret==Belos::Converged && !badRes; - - if (success) { - if (procVerbose) - std::cout << "End Result: TEST PASSED" << std::endl; - } else { - if (procVerbose) - 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); -} From 4569681779c03d82ae8ce7b155a4843f6f51c69a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Dutheillet-Lamonth=C3=A9zie?= Date: Fri, 22 Sep 2023 19:45:11 +0200 Subject: [PATCH 3/3] Belos: fix bad MPI call --- packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp b/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp index 57534b8d054a..de0f6bb502be 100644 --- a/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp +++ b/packages/belos/tpetra/test/TFQMR/test_tfqmr_diag.cpp @@ -267,8 +267,8 @@ IterativeInverseOperator::IterativeInverseOperator(int n_in, int { int nGlobal; - MPI_Allreduce(&n_in, &nGlobal, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); 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 ) );