forked from trilinos/Trilinos
-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
4 changed files
with
295 additions
and
0 deletions.
There are no files selected for viewing
246 changes: 246 additions & 0 deletions
246
packages/belos/tpetra/example/BiCGStab/BiCGStabTpetraExFile.cpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 ([email protected]) | ||
// | ||
// ************************************************************************ | ||
//@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 <Tpetra_Map.hpp> | ||
#include <Tpetra_Core.hpp> | ||
#include <Tpetra_Vector.hpp> | ||
#include <Tpetra_CrsMatrix.hpp> | ||
|
||
// Belos | ||
#include "BelosConfigDefs.hpp" | ||
#include "BelosLinearProblem.hpp" | ||
#include "BelosBiCGStabSolMgr.hpp" | ||
#include "BelosTpetraTestFramework.hpp" | ||
|
||
// Teuchos | ||
#include <Teuchos_RCP.hpp> | ||
#include <Teuchos_Comm.hpp> | ||
#include <Teuchos_CommHelpers.hpp> | ||
#include <Teuchos_DefaultComm.hpp> | ||
#include "Teuchos_ParameterList.hpp" | ||
#include "Teuchos_StandardCatchMacros.hpp" | ||
#include "Teuchos_CommandLineProcessor.hpp" | ||
|
||
|
||
template <typename ScalarType> | ||
int run(int argc, char *argv[]) { | ||
|
||
using ST = typename Tpetra::MultiVector<ScalarType>::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<ST>; | ||
using MT = typename SCT::magnitudeType; | ||
using MV = typename Tpetra::MultiVector<ST,LO,GO,NT>; | ||
using OP = typename Tpetra::Operator<ST,LO,GO,NT>; | ||
|
||
using tmap_t = Tpetra::Map<LO,GO,NT>; | ||
using tvector_t = Tpetra::Vector<ST,LO,GO,NT>; | ||
using tcrsmatrix_t = Tpetra::CrsMatrix<ST,LO,GO,NT>; | ||
|
||
using MVT = typename Belos::MultiVecTraits<ST,MV>; | ||
using OPT = typename Belos::OperatorTraits<ST,MV,OP>; | ||
|
||
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<tcrsmatrix_t> reader( Comm ); | ||
RCP<tcrsmatrix_t> A = reader.readFromFile( filename ); | ||
RCP<const tmap_t> Map = A->getDomainMap(); | ||
|
||
// Create initial vectors | ||
RCP<MV> 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<ParameterList> 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<Belos::LinearProblem<ST,MV,OP> > problem | ||
= rcp (new Belos::LinearProblem<ST,MV,OP> (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<Belos::SolverManager<ST,MV,OP> > solver | ||
= rcp (new Belos::BiCGStabSolMgr<ST,MV,OP> (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<ST> actual_resids (numrhs); | ||
std::vector<ST> 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<double>(argc,argv); | ||
// run<float>(argc,argv); // FAILS | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,3 +1,4 @@ | ||
|
||
ADD_SUBDIRECTORY(BiCGStab) | ||
ADD_SUBDIRECTORY(WrapTpetraSolver) | ||
ADD_SUBDIRECTORY(Orthog) |