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.
Merge pull request #168 from NexGenAnalytics/163-belos-add-LSQR-examples
#163: Belos: Provide `Tpetra` version of LSQR examples
- Loading branch information
Showing
4 changed files
with
694 additions
and
0 deletions.
There are no files selected for viewing
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 |
---|---|---|
|
@@ -2,3 +2,4 @@ | |
ADD_SUBDIRECTORY(WrapTpetraSolver) | ||
ADD_SUBDIRECTORY(Orthog) | ||
ADD_SUBDIRECTORY(GCRODR) | ||
ADD_SUBDIRECTORY(LSQR) |
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,31 @@ | ||
|
||
|
||
ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Triutils) | ||
ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Tpetra) | ||
|
||
IF (${PACKAGE_NAME}_ENABLE_Triutils AND ${PACKAGE_NAME}_ENABLE_Tpetra) | ||
TRIBITS_ADD_EXECUTABLE( | ||
LSQR_Tpetra_File_Ex | ||
SOURCES LSQRTpetraExFile.cpp | ||
COMM serial mpi | ||
) | ||
|
||
TRIBITS_COPY_FILES_TO_BINARY_DIR(Tpetra_CopyExampleLSQRFiles | ||
SOURCE_DIR ${Belos_SOURCE_DIR}/epetra/example/BlockGmres | ||
SOURCE_FILES orsirr1_scaled.hb | ||
EXEDEPS LSQR_Tpetra_File_Ex | ||
) | ||
|
||
# TODO: SOLVE CIRCULAR DEPENDENCY between Belos and Ifpack2 | ||
# ASSERT_DEFINED(${PACKAGE_NAME}_ENABLE_Ifpack2) | ||
# IF(${PACKAGE_NAME}_ENABLE_Ifpack2) | ||
|
||
# TRIBITS_ADD_EXECUTABLE( | ||
# PrecLSQR_Tpetra_File_Ex | ||
# SOURCES PrecLSQRTpetraExFile.cpp | ||
# COMM serial mpi | ||
# ) | ||
|
||
# ENDIF(${PACKAGE_NAME}_ENABLE_Ifpack2) | ||
|
||
ENDIF(${PACKAGE_NAME}_ENABLE_Tpetra) |
323 changes: 323 additions & 0 deletions
323
packages/belos/tpetra/example/LSQR/LSQRTpetraExFile.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,323 @@ | ||
//@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 | ||
// | ||
// This driver reads a problem from a file, which can be in Harwell-Boeing (*.hb), | ||
// Matrix Market (*.mtx), or triplet format (*.triU, *.triS). The right-hand side | ||
// from the problem, if it exists, will be used instead of multiple | ||
// random right-hand sides. The initial guesses are all set to zero. | ||
// | ||
// NOTE: No preconditioner is used in this example. | ||
// | ||
|
||
// Teuchos | ||
#include <Teuchos_CommandLineProcessor.hpp> | ||
#include <Teuchos_GlobalMPISession.hpp> | ||
#include <Teuchos_ParameterList.hpp> | ||
#include <Teuchos_StandardCatchMacros.hpp> | ||
|
||
// Tpetra | ||
#include <MatrixMarket_Tpetra.hpp> | ||
#include <Tpetra_Core.hpp> | ||
#include <Tpetra_CrsMatrix.hpp> | ||
|
||
// Belos | ||
#include "BelosConfigDefs.hpp" | ||
#include "BelosLSQRSolMgr.hpp" | ||
#include "BelosLinearProblem.hpp" | ||
#include "BelosTpetraAdapter.hpp" | ||
#include "BelosTpetraTestFramework.hpp" | ||
|
||
template <typename ScalarType> | ||
int run(int argc, char *argv[]) { | ||
using Teuchos::CommandLineProcessor; | ||
using Teuchos::GlobalMPISession; | ||
using Teuchos::ParameterList; | ||
using Teuchos::RCP; | ||
using Teuchos::rcp; | ||
using Teuchos::rcp_implicit_cast; | ||
using Teuchos::tuple; | ||
|
||
using ST = typename Tpetra::MultiVector<ScalarType>::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 MV = typename Tpetra::MultiVector<ST, LO, GO, NT>; | ||
using OP = typename Tpetra::Operator<ST, LO, GO, NT>; | ||
using MAP = typename Tpetra::Map<LO, GO, NT>; | ||
using MAT = typename Tpetra::CrsMatrix<ST, LO, GO, NT>; | ||
|
||
using MVT = typename Belos::MultiVecTraits<ST, MV>; | ||
using OPT = typename Belos::OperatorTraits<ST, MV, OP>; | ||
|
||
using MT = typename Teuchos::ScalarTraits<ST>::magnitudeType; | ||
using STM = Teuchos::ScalarTraits<MT>; | ||
using STS = Teuchos::ScalarTraits<ST>; | ||
|
||
using LinearProblem = typename Belos::LinearProblem<ST, MV, OP>; | ||
using Solver = ::Belos::LSQRSolMgr<ST, MV, OP>; | ||
|
||
Teuchos::GlobalMPISession session(&argc, &argv, NULL); | ||
RCP<const Teuchos::Comm<int>> comm = Tpetra::getDefaultComm(); | ||
|
||
bool verbose = false; | ||
bool success = true; | ||
|
||
try { | ||
bool procVerbose = false; | ||
bool debug = false; | ||
int frequency = -1; // frequency of status test output. | ||
int blockSize = 1; // blockSize | ||
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("orsirr1_scaled.hb"); | ||
std::string filenameRHS; // blank mean unset | ||
MT relResTol = STM::squareroot(STS::eps()); // relative residual tolerance | ||
// Like CG, LSQR is a short recurrence method that | ||
// does not have the "n" step convergence property in finite precision arithmetic. | ||
MT resGrowthFactor = 4.0; // In this example, warn if |resid| > resGrowthFactor * relResTol | ||
// With no preconditioner, this is only the difference between the "implicit" and the "explict | ||
// residual. | ||
|
||
MT relMatTol = 1.e-4; // relative Matrix error, default value sqrt(eps) | ||
MT maxCond = 1.e+8; // maximum condition number default value 1/eps | ||
MT damp = 0.; // regularization (or damping) parameter | ||
|
||
Teuchos::CommandLineProcessor cmdp(false, true); // e.g. ./a.out --tol=.1 --filename=foo.hb | ||
|
||
cmdp.setOption("verbose", "quiet", &verbose, "Print messages and results."); | ||
cmdp.setOption("debug", "nondebug", &debug, "Print debugging information from solver."); | ||
cmdp.setOption("frequency", &frequency, "Solvers frequency for printing residuals (#iters)."); | ||
cmdp.setOption( | ||
"filename", &filename, | ||
"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS"); | ||
cmdp.setOption("rhsFilename", &filenameRHS, | ||
"Filename for right-hand side. Acceptable file extension: *.mtx"); | ||
cmdp.setOption("lambda", &damp, "Regularization parameter"); | ||
cmdp.setOption("tol", &relResTol, "Relative residual tolerance"); | ||
cmdp.setOption("matrixTol", &relMatTol, "Relative error in Matrix"); | ||
cmdp.setOption("max-cond", &maxCond, "Maximum condition number"); | ||
cmdp.setOption("num-rhs", &numRHS, "Number of right-hand sides to be solved for."); | ||
cmdp.setOption("block-size", &blockSize, | ||
"Block size used by LSQR."); // must be one at this point | ||
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<MAT> reader(comm); | ||
RCP<MAT> A = reader.readFromFile(filename); | ||
RCP<const MAP> map = A->getDomainMap(); | ||
|
||
// Initialize vectors | ||
RCP<MV> vecB = rcp(new MV(map, numRHS)); | ||
RCP<MV> vecX = rcp(new MV(map, numRHS)); | ||
RCP<MV> B, X; | ||
|
||
// Rectangular matrices are embedded in square matrices. vecX := 0, vecB = A*randVec | ||
procVerbose = verbose && (comm->getRank() == 0); // Only print on the zero processor | ||
|
||
bool isRHS = false; | ||
if (filenameRHS != "") { | ||
isRHS = true; | ||
} | ||
|
||
// Check to see if the number of right-hand sides is the same as requested. | ||
if (numRHS > 1) { | ||
isRHS = false; // numRHS > 1 not yet supported | ||
X = rcp(new MV(map, numRHS)); | ||
B = rcp(new MV(map, numRHS)); | ||
X->randomize(); | ||
OPT::Apply(*A, *X, *B); // B := AX | ||
X->putScalar(0.0); // annihilate X | ||
} else { | ||
if (isRHS) { | ||
B->print(std::cout); | ||
B = Tpetra::MatrixMarket::Reader<MV>::readVectorFile(filenameRHS, comm, map); | ||
// std::cout << "rhs from input file " << std::endl; | ||
// B->print(std::cout); | ||
X = rcp(new MV(map, numRHS)); | ||
X->scale(0.0); | ||
} else { | ||
LO locNumCol = map->getMaxLocalIndex() + 1; // Create a known solution | ||
GO globNumCol = map->getMaxGlobalIndex() + 1; | ||
for (LO li = 0; li <= locNumCol; li++) { | ||
const auto gid = map->getGlobalElement(li); | ||
ST value = (ST)(globNumCol - 1 - gid); | ||
int numEntries = 1; | ||
vecX->replaceGlobalValue(numEntries, 0, value); | ||
} | ||
A->apply(*vecX, *vecB); // Create a consistent linear system | ||
|
||
// At this point, the initial guess is exact. | ||
bool goodInitGuess = true; // initial guess near solution | ||
bool zeroInitGuess = false; // annihilate initial guess | ||
if (goodInitGuess) { | ||
ST value = 1.e-2; // "Rel RHS Err" and "Rel Mat Err" apply to the residual equation, | ||
// LO numEntries = 1; // norm( b - A x_k ) ?<? relResTol norm( b- Axo). | ||
LO index = 0; // norm(b) is inaccessible to LSQR. | ||
vecX->sumIntoLocalValue(index, 0, value); | ||
} | ||
|
||
if (zeroInitGuess) { | ||
vecX->putScalar(0.0); | ||
} | ||
|
||
X = vecX; | ||
B = vecB; | ||
} | ||
} | ||
|
||
// Create parameter list for the LSQR solver manager | ||
const int numGlobalElements = B->getGlobalLength(); | ||
if (maxiters == -1) | ||
maxiters = numGlobalElements / blockSize - 1; // maximum number of iterations to run | ||
RCP<ParameterList> belosList = rcp(new ParameterList()); | ||
belosList->set("Block Size", blockSize); // LSQR blockSize, must be one | ||
belosList->set("Lambda", damp); // Regularization parameter | ||
belosList->set("Rel RHS Err", relResTol); // Relative convergence tolerance requested | ||
belosList->set("Rel Mat Err", relMatTol); // Maximum number of restarts allowed | ||
belosList->set("Condition Limit", maxCond); // upper bound for cond(A) | ||
belosList->set("Maximum Iterations", maxiters); // Maximum number of iterations allowed | ||
int verbosity = Belos::Errors + Belos::Warnings; | ||
if (verbose) { | ||
verbosity += Belos::TimingDetails + Belos::StatusTestDetails; | ||
if (frequency > 0) | ||
belosList->set("Output Frequency", frequency); | ||
} | ||
if (debug) { | ||
verbosity += Belos::Debug; | ||
} | ||
belosList->set("Verbosity", verbosity); | ||
|
||
// Construct an unpreconditioned linear problem instance. | ||
RCP<LinearProblem> problem = rcp(new LinearProblem(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; | ||
} | ||
// Apply Single Vector LSQR: | ||
// Create an iterative solver manager. | ||
RCP<Solver> newSolver = rcp(new Solver(problem, belosList)); | ||
|
||
if (procVerbose) { | ||
// Print a problem description | ||
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 iterations for a linear system: " << maxiters << std::endl; | ||
std::cout << "Relative residual tolerance: " << relResTol << std::endl; | ||
std::cout << std::endl; | ||
std::cout << "Solver's Description: " << std::endl; | ||
std::cout << newSolver->description() << std::endl; // visually verify the parameter list | ||
} | ||
|
||
// Perform solve | ||
Belos::ReturnType ret = newSolver->solve(); | ||
|
||
std::vector<ST> solNorm(numRHS); // get solution norm | ||
MVT::MvNorm(*X, solNorm); | ||
int numIters = newSolver->getNumIters(); // get number of solver iterations | ||
MT condNum = newSolver->getMatCondNum(); | ||
MT matrixNorm = newSolver->getMatNorm(); | ||
MT resNorm = newSolver->getResNorm(); | ||
MT lsResNorm = newSolver->getMatResNorm(); | ||
|
||
if (procVerbose) | ||
std::cout << "Number of iterations performed for this solve: " << numIters << std::endl | ||
<< "matrix condition number: " << condNum << std::endl | ||
<< "matrix norm: " << matrixNorm << std::endl | ||
<< "residual norm: " << resNorm << std::endl | ||
<< "solution norm: " << solNorm[0] << std::endl | ||
<< "least squares residual Norm: " << lsResNorm << std::endl; | ||
|
||
// Compute actual residuals. | ||
bool badRes = false; | ||
std::vector<ST> actualResids(numRHS); | ||
std::vector<ST> rhsNorm(numRHS); | ||
MV resid(map, numRHS); | ||
OPT::Apply(*A, *X, resid); | ||
MVT::MvAddMv(-1.0, resid, 1.0, *B, resid); | ||
MVT::MvNorm(resid, actualResids); | ||
MVT::MvNorm(*B, rhsNorm); | ||
if (procVerbose) { | ||
std::cout << "---------- Actual Residuals (normalized) ----------" << std::endl << std::endl; | ||
for (int i = 0; i < numRHS; i++) { | ||
ST actRes = actualResids[i] / rhsNorm[i]; | ||
std::cout << "Problem " << i << " : \t" << actRes << std::endl; | ||
if (actRes > relResTol * resGrowthFactor) { | ||
badRes = true; | ||
if (verbose) | ||
std::cout << "residual norm > " << relResTol * resGrowthFactor << std::endl; | ||
} | ||
} | ||
} | ||
|
||
if (ret != Belos::Converged || badRes) { | ||
success = false; | ||
if (procVerbose) | ||
std::cout << std::endl << "ERROR: Belos did not converge!" << std::endl; | ||
} else { | ||
success = true; | ||
if (procVerbose) | ||
std::cout << std::endl << "SUCCESS: Belos converged!" << std::endl; | ||
} | ||
} | ||
TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success); | ||
|
||
return success ? EXIT_SUCCESS : EXIT_FAILURE; | ||
} | ||
|
||
int main(int argc, char *argv[]) { | ||
return run<double>(argc, argv); | ||
// return run<float>(argc, argv); | ||
} // end LSQRTpetraExFile.cpp |
Oops, something went wrong.