diff --git a/packages/belos/tpetra/test/BlockGmres/CMakeLists.txt b/packages/belos/tpetra/test/BlockGmres/CMakeLists.txt index 59b4c9a1dbf5..6ff2f928de82 100644 --- a/packages/belos/tpetra/test/BlockGmres/CMakeLists.txt +++ b/packages/belos/tpetra/test/BlockGmres/CMakeLists.txt @@ -30,6 +30,24 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( COMM serial mpi ) +TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_resolve_gmres_hb + SOURCES test_resolve_gmres_hb.cpp + COMM serial mpi + ARGS + "--verbose --filename=orsirr1.hb" + "--verbose --filename=orsirr1.hb --pseudo" + STANDARD_PASS_OUTPUT + ) + +TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_pseudo_gmres_hb + SOURCES test_pseudo_gmres_hb.cpp + COMM serial mpi + ARGS "--verbose --filename=orsirr1.hb" + STANDARD_PASS_OUTPUT + ) + ASSERT_DEFINED(Belos_ENABLE_TSQR) IF (Belos_ENABLE_TSQR) TRIBITS_ADD_EXECUTABLE_AND_TEST( @@ -63,3 +81,9 @@ TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_BlockGMRES_hb_CopyFiles EXEDEPS Tpetra_BlockGMRES_hb_test ) +TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_BlockGMRES_pseudo_CopyFiles + SOURCE_DIR ${Belos_SOURCE_DIR}/tpetra/example/BlockGmres + SOURCE_FILES orsirr1.hb + EXEDEPS Tpetra_BlockGMRES_hb_test +) + diff --git a/packages/belos/tpetra/test/BlockGmres/Dummy b/packages/belos/tpetra/test/BlockGmres/Dummy deleted file mode 100644 index e69de29bb2d1..000000000000 diff --git a/packages/belos/tpetra/test/BlockGmres/test_pseudo_gmres_hb.cpp b/packages/belos/tpetra/test/BlockGmres/test_pseudo_gmres_hb.cpp new file mode 100644 index 000000000000..b9888874cd56 --- /dev/null +++ b/packages/belos/tpetra/test/BlockGmres/test_pseudo_gmres_hb.cpp @@ -0,0 +1,314 @@ +//@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. +// +// This test is for testing the deflation in the pseudo-block Gmres solver. +// One set of linear systems is solved and then augmented with additional +// linear systems and resolved. The already solved linear systems should be +// deflated immediately, leaving only the augmented systems to be solved. +// +// + +// Tpetra +#include +#include +#include +#include + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosTpetraTestFramework.hpp" +#include "BelosPseudoBlockGmresSolMgr.hpp" + + +template +int run(int argc, char *argv[]) { + + using ST = typename Tpetra::CrsMatrix::scalar_type; + using LO = typename Tpetra::CrsMatrix<>::local_ordinal_type; + using GO = typename Tpetra::CrsMatrix<>::global_ordinal_type; + using NT = typename Tpetra::CrsMatrix<>::node_type; + + using OP = typename Tpetra::Operator; + using MV = typename Tpetra::MultiVector; + using MT = typename Teuchos::ScalarTraits::magnitudeType; + + using tmap_t = Tpetra::Map; + using tcrsmatrix_t = Tpetra::CrsMatrix; + + using MVT = typename Belos::MultiVecTraits; + using OPT = typename Belos::OperatorTraits; + + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::ParameterList; + + Teuchos::GlobalMPISession session(&argc, &argv, nullptr); + + bool verbose = false; + bool success = true; + + try { + + const auto comm = Tpetra::getDefaultComm(); + const int myPID = comm->getRank(); + + bool procVerbose = false; + int frequency = -1; // how often residuals are printed by solver + int initNumRHS = 5; // how many right-hand sides get solved first + int augNumRHS = 10; // how many right-hand sides are augmented to the first group + int maxRestarts = 15; // number of restarts allowed + int length = 100; + int initBlockSize = 5;// blocksize used for the initial pseudo-block GMRES solve + int augBlockSize = 3; // blocksize used for the augmented pseudo-block GMRES solve + int maxIters = -1; // maximum iterations allowed + std::string filename("orsirr1.hb"); + MT tol = 1.0e-5; // relative residual tolerance + MT aug_tol = 1.0e-5; // relative residual tolerance for augmented system + + 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("filename",&filename,"Filename for Harwell-Boeing test matrix."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver."); + cmdp.setOption("aug-tol",&aug_tol,"Relative residual tolerance used by GMRES solver for augmented systems."); + cmdp.setOption("init-num-rhs",&initNumRHS,"Number of right-hand sides to be initially solved for."); + cmdp.setOption("aug-num-rhs",&augNumRHS,"Number of right-hand sides augmenting the initial solve."); + cmdp.setOption("max-restarts",&maxRestarts,"Maximum number of restarts allowed for GMRES solver."); + cmdp.setOption("block-size",&initBlockSize,"Block size used by GMRES for the initial solve."); + cmdp.setOption("aug-block-size",&augBlockSize,"Block size used by GMRES for the augmented solve."); + cmdp.setOption("max-iters",&maxIters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size)."); + cmdp.setOption("subspace-size",&length,"Dimension of Krylov subspace used by GMRES."); + if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (!verbose) + frequency = -1; // reset frequency if test is not verbose + + procVerbose = verbose && (myPID==0); /* Only print on zero processor */ + + // Get the problem + Belos::Tpetra::HarwellBoeingReader reader( comm ); + RCP A = reader.readFromFile( filename ); + RCP Map = A->getDomainMap(); + + // ********Other information used by block solver*********** + // *****************(can be user specified)****************** + + const int numGlobalElements = Map->getGlobalNumElements(); + if (maxIters == -1) + maxIters = numGlobalElements - 1; // maximum number of iterations to run + + ParameterList belosList; + belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization + belosList.set( "Block Size", initBlockSize ); // Blocksize to be used by iterative solver + belosList.set( "Maximum Iterations", maxIters ); // Maximum number of iterations allowed + belosList.set( "Maximum Restarts", maxRestarts ); // Maximum number of restarts allowed + belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + belosList.set( "Deflation Quorum", initBlockSize ); // Number of converged linear systems before deflation + belosList.set( "Timer Label", "Belos Init" ); // Label used by the timers in this solver + 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 initX = rcp( new MV(Map, initNumRHS) ); + RCP initB = rcp( new MV(Map, initNumRHS) ); + MVT::MvRandom( *initX ); + OPT::Apply( *A, *initX, *initB ); + initX->putScalar( 0.0 ); + Belos::LinearProblem initProblem( A, initX, initB ); + initProblem.setLabel("Belos Init"); + + bool set = initProblem.setProblem(); + + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Initial Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + + // ******************************************************************* + // *********************Perform initial solve************************* + // ******************************************************************* + + RCP< Belos::SolverManager > initSolver + = rcp( new Belos::PseudoBlockGmresSolMgr( rcp(&initProblem,false), rcp(&belosList,false) ) ); + + // Perform solve + Belos::ReturnType ret = initSolver->solve(); + + // Compute actual residuals + bool badRes = false; + std::vector actualResids( initNumRHS ); + std::vector rhsNorm( initNumRHS ); + MV initR( Map, initNumRHS ); + OPT::Apply( *A, *initX, initR ); + MVT::MvAddMv( -1.0, initR, 1.0, *initB, initR ); + MVT::MvNorm( initR, actualResids ); + MVT::MvNorm( *initB, rhsNorm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; + } + } + + if (ret != Belos::Converged || badRes==true) { + if (procVerbose) + std::cout << std::endl << "ERROR: Initial solve did not converge to solution!" << std::endl; + return -1; + } + + // ***************Construct augmented linear system**************** + + RCP augX = rcp( new MV(Map, initNumRHS+augNumRHS) ); + RCP augB = rcp( new MV(Map, initNumRHS+augNumRHS) ); + if (augNumRHS) { + MVT::MvRandom( *augX ); + OPT::Apply( *A, *augX, *augB ); + augX->putScalar( 0.0 ); + } + + // Copy previous linear system into + RCP tmpX = rcp( new MV( *augX ) ); + RCP tmpB = rcp( new MV( *augB ) ); + tmpX->scale( 1.0, *augX ); + tmpB->scale( 1.0, *augB ); + + Belos::LinearProblem augProblem( A, augX, augB ); + augProblem.setLabel("Belos Aug"); + + set = augProblem.setProblem(); + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Augmented Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + + // ******************************************************************* + // *******************Perform augmented solve************************* + // ******************************************************************* + + belosList.set( "Block Size", augBlockSize ); // Blocksize to be used by iterative solver + belosList.set( "Convergence Tolerance", aug_tol ); // Relative convergence tolerance requested + belosList.set( "Deflation Quorum", augBlockSize ); // Number of converged linear systems before deflation + belosList.set( "Timer Label", "Belos Aug" ); // Label used by the timers in this solver + belosList.set( "Implicit Residual Scaling", "Norm of RHS" ); // Implicit residual scaling for convergence + belosList.set( "Explicit Residual Scaling", "Norm of RHS" ); // Explicit residual scaling for convergence + RCP< Belos::SolverManager > augSolver + = rcp( new Belos::PseudoBlockGmresSolMgr( rcp(&augProblem,false), rcp(&belosList,false) ) ); + + // Perform solve + ret = augSolver->solve(); + + if (ret != Belos::Converged) { + if (procVerbose) + std::cout << std::endl << "ERROR: Augmented solver did not converge to solution!" << std::endl; + return -1; + } + + // **********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 initial right-hand sides: " << initNumRHS << std::endl; + std::cout << "Number of augmented right-hand sides: " << augNumRHS << std::endl; + std::cout << "Number of restarts allowed: " << maxRestarts << std::endl; + std::cout << "Length of block Arnoldi factorization: " << length < tol ) badRes = true; + } + } + if (ret!=Belos::Converged || badRes==true) { + success = false; + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + } else { + success = true; + if (procVerbose) + std::cout << "End Result: TEST PASSED" << std::endl; + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} // run + +int main(int argc, char *argv[]) { + // run with different ST + return run(argc,argv); + // run(argc,argv); // FAILS +} diff --git a/packages/belos/tpetra/test/BlockGmres/test_resolve_gmres_hb.cpp b/packages/belos/tpetra/test/BlockGmres/test_resolve_gmres_hb.cpp new file mode 100644 index 000000000000..61874c175645 --- /dev/null +++ b/packages/belos/tpetra/test/BlockGmres/test_resolve_gmres_hb.cpp @@ -0,0 +1,432 @@ +//@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 from the problem is being used instead of multiple +// random right-hand-sides. The initial guesses are all set to zero. +// +// NOTE: No preconditioner is used in this case. +// + +// Tpetra +#include +#include +#include + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosBlockGmresSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" +#include "BelosPseudoBlockGmresSolMgr.hpp" + +template +int run(int argc, char *argv[]) { + // + Teuchos::GlobalMPISession session(&argc, &argv, nullptr); + // + using ST = typename Tpetra::Vector::scalar_type; + using LO = typename Tpetra::Vector<>::local_ordinal_type; + using GO = typename Tpetra::Vector<>::global_ordinal_type; + using NT = typename Tpetra::Vector<>::node_type; + + using OP = typename Tpetra::Operator; + using MV = typename Tpetra::MultiVector; + using SCT = typename Teuchos::ScalarTraits; + using MT = typename SCT::magnitudeType; + + using tmap_t = Tpetra::Map; + using tcrsmatrix_t = Tpetra::CrsMatrix; + + using MVT = typename Belos::MultiVecTraits; + using OPT = typename Belos::OperatorTraits; + + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::ParameterList; + + bool success = false; + bool verbose = false; + + try { + + const auto comm = Tpetra::getDefaultComm(); + const int myPID = comm->getRank(); + + bool badRes = false; + bool procVerbose = false; + bool pseudo = false; // use pseudo block Gmres to solve this linear system. + int frequency = -1; + int blockSize = 1; + int numrhs = 1; + int maxRestarts = 15; // number of restarts allowed + int maxIters = -1; // maximum number of iterations allowed per linear system + std::string filename("orsirr1.hb"); + std::string ortho("DGKS"); + MT tol = 1.0e-5; // relative residual tolerance + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block Gmres to solve the linear systems."); + cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); + cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); + cmdp.setOption("ortho",&ortho,"Orthogonalization routine used by Gmres solver."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by Gmres solver."); + cmdp.setOption("max-restarts",&maxRestarts,"Maximum number of restarts allowed for Gmres solver."); + cmdp.setOption("blockSize",&blockSize,"Block size used by Gmres."); + cmdp.setOption("maxIters",&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 + + procVerbose = ( 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 ); + + // ********Other information used by block solver*********** + // *****************(can be user specified)****************** + + const int NumGlobalElements = B->getGlobalLength(); + if (maxIters == -1) + maxIters = NumGlobalElements/blockSize - 1; // maximum number of iterations to run + + ParameterList belosList; + belosList.set( "Num Blocks", maxIters ); // Maximum number of blocks in Krylov factorization + belosList.set( "Block Size", blockSize ); // Blocksize to be used by iterative solver + belosList.set( "Maximum Iterations", maxIters ); // Maximum number of iterations allowed + belosList.set( "Maximum Restarts", maxRestarts ); // Maximum number of restarts allowed + belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + belosList.set( "Orthogonalization", ortho ); // Orthogonalization routine + 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 an unpreconditioned linear problem instance. + Belos::LinearProblem problem( A, X, B ); + 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 block Gmres iteration************************* + // ******************************************************************* + + RCP< Belos::SolverManager > solver; + if (pseudo) + solver = rcp( new Belos::PseudoBlockGmresSolMgr( rcp(&problem,false), rcp(&belosList,false) ) ); + else + solver = rcp( new Belos::BlockGmresSolMgr( rcp(&problem,false), rcp(&belosList,false) ) ); + + // Perform solve + Belos::ReturnType ret = solver->solve(); + + if (ret!=Belos::Converged) { + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + return -1; + } + + // Get the number of iterations for this solve. + int numIters = solver->getNumIters(); + std::cout << "Number of iterations performed for this solve: " << numIters << std::endl; + + // Compute actual residuals. + std::vector actual_resids( numrhs ); + std::vector rhs_norm( numrhs ); + MV resid(Map, numrhs); + OPT::Apply( *A, *X, resid ); + MVT::MvAddMv( -1.0, resid, 1.0, *B, resid ); + MVT::MvNorm( resid, actual_resids ); + MVT::MvNorm( *B, rhs_norm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"<putScalar( 0.0 ); + solver->reset( Belos::Problem ); + + // Perform solve (again) + ret = solver->solve(); + + // Get the number of iterations for this solve. + numIters = solver->getNumIters(); + std::cout << "Number of iterations performed for this solve (manager reset): " << numIters << std::endl; + + if (ret!=Belos::Converged) { + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + return -1; + } + + // Compute actual residuals. + std::vector actual_resids2( numrhs ); + MV resid2(Map, numrhs); + OPT::Apply( *A, *X, resid2 ); + MVT::MvAddMv( -1.0, resid2, 1.0, *B, resid2 ); + MVT::MvNorm( resid2, actual_resids ); + MVT::MvAddMv( -1.0, resid2, 1.0, resid, resid2 ); + MVT::MvNorm( resid2, actual_resids2 ); + MVT::MvNorm( *B, rhs_norm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (manager reset) ----------"< SCT::prec() ) { + badRes = true; + std::cout << "Resolve residual vector is too different from first solve residual vector: " << actual_resids2[i] << std::endl; + } + } + } + + // ----------------------------------------------------------------- + // Resolve the first problem by resetting the linear problem. + // ----------------------------------------------------------------- + + X->putScalar( 0.0 ); + set = problem.setProblem(); + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + solver->setProblem( rcp(&problem,false) ); + + // Perform solve (again) + ret = solver->solve(); + + // Get the number of iterations for this solve. + numIters = solver->getNumIters(); + std::cout << "Number of iterations performed for this solve (manager setProblem()): " << numIters << std::endl; + + if (ret!=Belos::Converged) { + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + return -1; + } + + // Compute actual residuals. + OPT::Apply( *A, *X, resid2 ); + MVT::MvAddMv( -1.0, resid2, 1.0, *B, resid2 ); + MVT::MvNorm( resid2, actual_resids ); + MVT::MvAddMv( -1.0, resid2, 1.0, resid, resid2 ); + MVT::MvNorm( resid2, actual_resids2 ); + MVT::MvNorm( *B, rhs_norm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (manager setProblem()) ----------"< SCT::prec() ) { + badRes = true; + std::cout << "Resolve residual vector is too different from first solve residual vector: " << actual_resids2[i] << std::endl; + } + } + } + + // ---------------------------------------------------------------------------------- + // Resolve the first problem by resetting the solver manager and changing the labels. + // ---------------------------------------------------------------------------------- + + ParameterList belosList2; + belosList2.set( "Timer Label", "Belos Resolve w/ New Label" ); // Change timer labels. + solver->setParameters( rcp( &belosList2, false ) ); + + problem.setLabel( "Belos Resolve w/ New Label" ); + X->putScalar( 0.0 ); + solver->reset( Belos::Problem ); + + // Perform solve (again) + ret = solver->solve(); + + // Get the number of iterations for this solve. + numIters = solver->getNumIters(); + std::cout << "Number of iterations performed for this solve (label reset): " << numIters << std::endl; + + if (ret!=Belos::Converged) { + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + return -1; + } + + // Compute actual residuals. + OPT::Apply( *A, *X, resid2 ); + MVT::MvAddMv( -1.0, resid2, 1.0, *B, resid2 ); + MVT::MvNorm( resid2, actual_resids ); + MVT::MvAddMv( -1.0, resid2, 1.0, resid, resid2 ); + MVT::MvNorm( resid2, actual_resids2 ); + MVT::MvNorm( *B, rhs_norm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (label reset) ----------"< SCT::prec() ) { + badRes = true; + std::cout << "Resolve residual vector is too different from first solve residual vector: " << actual_resids2[i] << std::endl; + } + } + } + + // ------------------------------------------------------------- + // Construct a second unpreconditioned linear problem instance. + // ------------------------------------------------------------- + + RCP X2 = MVT::Clone(*X, numrhs); + MVT::MvInit( *X2, 0.0 ); + Belos::LinearProblem problem2( A, X2, B ); + problem2.setLabel("Belos Resolve"); + set = problem2.setProblem(); + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + return -1; + } + + // ******************************************************************* + // *************Start the block Gmres iteration************************* + // ******************************************************************* + + // Create the solver without either the problem or parameter list. + if (pseudo) + solver = rcp( new Belos::PseudoBlockGmresSolMgr() ); + else + solver = rcp( new Belos::BlockGmresSolMgr() ); + + // Set the problem after the solver construction. + solver->setProblem( rcp( &problem2, false ) ); + + // Get the valid list of parameters from the solver and print it. + RCP validList = solver->getValidParameters(); + if (procVerbose) { + if (pseudo) + std::cout << std::endl << "Valid parameters from the pseudo-block Gmres solver manager:" << std::endl; + else + std::cout << std::endl << "Valid parameters from the block Gmres solver manager:" << std::endl; + + std::cout << *validList << std::endl; + } + + // Set the parameter list after the solver construction. + belosList.set( "Timer Label", "Belos Resolve" ); // Set timer label to discern between the two solvers. + solver->setParameters( rcp( &belosList, false ) ); + + // Perform solve + ret = solver->solve(); + + // Get the number of iterations for this solve. + numIters = solver->getNumIters(); + std::cout << "Number of iterations performed for this solve (new solver): " << numIters << std::endl; + + // Compute actual residuals. + OPT::Apply( *A, *X2, resid2 ); + MVT::MvAddMv( -1.0, resid2, 1.0, *B, resid2 ); + MVT::MvNorm( resid2, actual_resids ); + MVT::MvAddMv( -1.0, resid2, 1.0, resid, resid2 ); + MVT::MvNorm( resid2, actual_resids2 ); + MVT::MvNorm( *B, rhs_norm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (new solver) ----------"< SCT::prec() ) { + badRes = true; + std::cout << "Resolve residual vector is too different from first solve residual vector: " << actual_resids2[i] << std::endl; + } + } + } + + // **********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 << "Block size used by solver: " << blockSize << std::endl; + std::cout << "Max number of Gmres iterations: " << maxIters << std::endl; + std::cout << "Relative residual tolerance: " << tol << std::endl; + std::cout << std::endl; + } + + 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 ); +} // run + +int main(int argc, char *argv[]) { + // run with different ST + return run(argc,argv); + // run(argc,argv); // FAILS +} diff --git a/packages/belos/tpetra/test/CMakeLists.txt b/packages/belos/tpetra/test/CMakeLists.txt index 53bd917be839..b32282e51c3f 100644 --- a/packages/belos/tpetra/test/CMakeLists.txt +++ b/packages/belos/tpetra/test/CMakeLists.txt @@ -7,6 +7,7 @@ ADD_SUBDIRECTORY(TFQMR) ADD_SUBDIRECTORY(MINRES) ADD_SUBDIRECTORY(FixedPoint) ADD_SUBDIRECTORY(LinearSolverFactory) +ADD_SUBDIRECTORY(MultiMatrixSolve) ADD_SUBDIRECTORY(MultipleSolves) ADD_SUBDIRECTORY(MVOPTester) ADD_SUBDIRECTORY(Native) diff --git a/packages/belos/tpetra/test/MultiMatrixSolve/BelosLinearMultiShiftProblem.hpp b/packages/belos/tpetra/test/MultiMatrixSolve/BelosLinearMultiShiftProblem.hpp new file mode 100644 index 000000000000..3a8e6c1a3475 --- /dev/null +++ b/packages/belos/tpetra/test/MultiMatrixSolve/BelosLinearMultiShiftProblem.hpp @@ -0,0 +1,358 @@ + +//@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 + +#ifndef BELOS_LINEAR_MULTIPROBLEM_HPP +#define BELOS_LINEAR_MULTIPROBLEM_HPP + +/// \file BelosLinearMultiShiftProblem.hpp +/// \brief Class which describes the linear problem to be solved by +/// the iterative solver. +#include "BelosLinearProblem.hpp" +#include "BelosMultiVecTraits.hpp" +#include "BelosOperatorTraits.hpp" +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_TimeMonitor.hpp" + +namespace Belos { + + //@} + + /// \class LinearMultiShiftProblem + /// \brief A linear system to solve multiple matrix and right-hand sides. + /// + /// This class encapsulates the general information needed for + /// solving a set of linear systems of equations using an iterative method. + /// A_i*x_i = b_i for i=1,...,N, where A_i is A+i*I + /// + /// \tparam ScalarType The type of the entries in the matrix and + /// vectors. + /// \tparam MV The (multi)vector type. + /// \tparam OP The operator type. (Operators are functions that + /// take a multivector as input and compute a multivector as + /// output.) + template + class LinearMultiShiftProblem : public LinearProblem { + public: + + //! @name Constructors/Destructor + //@{ + + /// \brief Default constructor. + /// + /// Creates an empty Belos::LinearMultiShiftProblem instance. The operator + /// A, left-hand-side X and right-hand-side B must be set using + /// the \c setOperator(), \c setLHS(), and \c setRHS() methods + /// respectively. + LinearMultiShiftProblem (void); + + /// \brief Unpreconditioned linear system constructor. + /// + /// Creates an unpreconditioned LinearMultiShiftProblem instance with the + /// operator (\c A), initial guess (\c X), and right hand side (\c + /// B). Preconditioners can be set using the \c setLeftPrec() and + /// \c setRightPrec() methods. + LinearMultiShiftProblem (const Teuchos::RCP &A, + const Teuchos::RCP &X, + const Teuchos::RCP &B); + + //! Destructor (declared virtual for memory safety of derived classes). + virtual ~LinearMultiShiftProblem (void) {} + + //@} + + //! @name Set / Reset method + //@{ + + /// \brief Set up the linear problem manager. + /// + /// Call this method if you want to solve the linear system with a + /// different left- or right-hand side, or if you want to prepare + /// the linear problem to solve the linear system that was already + /// passed in. (In the latter case, call this method with the + /// default arguments.) The internal flags will be set as if the + /// linear system manager was just initialized, and the initial + /// residual will be computed. + /// + /// Many of Belos' solvers require that this method has been + /// called on the linear problem, before they can solve it. + /// + /// \param newX [in/out] If you want to solve the linear system + /// with a different left-hand side, pass it in here. + /// Otherwise, set this to null (the default value). + /// + /// \param newB [in] If you want to solve the linear system with a + /// different right-hand side, pass it in here. Otherwise, set + /// this to null (the default value). + /// + /// \return Whether the linear problem was successfully set up. + /// Successful setup requires at least that the matrix operator + /// A, the left-hand side X, and the right-hand side B all be + /// non-null. + bool + setProblem (const Teuchos::RCP &newX = Teuchos::null, + const Teuchos::RCP &newB = Teuchos::null); + + + void + setShift( bool addShift ) + { + addShift_ = addShift; + } + + //@} + + //! @name Apply / Compute methods + //@{ + + /// \brief Apply ONLY the operator to \c x, returning \c y. + /// + /// This method only applies the linear problem's operator, + /// without any preconditioners that may have been defined. + /// Flexible variants of Krylov methods will use this method. If + /// no operator has been defined, this method just copies x into + /// y. + void applyOp( const MV& x, MV& y ) const; + + //@} + + private: + bool addShift_; + + using LinearProblem::A_; + using LinearProblem::X_; + using LinearProblem::curX_; + using LinearProblem::B_; + using LinearProblem::curB_; + using LinearProblem::R0_; + using LinearProblem::PR0_; + using LinearProblem::R0_user_; + using LinearProblem::PR0_user_; + using LinearProblem::LP_; + using LinearProblem::RP_; + using LinearProblem::timerOp_; + using LinearProblem::timerPrec_; + using LinearProblem::blocksize_; + using LinearProblem::num2Solve_; + using LinearProblem::rhsIndex_; + using LinearProblem::lsNum_; + using LinearProblem::isSet_; + using LinearProblem::solutionUpdated_; + using LinearProblem::label_; + + typedef MultiVecTraits MVT; + typedef OperatorTraits OPT; + }; + + //-------------------------------------------- + // Constructor Implementations + //-------------------------------------------- + + template + LinearMultiShiftProblem::LinearMultiShiftProblem(void) + : LinearProblem(), + addShift_(true) + {} + + template + LinearMultiShiftProblem::LinearMultiShiftProblem(const Teuchos::RCP &A, + const Teuchos::RCP &X, + const Teuchos::RCP &B + ) + : LinearProblem(A, X, B), + addShift_(true) + { + int numVecs = MVT::GetNumberVecs( *B ); + rhsIndex_.resize(numVecs); + for (int i=0; i + bool + LinearMultiShiftProblem:: + setProblem (const Teuchos::RCP &newX, + const Teuchos::RCP &newB) + { + // Create timers if the haven't been created yet. + if (timerOp_ == Teuchos::null) { + std::string opLabel = label_ + ": Operation Op*x"; +#ifdef BELOS_TEUCHOS_TIME_MONITOR + timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel ); +#endif + } + if (timerPrec_ == Teuchos::null) { + std::string precLabel = label_ + ": Operation Prec*x"; +#ifdef BELOS_TEUCHOS_TIME_MONITOR + timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel ); +#endif + } + + // Set the linear system using the arguments newX and newB + if (newX != Teuchos::null) + X_ = newX; + if (newB != Teuchos::null) + B_ = newB; + + // Invalidate the current linear system indices and multivectors. + curX_ = Teuchos::null; + curB_ = Teuchos::null; + + // Temporarily set the rhsIndex_ to all the vectors of the linear system. + // This is necessary for the applyOp to perform correctly in computing the + // initial residual. It will be reset at the end of this method. + int numVecs = MVT::GetNumberVecs( *B_ ); + rhsIndex_.resize(numVecs); + for (int i=0; icomputeCurrResVec( &*R0_, &*X_, &*B_ ); + + if (LP_!=Teuchos::null) { + if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) { + PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) ); + } + this->applyLeftPrec( *R0_, *PR0_ ); + } + else { + PR0_ = R0_; + } + } + else { // User specified the residuals + // If the user did not specify the right sized residual, create one and set R0_user_ to null. + if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) { + Teuchos::RCP helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) ); + this->computeCurrResVec( &*helper, &*X_, &*B_ ); + R0_user_ = Teuchos::null; + R0_ = helper; + } + + if (LP_!=Teuchos::null) { + // If the user provided preconditioned residual is the wrong size or pointing at + // the wrong object, create one and set the PR0_user_ to null. + if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (PR0_user_==R0_user_) + || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) { + Teuchos::RCP helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) ); + + // Get the initial residual from getInitResVec because R0_user_ may be null from above. + this->applyLeftPrec( *(this->getInitResVec()), *helper ); + PR0_user_ = Teuchos::null; + PR0_ = helper; + } + } + else { + // The preconditioned initial residual vector is the residual vector. + // NOTE: R0_user_ could be null if incompatible. + if (R0_user_!=Teuchos::null) + { + PR0_user_ = R0_user_; + } + else + { + PR0_user_ = Teuchos::null; + PR0_ = R0_; + } + } + } + + // Resize the rhsIndex_ back to zero length to correspond with null pointers for curX_ and curB_. + rhsIndex_.resize(0); + + // The problem has been set and is ready for use. + isSet_ = true; + + // Return isSet. + return isSet_; + } + + template + void LinearMultiShiftProblem::applyOp( const MV& x, MV& y ) const { + + if (A_.get()) + { + { +#ifdef BELOS_TEUCHOS_TIME_MONITOR + Teuchos::TimeMonitor OpTimer(*timerOp_); +#endif + OPT::Apply( *A_,x, y); + } + // Now add in shift based on rhsIndex_ + if (addShift_) + { + std::vector newIdx( 1 ); + for (unsigned int i=0; i tmp_x = MVT::CloneView( x, newIdx ); + Teuchos::RCP tmp_y = MVT::CloneViewNonConst( y, newIdx ); + MVT::MvAddMv( rhsIndex_[i]*Teuchos::ScalarTraits::one(), *tmp_x, + Teuchos::ScalarTraits::one(), *tmp_y, *tmp_y ); + } + } + } + else { + MVT::MvAddMv( Teuchos::ScalarTraits::one(), x, + Teuchos::ScalarTraits::zero(), x, y ); + } + } + +} // end Belos namespace + +#endif /* BELOS_LINEAR_MULTIPROBLEM_HPP */ + + diff --git a/packages/belos/tpetra/test/MultiMatrixSolve/CMakeLists.txt b/packages/belos/tpetra/test/MultiMatrixSolve/CMakeLists.txt new file mode 100644 index 000000000000..acb6320245a3 --- /dev/null +++ b/packages/belos/tpetra/test/MultiMatrixSolve/CMakeLists.txt @@ -0,0 +1,19 @@ + +ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) +IF (${PACKAGE_NAME}_ENABLE_Triutils) + + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_pseudo_gmres_multi_hb + SOURCES test_pseudo_gmres_multi_hb.cpp BelosLinearMultiShiftProblem.hpp + COMM serial mpi + ARGS "--verbose --filename=orsirr1.hb" + STANDARD_PASS_OUTPUT + ) + + TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyTestGmresMultiFiles + SOURCE_DIR ${Belos_SOURCE_DIR}/tpetra/example/BlockGmres + SOURCE_FILES orsirr1.hb + EXEDEPS Tpetra_pseudo_gmres_multi_hb + ) + +ENDIF(${PACKAGE_NAME}_ENABLE_Triutils) diff --git a/packages/belos/tpetra/test/MultiMatrixSolve/test_pseudo_gmres_multi_hb.cpp b/packages/belos/tpetra/test/MultiMatrixSolve/test_pseudo_gmres_multi_hb.cpp new file mode 100644 index 000000000000..74c2e243190c --- /dev/null +++ b/packages/belos/tpetra/test/MultiMatrixSolve/test_pseudo_gmres_multi_hb.cpp @@ -0,0 +1,328 @@ +//@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. +// +// This test is for testing the deflation in the pseudo-block Gmres solver. +// One set of linear systems is solved and then augmented with additional +// linear systems and resolved. The already solved linear systems should be +// deflated immediately, leaving only the augmented systems to be solved. +// +// + +// Tpetra +#include +#include +#include +#include + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosTpetraTestFramework.hpp" +#include "BelosPseudoBlockGmresSolMgr.hpp" +#include "BelosLinearMultiShiftProblem.hpp" + +template +int run(int argc, char *argv[]) { + + using ST = typename Tpetra::CrsMatrix::scalar_type; + using LO = typename Tpetra::CrsMatrix<>::local_ordinal_type; + using GO = typename Tpetra::CrsMatrix<>::global_ordinal_type; + using NT = typename Tpetra::CrsMatrix<>::node_type; + + using OP = typename Tpetra::Operator; + using MV = typename Tpetra::MultiVector; + using MT = typename Teuchos::ScalarTraits::magnitudeType; + using MVT = typename Belos::MultiVecTraits; + + using tmap_t = Tpetra::Map; + using tcrsmatrix_t = Tpetra::CrsMatrix; + + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::ParameterList; + + Teuchos::GlobalMPISession session(&argc, &argv, nullptr); + + bool verbose = false; + bool success = true; + + try { + + const auto comm = Tpetra::getDefaultComm(); + const int myPID = comm->getRank(); + + bool procVerbose = false; + bool useShift = true; + int frequency = -1; // how often residuals are printed by solver + int initNumRHS = 3; // how many right-hand sides get solved first + int augNumRHS = 2; // how many right-hand sides are augmented to the first group + int maxRestarts = 15; // number of restarts allowed + int length = 100; + int initBlockSize = 2;// blocksize used for the initial pseudo-block GMRES solve + int augBlockSize = 1; // blocksize used for the augmented pseudo-block GMRES solve + int maxIters = -1; // maximum iterations allowed + std::string filename("orsirr1.hb"); + MT tol = 1.0e-5; // relative residual tolerance + MT aug_tol = 1.0e-5; // relative residual tolerance for augmented system + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption("verbose","quiet",&verbose,"Print messages and results."); + cmdp.setOption("use-shift", "no-shift", &useShift, "Whether shift should be used."); + cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)."); + cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix."); + cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver."); + cmdp.setOption("aug-tol",&aug_tol,"Relative residual tolerance used by GMRES solver for augmented systems."); + cmdp.setOption("init-num-rhs",&initNumRHS,"Number of right-hand sides to be initially solved for."); + cmdp.setOption("aug-num-rhs",&augNumRHS,"Number of right-hand sides augmenting the initial solve."); + cmdp.setOption("max-restarts",&maxRestarts,"Maximum number of restarts allowed for GMRES solver."); + cmdp.setOption("block-size",&initBlockSize,"Block size used by GMRES for the initial solve."); + cmdp.setOption("aug-block-size",&augBlockSize,"Block size used by GMRES for the augmented solve."); + cmdp.setOption("max-iters",&maxIters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size)."); + cmdp.setOption("subspace-size",&length,"Dimension of Krylov subspace used by GMRES."); + if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (!verbose) { + frequency = -1; // reset frequency if test is not verbose + } + procVerbose = verbose && (myPID==0); /* Only print on zero processor */ + + // Get the problem + Belos::Tpetra::HarwellBoeingReader reader( comm ); + RCP A = reader.readFromFile( filename ); + RCP Map = A->getDomainMap(); + + // ********Other information used by block solver*********** + // *****************(can be user specified)****************** + + const int numGlobalElements = Map->getGlobalNumElements(); + if (maxIters == -1) + maxIters = numGlobalElements - 1; // maximum number of iterations to run + + ParameterList belosList; + belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization + belosList.set( "Block Size", initBlockSize ); // Blocksize to be used by iterative solver + belosList.set( "Maximum Iterations", maxIters ); // Maximum number of iterations allowed + belosList.set( "Maximum Restarts", maxRestarts ); // Maximum number of restarts allowed + belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested + belosList.set( "Deflation Quorum", initBlockSize ); // Number of converged linear systems before deflation + belosList.set( "Timer Label", "Belos Init" ); // Label used by the timers in this solver + 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 initX = rcp( new MV(Map, initNumRHS) ); + RCP initB = rcp( new MV(Map, initNumRHS) ); + + Belos::LinearMultiShiftProblem initProblem( A, initX, initB ); + initProblem.setLabel("Belos Init"); + initProblem.setShift( useShift ); + bool set = initProblem.setProblem(); + + MVT::MvRandom( *initX ); + initProblem.applyOp( *initX, *initB ); + initX->putScalar( 0.0 ); + set = initProblem.setProblem( initX, initB ); + + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Initial Belos::LinearMultiShiftProblem failed to set up correctly!" << std::endl; + return -1; + } + + // ******************************************************************* + // *********************Perform initial solve************************* + // ******************************************************************* + + RCP< Belos::SolverManager > initSolver + = rcp( new Belos::PseudoBlockGmresSolMgr( rcp(&initProblem,false), rcp(&belosList,false) ) ); + + // Perform solve + Belos::ReturnType ret = initSolver->solve(); + + // Compute actual residuals. + std::vector initIndex( initNumRHS ); + for (int i=0; i actualResids( initNumRHS ); + std::vector rhsNorm( initNumRHS ); + MV initR( Map, initNumRHS ); + initProblem.applyOp( *initX, initR ); + MVT::MvAddMv( -1.0, initR, 1.0, *initB, initR ); + MVT::MvNorm( initR, actualResids ); + MVT::MvNorm( *initB, rhsNorm ); + if (procVerbose) { + std::cout<< "---------- Actual Residuals (normalized) ----------"< tol) badRes = true; + } + } + + if (ret != Belos::Converged || badRes==true) { + if (procVerbose) + std::cout << std::endl << "ERROR: Initial solve did not converge to solution!" << std::endl; + return -1; + } + + // ***************Construct augmented linear system**************** + + RCP augX = rcp( new MV(Map, initNumRHS+augNumRHS) ); + RCP augB = rcp( new MV(Map, initNumRHS+augNumRHS) ); + Belos::LinearMultiShiftProblem augProblem( A, augX, augB ); + augProblem.setLabel("Belos Aug"); + augProblem.setShift( useShift ); + set = augProblem.setProblem(); + + if (augNumRHS) { + MVT::MvRandom( *augX ); + augProblem.applyOp( *augX, *augB ); + augX->putScalar( 0.0 ); + } + + // Copy previous linear system into + RCP tmpX = rcp( new MV( *augX ) ); + RCP tmpB = rcp( new MV( *augB ) ); + tmpX->scale( 1.0, *augX ); + tmpB->scale( 1.0, *augB ); + + set = augProblem.setProblem( augX, augB ); + if (set == false) { + if (procVerbose) + std::cout << std::endl << "ERROR: Augmented Belos::LinearMultiShiftProblem failed to set up correctly!" << std::endl; + return -1; + } + + // ******************************************************************* + // *******************Perform augmented solve************************* + // ******************************************************************* + + belosList.set( "Block Size", augBlockSize ); // Blocksize to be used by iterative solver + belosList.set( "Convergence Tolerance", aug_tol ); // Relative convergence tolerance requested + belosList.set( "Deflation Quorum", augBlockSize ); // Number of converged linear systems before deflation + belosList.set( "Timer Label", "Belos Aug" ); // Label used by the timers in this solver + belosList.set( "Implicit Residual Scaling", "Norm of RHS" ); // Implicit residual scaling for convergence + belosList.set( "Explicit Residual Scaling", "Norm of RHS" ); // Explicit residual scaling for convergence + RCP< Belos::SolverManager > augSolver + = rcp( new Belos::PseudoBlockGmresSolMgr( rcp(&augProblem,false), rcp(&belosList,false) ) ); + + // Perform solve + ret = augSolver->solve(); + + if (ret != Belos::Converged) { + if (procVerbose) + std::cout << std::endl << "ERROR: Augmented solver did not converge to solution!" << std::endl; + return -1; + } + + // **********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 initial right-hand sides: " << initNumRHS << std::endl; + std::cout << "Number of augmented right-hand sides: " << augNumRHS << std::endl; + std::cout << "Number of restarts allowed: " << maxRestarts << std::endl; + std::cout << "Length of Krylov subspace: " << length < totalIndex( totalNumRHS ); + for (int i=0; i tol ) badRes = true; + } + } + if (ret!=Belos::Converged || badRes==true) { + success = false; + if (procVerbose) + std::cout << "End Result: TEST FAILED" << std::endl; + } else { + success = true; + if (procVerbose) + std::cout << "End Result: TEST PASSED" << std::endl; + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} // run + +int main(int argc, char *argv[]) { + // run with different ST + return run(argc,argv); + // run(argc,argv); // FAILS +}