Skip to content

Commit

Permalink
Merge branch 'develop' into adelus-add-old-factorsolve-interfaces
Browse files Browse the repository at this point in the history
  • Loading branch information
@vqd8a
vqd8a committed Jan 3, 2023
2 parents 4366811 + d231f37 commit 340dded
Show file tree
Hide file tree
Showing 108 changed files with 646 additions and 370 deletions.
1 change: 0 additions & 1 deletion packages/belos/kokkos/example/BlockCGKokkosExFile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,6 @@ bool success = true;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Belos::KokkosMultiVec<ST, EXSP> MV;
typedef Belos::KokkosCrsOperator<ST, OT, EXSP> OP;
typedef Belos::MultiVec<ST> KMV;
typedef Belos::Operator<ST> KOP;
typedef Belos::MultiVecTraits<ST,KMV> MVT;
Expand Down
1 change: 0 additions & 1 deletion packages/belos/kokkos/example/BlockGmresKokkosExFile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,6 @@ bool success = true;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Belos::KokkosMultiVec<ST, EXSP> MV;
typedef Belos::KokkosCrsOperator<ST, OT, EXSP> OP;
typedef Belos::MultiVec<ST> KMV;
typedef Belos::Operator<ST> KOP;
typedef Belos::MultiVecTraits<ST,KMV> MVT;
Expand Down
3 changes: 1 addition & 2 deletions packages/belos/src/BelosMinresSolMgr.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ namespace Belos {
typedef OperatorTraits<ScalarType,MV,OP> OPT;
typedef Teuchos::ScalarTraits<ScalarType> SCT;
typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
typedef Teuchos::ScalarTraits< MagnitudeType > MT;
typedef Teuchos::ScalarTraits< MagnitudeType > MST;

public:

Expand Down Expand Up @@ -408,7 +408,6 @@ namespace Belos {
using Teuchos::rcpFromRef;
using Teuchos::EnhancedNumberValidator;
typedef MagnitudeType MT;
typedef Teuchos::ScalarTraits<MT> MST;

// List of parameters accepted by MINRES, and their default values.
RCP<ParameterList> pl = parameterList ("MINRES");
Expand Down
2 changes: 0 additions & 2 deletions packages/belos/tpetra/src/BelosTpetraTestFramework.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,7 +679,6 @@ namespace Belos {
typedef local_ordinal_type LO;
typedef global_ordinal_type GO;
typedef node_type NT;
typedef ::Tpetra::Map<LO, GO, NT> map_type;

// For a square matrix, we only need a Map for the range of the matrix.
RCP<const map_type> pRangeMap =
Expand Down Expand Up @@ -730,7 +729,6 @@ namespace Belos {
typedef typename SparseMatrixType::local_ordinal_type LO;
typedef typename SparseMatrixType::global_ordinal_type GO;
typedef typename SparseMatrixType::node_type NT;
typedef ::Tpetra::Map<LO, GO, NT> map_type;

// For a square matrix, we only need a Map for the range of the matrix.
RCP<const map_type> pRangeMap =
Expand Down
46 changes: 23 additions & 23 deletions packages/belos/tpetra/src/solvers/Belos_Tpetra_Cg.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,19 +74,19 @@ class Cg: public Krylov <SC, MV, OP> {
const SolverInput<SC>& input)
{
using std::endl;
using device_type = typename MV::device_type;
using dev_type = typename MV::device_type;
using val_type = typename MV::dot_type;
using STS = Kokkos::ArithTraits<val_type>;
using mag_type = typename STS::mag_type;
using STM = Kokkos::ArithTraits<mag_type>;
using ATS = Kokkos::ArithTraits<val_type>;
using magnitude_type = typename ATS::mag_type;
using ATM = Kokkos::ArithTraits<magnitude_type>;

const auto ONE = STS::one ();
const auto ONE = ATS::one ();

SolverOutput<SC> output {};

// scalars
mag_type PAP;
mag_type alpha;
magnitude_type PAP;
magnitude_type alpha;

vec_type P (B.getMap ());
vec_type AP (B.getMap ());
Expand All @@ -95,17 +95,17 @@ class Cg: public Krylov <SC, MV, OP> {
vec_type MR = * (R_MR.getVectorNonConst (1));

// compute initial residual
Kokkos::View<val_type*, device_type> r_beta ("results[numVecs]", 2);
mag_type beta_old = STM::zero ();
mag_type r_norm = STM::zero ();
mag_type r_norm_orig = STM::zero ();
Kokkos::View<val_type*, dev_type> r_beta ("results[numVecs]", 2);
magnitude_type beta_old = ATM::zero ();
magnitude_type r_norm = ATM::zero ();
magnitude_type r_norm_orig = ATM::zero ();

A.apply (X, R);
R.update (ONE, B, -ONE);

if (input.precoSide == "none") { // no preconditioner
Tpetra::deep_copy (P, R);
beta_old = STS::real (R.dot (R));
beta_old = ATS::real (R.dot (R));
r_norm = beta_old;
}
else {
Expand All @@ -117,29 +117,29 @@ class Cg: public Krylov <SC, MV, OP> {
//TODO: idot is used for now.
auto req = Tpetra::idot (r_beta, R_MR, R);
req->wait ();
r_norm = STS::real (r_beta(0));
beta_old = STS::real (r_beta(1));
r_norm = ATS::real (r_beta(0));
beta_old = ATS::real (r_beta(1));
}
r_norm = std::sqrt (r_norm);
r_norm_orig = r_norm;

// quick-return
mag_type metric =
magnitude_type metric =
this->getConvergenceMetric (r_norm, r_norm_orig, input);
if (metric <= input.tol) {
if (outPtr != nullptr) {
*outPtr << "Initial guess' residual norm " << r_norm
<< " meets tolerance " << input.tol << endl;
}
output.absResid = r_norm;
output.relResid = STM::one ();
output.relResid = ATM::one ();
output.numIters = 0;
output.converged = true;
return output;
}

// main loop
mag_type beta_new = STM::zero ();
magnitude_type beta_new = ATM::zero ();
for (int iter = 0; iter < input.maxNumIters; ++iter) {
if (outPtr != nullptr) {
*outPtr << "Iteration " << (iter+1) << " of " << input.maxNumIters << ":" << endl;
Expand All @@ -148,13 +148,13 @@ class Cg: public Krylov <SC, MV, OP> {
}

A.apply (P, AP); // AP = A*P
PAP = STS::real (P.dot (AP));
PAP = ATS::real (P.dot (AP));

if (outPtr != nullptr) {
*outPtr << "PAP: " << PAP << endl;
}
TEUCHOS_TEST_FOR_EXCEPTION
(STS::real (PAP) <= STM::zero (), std::runtime_error,
(ATS::real (PAP) <= ATM::zero (), std::runtime_error,
"At iteration " << (iter+1) << " out of " << input.maxNumIters
<< ", P.dot(AP) = " << PAP << " <= 0. This usually means that "
"the matrix A is not symmetric (Hermitian) positive definite.");
Expand All @@ -168,16 +168,16 @@ class Cg: public Krylov <SC, MV, OP> {
R.update (static_cast<SC> (-alpha), AP, ONE); // R = R - alpha*AP

if (input.precoSide == "none") { // no preconditioner
beta_new = STS::real (R.dot (R));
beta_new = ATS::real (R.dot (R));
r_norm = beta_new;
}
else {
M.apply (R, MR);
//TODO: idot is used to compute [MR, R]'*[R] for now.
auto req = Tpetra::idot (r_beta, R_MR, R);
req->wait ();
r_norm = STS::real (r_beta(0));
beta_new = STS::real (r_beta(1));
r_norm = ATS::real (r_beta(0));
beta_new = ATS::real (r_beta(1));
}
r_norm = std::sqrt (r_norm);
metric = this->getConvergenceMetric (r_norm, r_norm_orig, input);
Expand All @@ -195,7 +195,7 @@ class Cg: public Krylov <SC, MV, OP> {
return output;
}
else if (iter + 1 < input.maxNumIters) { // not last iteration
const mag_type beta = beta_new / beta_old;
const magnitude_type beta = beta_new / beta_old;
if (input.precoSide == "none") {
P.update (ONE, R, static_cast<SC> (beta)); // P += beta*R
}
Expand Down
36 changes: 18 additions & 18 deletions packages/belos/tpetra/src/solvers/Belos_Tpetra_CgPipeline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,13 +74,13 @@ class CgPipeline : public Krylov<SC, MV, OP> {
const SolverInput<SC>& input) override
{
using std::endl;
using device_type = typename MV::device_type;
using dev_type = typename MV::device_type;
using val_type = typename MV::dot_type;
using STS = Kokkos::ArithTraits<val_type>;
using mag_type = typename STS::mag_type;
using STM = Kokkos::ArithTraits<mag_type>;
using ATS = Kokkos::ArithTraits<val_type>;
using magnitude_type = typename ATS::mag_type;
using ATM = Kokkos::ArithTraits<magnitude_type>;

const auto ONE = STS::one ();
const auto ONE = ATS::one ();

SolverOutput<SC> output {};

Expand All @@ -96,7 +96,7 @@ class CgPipeline : public Krylov<SC, MV, OP> {
vec_type R2 = * (R_R.getVectorNonConst (1));

// results of [R R]'*[R AR]
Kokkos::View<val_type*, device_type> RR_RAR ("results[numVecs]", 2);
Kokkos::View<val_type*, dev_type> RR_RAR ("results[numVecs]", 2);
auto RR_RAR_host = Kokkos::create_mirror_view (RR_RAR);
vec_type P (R, Teuchos::Copy);
vec_type AP (P.getMap ());
Expand All @@ -108,12 +108,12 @@ class CgPipeline : public Krylov<SC, MV, OP> {

// local vars
val_type RAR, PAP;
mag_type alpha = STM::zero ();
mag_type beta = STM::zero ();
mag_type r_norm = STM::zero ();
mag_type r_norm_orig = STM::zero ();
mag_type beta_new = STM::zero ();
mag_type beta_old = STM::zero ();
magnitude_type alpha = ATM::zero ();
magnitude_type beta = ATM::zero ();
magnitude_type r_norm = ATM::zero ();
magnitude_type r_norm_orig = ATM::zero ();
magnitude_type beta_new = ATM::zero ();
magnitude_type beta_old = ATM::zero ();

// Initial step
if (input.precoSide == "none") {
Expand Down Expand Up @@ -162,13 +162,13 @@ class CgPipeline : public Krylov<SC, MV, OP> {
req->wait ();
Kokkos::deep_copy (RR_RAR_host, RR_RAR);
RAR = RR_RAR_host(1);
beta_new = STS::real (RR_RAR_host(0));
beta_new = ATS::real (RR_RAR_host(0));

r_norm = std::sqrt( beta_new );
if (iter == 0) {
r_norm_orig = r_norm;
}
const mag_type metric =
const magnitude_type metric =
this->getConvergenceMetric (r_norm, r_norm_orig, input);
if (outPtr != nullptr) {
*outPtr << ", RAR: " << RAR << ", r_norm: " << r_norm << ", metric: " << metric;
Expand All @@ -187,8 +187,8 @@ class CgPipeline : public Krylov<SC, MV, OP> {
}

if (iter == 0) {
alpha = beta_new / STS::real (RAR);
beta = STM::zero ();
alpha = beta_new / ATS::real (RAR);
beta = ATM::zero ();
}
else {
// beta
Expand All @@ -199,14 +199,14 @@ class CgPipeline : public Krylov<SC, MV, OP> {
// PAP
PAP = RAR - beta_new * (beta / alpha);
TEUCHOS_TEST_FOR_EXCEPTION
(STS::real (RAR) <= STM::zero (), std::runtime_error,
(ATS::real (RAR) <= ATM::zero (), std::runtime_error,
"At iteration " << (iter+1)
<< " out of " << input.maxNumIters << ", R.dot(AR) = " << RAR <<
" <= 0. This usually means that the matrix A is not symmetric "
"(Hermitian) positive definite.");

// alpha
alpha = beta_new / STS::real (PAP);
alpha = beta_new / ATS::real (PAP);
}
if (outPtr != nullptr) {
*outPtr << ", alpha: " << alpha << endl;
Expand Down
48 changes: 24 additions & 24 deletions packages/belos/tpetra/src/solvers/Belos_Tpetra_CgSingleReduce.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,33 +74,33 @@ class CgSingleReduce: public Krylov<SC, MV, OP> {
const SolverInput<SC>& input) override
{
using std::endl;
using device_type = typename MV::device_type;
using STS = Teuchos::ScalarTraits<SC>;
using mag_type = typename STS::magnitudeType;
using STM = Teuchos::ScalarTraits<mag_type>;
using dev_type = typename MV::device_type;
using ATS = Kokkos::ArithTraits<SC>;
using magnitude_type = typename ATS::mag_type;
using ATM = Kokkos::ArithTraits<magnitude_type>;
using dot_type = typename MV::dot_type;

const SC ONE = STS::one ();
const SC ONE = ATS::one ();
SolverOutput<SC> output {};

// compute initial residual
vec_type MR (B.getMap ());
mag_type beta_old = STM::zero ();
magnitude_type beta_old = ATM::zero ();
if (input.precoSide == "none") {
beta_old = STS::real (B.dot (B));
beta_old = ATS::real (B.dot (B));
}
else {
M.apply (B, MR);
beta_old = STS::real (B.dot (MR));
beta_old = ATS::real (B.dot (MR));
}
mag_type r_norm = std::sqrt (beta_old);
mag_type r_norm_orig = r_norm;
magnitude_type r_norm = std::sqrt (beta_old);
magnitude_type r_norm_orig = r_norm;

// quick return
mag_type metric = this->getConvergenceMetric (r_norm, r_norm_orig, input);
magnitude_type metric = this->getConvergenceMetric (r_norm, r_norm_orig, input);
if (metric <= input.tol) {
output.absResid = r_norm;
output.relResid = STM::one ();
output.relResid = ATM::one ();
output.numIters = 0;
output.converged = true;
// R doesn't exist yet, so we don't have to copy R back to B
Expand All @@ -115,18 +115,18 @@ class CgSingleReduce: public Krylov<SC, MV, OP> {
vec_type AR = * (R_AR.getVectorNonConst (1));
Tpetra::deep_copy (R, B);
// results of [R AR]'*R
mag_type RAR;
Kokkos::View<dot_type*, device_type> RR_RAR ("results[numVecs]", 2);
magnitude_type RAR;
Kokkos::View<dot_type*, dev_type> RR_RAR ("results[numVecs]", 2);
auto RR_RAR_host = Kokkos::create_mirror_view (RR_RAR);
vec_type P (R, Teuchos::Copy);
vec_type AP (P.getMap ());

// Initial step
// AR = A*R
mag_type PAP;
magnitude_type PAP;
if (input.precoSide == "none") {
A.apply (R, AR);
PAP = STS::real (R.dot (AR));
PAP = ATS::real (R.dot (AR));
}
else {
M.apply (R, MR);
Expand All @@ -140,14 +140,14 @@ class CgSingleReduce: public Krylov<SC, MV, OP> {
req->wait ();

Kokkos::deep_copy (RR_RAR_host, RR_RAR);
beta_old = STS::real (RR_RAR_host(0));
PAP = STS::real (RR_RAR_host(1));
beta_old = ATS::real (RR_RAR_host(0));
PAP = ATS::real (RR_RAR_host(1));

r_norm = std::sqrt (beta_old);
}
mag_type alpha = beta_old / PAP;
mag_type beta = STM::zero ();
mag_type beta_new = STM::zero ();
magnitude_type alpha = beta_old / PAP;
magnitude_type beta = ATM::zero ();
magnitude_type beta_new = ATM::zero ();
// main loop
for (int iter = 0; iter < input.maxNumIters; ++iter) {
if (outPtr != nullptr) {
Expand Down Expand Up @@ -191,8 +191,8 @@ class CgSingleReduce: public Krylov<SC, MV, OP> {
}
// * all-reduce *
Kokkos::deep_copy (RR_RAR_host, RR_RAR);
beta_new = STS::real (RR_RAR_host(0));
RAR = STS::real (RR_RAR_host(1));
beta_new = ATS::real (RR_RAR_host(0));
RAR = ATS::real (RR_RAR_host(1));

// convergence check
r_norm = std::sqrt( beta_new );
Expand Down Expand Up @@ -221,7 +221,7 @@ class CgSingleReduce: public Krylov<SC, MV, OP> {
// PAP
PAP = RAR - beta_new * (beta /alpha);
TEUCHOS_TEST_FOR_EXCEPTION
(RAR <= STM::zero (), std::runtime_error, "At iteration " << (iter+1)
(RAR <= ATM::zero (), std::runtime_error, "At iteration " << (iter+1)
<< " out of " << input.maxNumIters << ", R.dot(AR) = " << RAR <<
" <= 0. This usually means that the matrix A is not symmetric "
"(Hermitian) positive definite.");
Expand Down
3 changes: 1 addition & 2 deletions packages/ifpack2/src/Ifpack2_BlockRelaxation_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -698,9 +698,8 @@ initialize ()
// only needed when Schwarz combine mode is ADD as opposed to ZERO (which is RAS)

if (schwarzCombineMode_ == "ADD") {
typedef Teuchos::RCP<Tpetra::Import<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> const > import_type;
typedef Tpetra::MultiVector< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type,typename MatrixType::node_type> scMV;
import_type theImport = A_->getGraph()->getImporter();
Teuchos::RCP<const import_type> theImport = A_->getGraph()->getImporter();
if (!theImport.is_null()) {
scMV nonOverLapW(theImport->getSourceMap(), 1, false);
Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
Expand Down
2 changes: 1 addition & 1 deletion packages/ifpack2/src/Ifpack2_BlockTriDiContainer_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ namespace Ifpack2 {
impl_ = Teuchos::rcp(new BlockTriDiContainerDetails::ImplObject<MatrixType>());

using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
using block_crs_matrix_type = typename impl_type::tpetra_block_crs_matrix_type;
// using block_crs_matrix_type = typename impl_type::tpetra_block_crs_matrix_type;

impl_->A = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(matrix);
TEUCHOS_TEST_FOR_EXCEPT_MSG
Expand Down
Loading

0 comments on commit 340dded

Please sign in to comment.