diff --git a/packages/belos/tpetra/example/BiCGStab/BiCGStabTpetraExFile.cpp b/packages/belos/tpetra/example/BiCGStab/BiCGStabTpetraExFile.cpp new file mode 100644 index 000000000000..f1580f19d15e --- /dev/null +++ b/packages/belos/tpetra/example/BiCGStab/BiCGStabTpetraExFile.cpp @@ -0,0 +1,246 @@ +//@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 +// +// Update of belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp +// +// This driver reads a problem from a Harwell-Boeing (HB) file. +// The right-hand-side corresponds to a randomly generated solution. +// The initial guesses are all set to zero. +// +// NOTE: No preconditioner is used in this example. +// + +// Tpetra +#include +#include +#include +#include + +// Belos +#include "BelosConfigDefs.hpp" +#include "BelosLinearProblem.hpp" +#include "BelosBiCGStabSolMgr.hpp" +#include "BelosTpetraTestFramework.hpp" + +// Teuchos +#include +#include +#include +#include +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_StandardCatchMacros.hpp" +#include "Teuchos_CommandLineProcessor.hpp" + + +template +int run(int argc, char *argv[]) { + + 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 SCT = typename Teuchos::ScalarTraits; + using MT = typename SCT::magnitudeType; + using MV = typename Tpetra::MultiVector; + using OP = typename Tpetra::Operator; + + using tmap_t = Tpetra::Map; + using tvector_t = Tpetra::Vector; + 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 mpiSession (&argc, &argv, &std::cout); + const auto Comm = Tpetra::getDefaultComm(); + const int MyPID = Comm->getRank(); + + bool verbose = false; + bool success = true; + + // This "try" relates to TEUCHOS_STANDARD_CATCH_STATEMENTS near the + // bottom of main(). That macro has the corresponding "catch". + try { + bool proc_verbose = false; + int frequency = -1; // frequency of status test output. + int numrhs = 1; // number of right-hand sides to solve for + int maxiters = -1; // maximum number of iterations allowed per linear system + std::string filename ("cage4.hb"); + MT tol = 1.0e-5; // relative residual tolerance + + Teuchos::CommandLineProcessor cmdp(false,true); + cmdp.setOption ("verbose", "quiet", &verbose, "Whether to print messages " + "and results."); + cmdp.setOption ("frequency", &frequency, "Frequency of solver output " + "(1 means every iteration; -1 means never)."); + cmdp.setOption ("filename", &filename, "Test matrix filename. " + "Allowed file extensions: *.hb, *.mtx, *.triU, *.triS"); + cmdp.setOption ("tol", &tol, "Relative residual tolerance for solver."); + cmdp.setOption ("num-rhs", &numrhs, "Number of right-hand sides to solve."); + cmdp.setOption ("max-iters", &maxiters, "Maximum number of iterations per " + "linear system (-1 = adapted to problem/block size)."); + if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) { + return -1; + } + if (!verbose) + frequency = -1; // reset frequency if test is not verbose + + proc_verbose = ( verbose && (MyPID==0) ); + + if (proc_verbose) { + std::cout << Belos::Belos_Version() << std::endl << std::endl; + } + // + // Get the problem + // + Belos::Tpetra::HarwellBoeingReader reader( Comm ); + RCP A = reader.readFromFile( filename ); + RCP Map = A->getDomainMap(); + + // Create initial vectors + RCP B, X; + X = rcp( new MV(Map,numrhs) ); + MVT::MvRandom( *X ); + B = rcp( new MV(Map,numrhs) ); + OPT::Apply( *A, *X, *B ); + MVT::MvInit( *X, 0.0 ); + + // + // Create parameter list for the Belos solver + // + const int NumGlobalElements = B->getGlobalLength (); + if (maxiters == -1) { + maxiters = NumGlobalElements - 1; // maximum number of iterations to run + } + RCP belosList (new ParameterList ("Belos")); + belosList->set ("Maximum Iterations", maxiters); + belosList->set ("Convergence Tolerance", tol); + if (verbose) { + belosList->set ("Verbosity", Belos::Errors + Belos::Warnings + + Belos::TimingDetails + Belos::StatusTestDetails); + if (frequency > 0) { + belosList->set ("Output Frequency", frequency); + } + } + else { + belosList->set ("Verbosity", Belos::Errors + Belos::Warnings + + Belos::FinalSummary); + } + // + // Construct a preconditioned linear problem + // + RCP > problem + = rcp (new Belos::LinearProblem (A, X, B)); + bool set = problem->setProblem (); + if (! set) { + if (proc_verbose) { + std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; + } + return -1; + } + + // Create a Belos solver. + RCP > solver + = rcp (new Belos::BiCGStabSolMgr (problem, belosList)); + + if (proc_verbose) { + std::cout << std::endl << std::endl; + std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl; + std::cout << "Number of right-hand sides: " << numrhs << std::endl; + std::cout << "Max number of CG iterations: " << maxiters << std::endl; + std::cout << "Relative residual tolerance: " << tol << std::endl; + std::cout << std::endl; + } + + // Ask Belos to solve the linear system. + Belos::ReturnType ret = solver->solve(); + + // + // Compute actual residuals. + // + bool badRes = false; + std::vector actual_resids (numrhs); + std::vector rhs_norm (numrhs); + MV resid (Map, numrhs); + OPT::Apply (*A, *X, resid); + MVT::MvAddMv (-1.0, resid, 1.0, *B, resid); + MVT::MvNorm (resid, actual_resids); + MVT::MvNorm (*B, rhs_norm); + if (proc_verbose) { + std::cout << "---------- Actual Residuals (normalized) ----------" << std::endl << std::endl; + for (int i = 0; i < numrhs; ++i) { + ST actRes = actual_resids[i] / rhs_norm[i]; + std::cout <<"Problem " << i << " : \t" << actRes << std::endl; + if (actRes > tol) { + badRes = true; + } + } + } + + if (ret != Belos::Converged || badRes) { + success = false; + if (proc_verbose) { + std::cout << std::endl << "ERROR: Belos did not converge!" << std::endl + << "End Result: TEST FAILED" << std::endl; + } + } else { + success = true; + if (proc_verbose) { + std::cout << std::endl << "SUCCESS: Belos converged!" << std::endl + << "End Result: TEST PASSED" << std::endl; + } + } + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); + + return success ? EXIT_SUCCESS : EXIT_FAILURE; +} + +int main(int argc, char *argv[]) { + // run with different ST + return run(argc,argv); + // run(argc,argv); // FAILS +} diff --git a/packages/belos/tpetra/example/BiCGStab/CMakeLists.txt b/packages/belos/tpetra/example/BiCGStab/CMakeLists.txt new file mode 100644 index 000000000000..b97c4fc86e86 --- /dev/null +++ b/packages/belos/tpetra/example/BiCGStab/CMakeLists.txt @@ -0,0 +1,24 @@ + +ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) + + +# Dependency on Triutils fixes Bug 6449. +IF (${PACKAGE_NAME}_ENABLE_Triutils) + + TRIBITS_ADD_EXECUTABLE_AND_TEST( + Tpetra_BiCGStab_file + SOURCES BiCGStabTpetraExFile.cpp + COMM serial mpi + ARGS "--verbose --max-iters=10" + "--verbose --max-iters=10 --tol=1e-10" + "--verbose --max-iters=10 --num-rhs=2" + STANDARD_PASS_OUTPUT + ) + + TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyExBiCGSTABFiles + SOURCE_DIR ${${PACKAGE_NAME}_SOURCE_DIR}/tpetra/example/BiCGStab + SOURCE_FILES cage4.hb + EXEDEPS BiCGStab_file + ) + +ENDIF() diff --git a/packages/belos/tpetra/example/BiCGStab/cage4.hb b/packages/belos/tpetra/example/BiCGStab/cage4.hb new file mode 100644 index 000000000000..81e5800e5534 --- /dev/null +++ b/packages/belos/tpetra/example/BiCGStab/cage4.hb @@ -0,0 +1,24 @@ +vanHeukelum/cage4; 2003; A. van Heukelum; ed: T. Davis |905 + 20 1 2 17 +rua 9 9 49 0 +(26I3) (40I2) (3E23.15) + 1 6 11 16 21 27 33 39 45 50 + 1 2 4 5 8 1 2 3 5 6 2 3 4 6 7 1 3 4 7 8 1 2 5 6 7 9 2 3 5 6 8 9 3 4 5 7 8 9 1 4 + 6 7 8 9 5 6 7 8 9 + 0.750000000000000E+00 0.750276671145870E-01 0.916389995520797E-01 + 0.375138335572935E-01 0.458194997760398E-01 0.137458499328119E+00 + 0.687569167786467E+00 0.916389995520797E-01 0.375138335572935E-01 + 0.458194997760398E-01 0.112541500671881E+00 0.666666666666667E+00 + 0.137458499328120E+00 0.458194997760398E-01 0.375138335572935E-01 + 0.112541500671881E+00 0.750276671145870E-01 0.729097498880199E+00 + 0.375138335572935E-01 0.458194997760398E-01 0.137458499328119E+00 + 0.750276671145870E-01 0.537513833557293E+00 0.750276671145870E-01 + 0.916389995520797E-01 0.833333333333333E-01 0.112541500671881E+00 + 0.916389995520797E-01 0.137458499328120E+00 0.445874834005214E+00 + 0.137458499328120E+00 0.750276671145870E-01 0.750276671145870E-01 + 0.137458499328120E+00 0.112541500671881E+00 0.470791832661453E+00 + 0.112541500671881E+00 0.916389995520797E-01 0.112541500671881E+00 + 0.916389995520797E-01 0.750276671145870E-01 0.916389995520797E-01 + 0.545819499776040E+00 0.833333333333333E-01 0.250000000000000E+00 + 0.150055334229174E+00 0.183277999104159E+00 0.250000000000000E+00 + 0.166666666666667E+00 diff --git a/packages/belos/tpetra/example/CMakeLists.txt b/packages/belos/tpetra/example/CMakeLists.txt index 85e9ddf20f5f..70dbfa9a89bc 100644 --- a/packages/belos/tpetra/example/CMakeLists.txt +++ b/packages/belos/tpetra/example/CMakeLists.txt @@ -1,3 +1,4 @@ +ADD_SUBDIRECTORY(BiCGStab) ADD_SUBDIRECTORY(WrapTpetraSolver) ADD_SUBDIRECTORY(Orthog)