43 #ifndef IFPACK2_DATABASESCHWARZ_DEF_HPP
44 #define IFPACK2_DATABASESCHWARZ_DEF_HPP
46 #include "Ifpack2_Parameters.hpp"
48 #include "Tpetra_CrsMatrix.hpp"
49 #include "Teuchos_FancyOStream.hpp"
50 #include "Teuchos_oblackholestream.hpp"
59 template<
class MatrixType>
63 IsInitialized_(false),
74 PatchTolerance_(1e-3),
78 this->setObjectLabel(
"Ifpack2::DatabaseSchwarz");
82 template<
class MatrixType>
87 IsInitialized_(false),
98 PatchTolerance_(1e-3),
102 this->setObjectLabel(
"Ifpack2::DatabaseSchwarz");
107 template<
class MatrixType>
112 template<
class MatrixType>
117 IsInitialized_ =
false;
124 template<
class MatrixType>
129 this->setParametersImpl(const_cast<Teuchos::ParameterList&>(params));
133 template<
class MatrixType>
139 const int defaultPatchSize = 9;
140 const double defaultPatchTolerance = 1e-3;
141 const bool defaultSkipDatabase =
false;
142 const bool defaultVerbose =
false;
145 PatchSize_ = params.
get<
int>(
"database schwarz: patch size",defaultPatchSize);
148 PatchSize_ < 0, std::invalid_argument,
149 "Ifpack2::DatabaseSchwarz::setParameters: \"database schwarz: patch size\" parameter "
150 "must be a nonnegative integer. You gave a value of " << PatchSize_ <<
".");
153 PatchTolerance_ = params.
get(
"database schwarz: patch tolerance", defaultPatchTolerance);
156 PatchTolerance_ <= 0, std::invalid_argument,
157 "Ifpack2::DatabaseSchwarz::setParameters: \"database schwarz: patch tolerance\" parameter "
158 "must be a positive double. You gave a value of " << PatchTolerance_ <<
".");
161 SkipDatabase_ = params.
get<
bool>(
"database schwarz: skip database",defaultSkipDatabase);
164 Verbose_ = params.
get<
bool>(
"database schwarz: print database summary",defaultVerbose);
168 template<
class MatrixType>
172 (void) zeroStartingSolution;
175 template<
class MatrixType>
181 A.
is_null(), std::runtime_error,
"Ifpack2::DatabaseSchwarz::getComm: The input "
182 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
183 "before calling this method.");
184 return A->getRowMap()->getComm();
188 template<
class MatrixType>
196 template<
class MatrixType>
197 Teuchos::RCP<
const Tpetra::CrsMatrix<
typename MatrixType::scalar_type,
198 typename MatrixType::local_ordinal_type,
199 typename MatrixType::global_ordinal_type,
200 typename MatrixType::node_type> >
205 return Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (getMatrix());
209 template<
class MatrixType>
216 A.
is_null(), std::runtime_error,
"Ifpack2::DatabaseSchwarz::getDomainMap: The "
217 "input matrix A is null. Please call setMatrix() with a nonnull input "
218 "matrix before calling this method.");
219 return A->getDomainMap();
223 template<
class MatrixType>
230 A.
is_null(), std::runtime_error,
"Ifpack2::DatabaseSchwarz::getRangeMap: The "
231 "input matrix A is null. Please call setMatrix() with a nonnull input "
232 "matrix before calling this method.");
233 return A->getRangeMap();
237 template<
class MatrixType>
243 template<
class MatrixType>
245 return NumInitialize_;
249 template<
class MatrixType>
255 template<
class MatrixType>
261 template<
class MatrixType>
263 return InitializeTime_;
267 template<
class MatrixType>
273 template<
class MatrixType>
279 template<
class MatrixType>
281 return ComputeFlops_;
285 template<
class MatrixType>
290 template<
class MatrixType>
294 A.
is_null(), std::runtime_error,
"Ifpack2::DatabaseSchwarz::getNodeSmootherComplexity: "
295 "The input matrix A is null. Please call setMatrix() with a nonnull "
296 "input matrix, then call compute(), before calling this method.");
298 return A->getLocalNumRows() + A->getLocalNumEntries();
303 template<
class MatrixType>
306 apply(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
307 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
312 const std::string timerName (
"Ifpack2::DatabaseSchwarz::apply");
318 double startTime = timer->
wallTime();
327 ! isComputed(), std::runtime_error,
328 "Ifpack2::DatabaseSchwarz::apply(): You must call the compute() method before "
329 "you may call apply().");
331 X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
332 "Ifpack2::DatabaseSchwarz::apply(): X and Y must have the same number of "
333 "columns. X.getNumVectors() = " << X.getNumVectors() <<
" != "
334 <<
"Y.getNumVectors() = " << Y.getNumVectors() <<
".");
340 auto X_view = X.getLocalViewHost(Tpetra::Access::ReadOnly);
341 auto Y_view = Y.getLocalViewHost(Tpetra::Access::ReadWrite);
345 for(
unsigned int ipatch=0; ipatch<NumPatches_; ipatch++) {
346 int idatabase = DatabaseIndices_[ipatch];
350 for(
unsigned int c=0; c<X_view.extent(1); ++c) {
351 for(
unsigned int i=0; i<x_patch.
size(); ++i) {
352 x_patch[i] = X_view(PatchIndices_[ipatch][i],c);
361 int* ipiv = &ipiv_[idatabase*PatchSize_];
363 lapack.
GETRS(
'N', DatabaseMatrices_[idatabase]->numRows(), numRhs,
364 DatabaseMatrices_[idatabase]->values(), DatabaseMatrices_[idatabase]->numRows(),
366 DatabaseMatrices_[idatabase]->numRows(), &INFO);
370 INFO < 0, std::logic_error,
"Ifpack2::DatabaseSchwarz::compute: "
371 "LAPACK's _GETRF (LU factorization with partial pivoting) was called "
372 "incorrectly. INFO = " << INFO <<
" < 0. "
373 "Please report this bug to the Ifpack2 developers.");
378 INFO > 0, std::runtime_error,
"Ifpack2::DatabaseSchwarz::compute: "
379 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the "
380 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") "
381 "(one-based index i) is exactly zero. This probably means that the input "
382 "patch is singular.");
385 for(
unsigned int c=0; c<Y_view.extent(1); ++c) {
386 for(
unsigned int i=0; i<x_patch.
size(); ++i) {
387 Y_view(PatchIndices_[ipatch][i],c) += alpha*Weights_[PatchIndices_[ipatch][i]]*x_patch[i];
393 ApplyTime_ += (timer->
wallTime() - startTime);
397 template<
class MatrixType>
400 applyMat (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
401 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
405 X.getNumVectors() != Y.getNumVectors(), std::invalid_argument,
406 "Ifpack2::DatabaseSchwarz::applyMat: X.getNumVectors() != Y.getNumVectors().");
410 A.is_null(), std::runtime_error,
"Ifpack2::DatabaseSchwarz::applyMat: The input "
411 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
412 "before calling this method.");
414 A->apply (X, Y, mode);
418 template<
class MatrixType>
423 const std::string timerName (
"Ifpack2::DatabaseSchwarz::initialize");
428 IsInitialized_ =
true;
433 template<
class MatrixType>
436 const std::string timerName (
"Ifpack2::DatabaseSchwarz::compute");
442 double startTime = timer->
wallTime();
447 if (!isInitialized()) {
451 const int maxNnzPerRow = A_->getGlobalMaxNumRowEntries();
454 PatchIndices_.resize(0);
458 size_t num_entries = A_->getNumEntriesInLocalRow(irow);
463 typename row_matrix_type::nonconst_local_inds_host_view_type row_inds(
"row indices", maxNnzPerRow);
464 typename row_matrix_type::nonconst_values_host_view_type row_vals(
"row values", maxNnzPerRow);
465 A_->getLocalRowCopy(irow, row_inds, row_vals, num_entries);
468 bool is_new_patch =
true;
469 for(
size_t ipatch=0; ipatch<PatchIndices_.size(); ++ipatch) {
470 for(
size_t i=0; i<PatchIndices_[ipatch].size(); ++i) {
471 if(PatchIndices_[ipatch][i] == irow) {
472 is_new_patch =
false;
473 ipatch=PatchIndices_.size();
482 std::vector<typename row_matrix_type::local_ordinal_type> indices_vector = Teuchos::createVector(indices_array_view);
483 PatchIndices_.push_back(indices_vector);
492 DatabaseIndices_.resize(NumPatches_,-1);
493 Weights_.resize(A_->getLocalNumRows(),0);
494 std::vector<double> index_count(A_->getLocalNumRows(),0);
496 for(
unsigned int ipatch=0; ipatch< NumPatches_; ++ipatch) {
498 DenseMatRCP patch_matrix =
Teuchos::rcp(
new DenseMatType(PatchSize_, PatchSize_));
499 auto indices_vector = PatchIndices_[ipatch];
502 index_count[indices_vector[i]]++;
504 typename row_matrix_type::nonconst_local_inds_host_view_type row_inds(
"row indices", maxNnzPerRow);
505 typename row_matrix_type::nonconst_values_host_view_type row_vals(
"row values", maxNnzPerRow);
507 A_->getLocalRowCopy(indices_vector[i], row_inds, row_vals, num_entries);
509 for(
size_t k=0; k<num_entries; ++k) {
510 if(row_inds(k) == indices_vector[j]) {
511 (*patch_matrix)(i,j) = row_vals(k);
519 bool found_match =
false;
521 for(
size_t idatabase=0; idatabase<DatabaseMatrices_.size(); ++idatabase) {
527 DenseMatRCP database_candidate = DatabaseMatrices_[idatabase];
537 DatabaseIndices_[ipatch] = idatabase;
546 DatabaseMatrices_.push_back(patch_matrix);
547 DatabaseIndices_[ipatch] = DatabaseMatrices_.size()-1;
549 "Ifpack2::DatabaseSchwarz::compute: A matrix was added to the database, but appears to be null!"
550 "Please report this bug to the Ifpack2 developers.");
555 for(
unsigned int i=0; i<index_count.size(); ++i) {
557 "Ifpack2::DatabaseSchwarz::compute: DOF " << i <<
" has no corresponding patch! "
558 "All DOFs must be able to locate in a patch to use this method! "
559 "Have you considered adjusting the patch size? Currently it is " << PatchSize_ <<
".");
560 Weights_[i] = 1./index_count[i];
562 DatabaseSize_ = DatabaseMatrices_.size();
565 std::vector<int> database_counts(DatabaseSize_,0);
566 for(
unsigned int ipatch=0; ipatch<NumPatches_; ++ipatch) {
567 database_counts[DatabaseIndices_[ipatch]]++;
573 ipiv_.resize(DatabaseSize_*PatchSize_);
574 std::fill(ipiv_.begin (), ipiv_.end (), 0);
575 for(
unsigned int idatabase=0; idatabase<DatabaseSize_; idatabase++) {
576 int* ipiv = &ipiv_[idatabase*PatchSize_];
578 lapack.
GETRF(DatabaseMatrices_[idatabase]->numRows(),
579 DatabaseMatrices_[idatabase]->numCols(),
580 DatabaseMatrices_[idatabase]->values(),
581 DatabaseMatrices_[idatabase]->stride(),
586 INFO < 0, std::logic_error,
"Ifpack2::DatabaseSchwarz::compute: "
587 "LAPACK's _GETRF (LU factorization with partial pivoting) was called "
588 "incorrectly. INFO = " << INFO <<
" < 0. "
589 "Please report this bug to the Ifpack2 developers.");
594 std::cout <<
"SINGULAR LOCAL MATRIX, COUNT=" << database_counts[idatabase] << std::endl;
597 INFO > 0, std::runtime_error,
"Ifpack2::DatabaseSchwarz::compute: "
598 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the "
599 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") "
600 "(one-based index i) is exactly zero. This probably means that the input "
601 "patch is singular.");
607 ComputeTime_ += (timer->
wallTime() - startTime);
611 std::cout <<
"Ifpack2::DatabaseSchwarz()::Compute() summary\n";
612 std::cout <<
"Found " << NumPatches_ <<
" patches of size " << PatchSize_ <<
" in matrix A\n";
613 std::cout <<
"Database tol = " << PatchTolerance_ <<
"\n";
614 std::cout <<
"Database size = " << DatabaseSize_ <<
" patches\n";
619 template <
class MatrixType>
621 std::ostringstream out;
626 out <<
"\"Ifpack2::DatabaseSchwarz\": {";
627 out <<
"Initialized: " << (isInitialized() ?
"true" :
"false")
628 <<
", Computed: " << (isComputed() ?
"true" :
"false")
629 <<
", patch size: " << PatchSize_
630 <<
", patch tolerance: " << PatchTolerance_
631 <<
", skip database: " << (SkipDatabase_ ?
"true" :
"false")
632 <<
", print database summary: " << (Verbose_ ?
"true" :
"false");
634 if (getMatrix().is_null()) {
635 out <<
"Matrix: null";
638 out <<
"Global matrix dimensions: ["
639 << getMatrix()->getGlobalNumRows() <<
", "
640 << getMatrix()->getGlobalNumCols() <<
"]"
641 <<
", Global nnz: " << getMatrix()->getGlobalNumEntries();
649 template <
class MatrixType>
672 const int myRank = this->getComm()->getRank();
676 out <<
"\"Ifpack2::DatabaseSchwarz\":" << endl;
681 out <<
"Template parameters:" << endl;
684 out <<
"Scalar: " << TypeNameTraits<scalar_type>::name() << endl
685 <<
"LocalOrdinal: " << TypeNameTraits<local_ordinal_type>::name() << endl
686 <<
"GlobalOrdinal: " << TypeNameTraits<global_ordinal_type>::name() << endl
687 <<
"Device: " << TypeNameTraits<device_type>::name() << endl;
689 out <<
"Initialized: " << (isInitialized() ?
"true" :
"false") << endl
690 <<
"Computed: " << (isComputed() ?
"true" :
"false") << endl;
692 if (getMatrix().is_null()) {
693 out <<
"Matrix: null" << endl;
696 out <<
"Global matrix dimensions: ["
697 << getMatrix()->getGlobalNumRows() <<
", "
698 << getMatrix()->getGlobalNumCols() <<
"]" << endl
699 <<
"Global nnz: " << getMatrix()->getGlobalNumEntries() << endl;
708 #define IFPACK2_DATABASESCHWARZ_INSTANT(S,LO,GO,N) \
709 template class Ifpack2::DatabaseSchwarz< Tpetra::RowMatrix<S, LO, GO, N> >;
711 #endif // IFPACK2_DATABASESCHWARZ_DEF_HPP
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:190
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_DatabaseSchwarz_decl.hpp:118
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, returning the result in Y. Y = alpha*Op(A)*X + beta*Y.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:306
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:226
double getApplyTime() const
The total time spent in all calls to apply().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:274
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:291
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:256
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:419
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:238
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:262
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to a Teuchos::FancyOStream.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:651
void setParameters(const Teuchos::ParameterList ¶ms)
Set (or reset) parameters.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:126
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:177
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:212
virtual ~DatabaseSchwarz()
Destructor.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:108
double getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:268
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:286
void GETRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, OrdinalType *info) const
void compute()
(Re)compute the left scaling, and (if applicable) estimate max and min eigenvalues of D_inv * A...
Definition: Ifpack2_DatabaseSchwarz_def.hpp:434
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
DatabaseSchwarz(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:61
Teuchos::RCP< const Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getCrsMatrix() const
Attempt to return the matrix A as a Tpetra::CrsMatrix.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:202
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_DatabaseSchwarz_decl.hpp:127
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:113
void setZeroStartingSolution(bool zeroStartingSolution)
Set this preconditioner's parameters.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:170
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:280
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_DatabaseSchwarz_def.hpp:620
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:244
static magnitudeType magnitude(T a)
Overlapping Schwarz where redundant patches are not stored explicitly.
Definition: Ifpack2_DatabaseSchwarz_decl.hpp:97
void applyMat(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Compute Y = Op(A)*X, where Op(A) is either A, , or .
Definition: Ifpack2_DatabaseSchwarz_def.hpp:400
void GETRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
TypeTo as(const TypeFrom &t)
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_DatabaseSchwarz_decl.hpp:115
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_DatabaseSchwarz_decl.hpp:121
int getNumCompute() const
The total number of successful calls to compute().
Definition: Ifpack2_DatabaseSchwarz_def.hpp:250