43 #ifndef IFPACK2_TRIDICONTAINER_DEF_HPP 44 #define IFPACK2_TRIDICONTAINER_DEF_HPP 47 #include "Teuchos_LAPACK.hpp" 51 # include "Teuchos_DefaultMpiComm.hpp" 53 # include "Teuchos_DefaultSerialComm.hpp" 59 template<
class MatrixType,
class LocalScalarType>
62 const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
63 Container<MatrixType> (matrix, localRows),
64 numRows_ (localRows.size ()),
65 diagBlock_ (numRows_, numRows_),
69 using Teuchos::ArrayView;
72 using Teuchos::toString;
73 typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
74 typedef typename ArrayView<const local_ordinal_type>::size_type size_type;
75 TEUCHOS_TEST_FOR_EXCEPTION(
76 ! matrix->hasColMap (), std::invalid_argument,
"Ifpack2::TriDiContainer: " 77 "The constructor's input matrix must have a column Map.");
80 const map_type& rowMap = * (matrix->getRowMap ());
81 const size_type numRows = localRows.size ();
82 bool rowIndicesValid =
true;
83 Array<local_ordinal_type> invalidLocalRowIndices;
84 for (size_type i = 0; i < numRows; ++i) {
85 if (! rowMap.isNodeLocalElement (localRows[i])) {
86 rowIndicesValid =
false;
87 invalidLocalRowIndices.push_back (localRows[i]);
91 TEUCHOS_TEST_FOR_EXCEPTION(
92 ! rowIndicesValid, std::invalid_argument,
"Ifpack2::TriDiContainer: " 93 "On process " << rowMap.getComm ()->getRank () <<
" of " 94 << rowMap.getComm ()->getSize () <<
", in the given set of local row " 95 "indices localRows = " <<
toString (localRows) <<
", the following " 96 "entries are not valid local row indices on the calling process: " 97 <<
toString (invalidLocalRowIndices) <<
".");
100 RCP<const Teuchos::Comm<int> > localComm =
101 rcp (
new Teuchos::MpiComm<int> (MPI_COMM_SELF));
103 RCP<const Teuchos::Comm<int> > localComm =
104 rcp (
new Teuchos::SerialComm<int> ());
110 localMap_ = rcp (
new map_type (numRows_, indexBase, localComm));
113 template<
class MatrixType,
class LocalScalarType>
116 const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
117 Container<MatrixType> (matrix, localRows)
119 TEUCHOS_TEST_FOR_EXCEPTION
120 (
true, std::logic_error,
"Ifpack2::TriDiContainer: Not implemented for " 121 "LocalScalarType = " << Teuchos::TypeNameTraits<LocalScalarType>::name ()
125 template<
class MatrixType,
class LocalScalarType>
129 template<
class MatrixType,
class LocalScalarType>
134 template<
class MatrixType,
class LocalScalarType>
140 template<
class MatrixType,
class LocalScalarType>
147 template<
class MatrixType,
class LocalScalarType>
150 return IsInitialized_;
153 template<
class MatrixType,
class LocalScalarType>
156 return IsInitialized_;
159 template<
class MatrixType,
class LocalScalarType>
165 template<
class MatrixType,
class LocalScalarType>
171 template<
class MatrixType,
class LocalScalarType>
178 template<
class MatrixType,
class LocalScalarType>
185 template<
class MatrixType,
class LocalScalarType>
193 IsInitialized_ =
false;
197 diagBlock_.putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
198 std::fill (ipiv_.begin (), ipiv_.end (), 0);
200 IsInitialized_ =
true;
203 template<
class MatrixType,
class LocalScalarType>
208 template<
class MatrixType,
class LocalScalarType>
211 TEUCHOS_TEST_FOR_EXCEPTION(
212 static_cast<size_t> (ipiv_.size ()) != numRows_, std::logic_error,
213 "Ifpack2::TriDiContainer::compute: ipiv_ array has the wrong size. " 214 "Please report this bug to the Ifpack2 developers.");
217 if (! this->isInitialized ()) {
222 extract (this->getMatrix ());
228 template<
class MatrixType,
class LocalScalarType>
233 template<
class MatrixType,
class LocalScalarType>
236 Teuchos::LAPACK<int, local_scalar_type> lapack;
238 lapack.GTTRF (diagBlock_.numRowsCols (),
243 ipiv_.getRawPtr (), &INFO);
245 TEUCHOS_TEST_FOR_EXCEPTION(
246 INFO < 0, std::logic_error,
"Ifpack2::TriDiContainer::factor: " 247 "LAPACK's _GTTRF (LU factorization with partial pivoting) was called " 248 "incorrectly. INFO = " << INFO <<
" < 0. " 249 "Please report this bug to the Ifpack2 developers.");
253 TEUCHOS_TEST_FOR_EXCEPTION(
254 INFO > 0, std::runtime_error,
"Ifpack2::TriDiContainer::factor: " 255 "LAPACK's _GTTRF (LU factorization with partial pivoting) reports that the " 256 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") " 257 "(one-based index i) is exactly zero. This probably means that the input " 258 "matrix has a singular diagonal block.");
261 template<
class MatrixType,
class LocalScalarType>
266 template<
class MatrixType,
class LocalScalarType>
270 Teuchos::ETransp mode,
271 LocalScalarType alpha,
272 LocalScalarType beta)
const 274 using Teuchos::ArrayRCP;
277 using Teuchos::rcpFromRef;
279 TEUCHOS_TEST_FOR_EXCEPTION(
280 X.getLocalLength () != Y.getLocalLength (),
281 std::logic_error,
"Ifpack2::TriDiContainer::applyImpl: X and Y have " 282 "incompatible dimensions (" << X.getLocalLength () <<
" resp. " 283 << Y.getLocalLength () <<
"). Please report this bug to " 284 "the Ifpack2 developers.");
285 TEUCHOS_TEST_FOR_EXCEPTION(
286 localMap_->getNodeNumElements () != X.getLocalLength (),
287 std::logic_error,
"Ifpack2::TriDiContainer::applyImpl: The inverse " 288 "operator and X have incompatible dimensions (" <<
289 localMap_->getNodeNumElements () <<
" resp. " 290 << X.getLocalLength () <<
"). Please report this bug to " 291 "the Ifpack2 developers.");
292 TEUCHOS_TEST_FOR_EXCEPTION(
293 localMap_->getNodeNumElements () != Y.getLocalLength (),
294 std::logic_error,
"Ifpack2::TriDiContainer::applyImpl: The inverse " 295 "operator and Y have incompatible dimensions (" <<
296 localMap_->getNodeNumElements () <<
" resp. " 297 << Y.getLocalLength () <<
"). Please report this bug to " 298 "the Ifpack2 developers.");
299 TEUCHOS_TEST_FOR_EXCEPTION(
300 X.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRowsCols()),
301 std::logic_error,
"Ifpack2::TriDiContainer::applyImpl: The input " 302 "multivector X has incompatible dimensions from those of the " 303 "inverse operator (" << X.getLocalLength () <<
" vs. " 304 << (mode == Teuchos::NO_TRANS ? diagBlock_.numRowsCols () : diagBlock_.numRowsCols())
305 <<
"). Please report this bug to the Ifpack2 developers.");
306 TEUCHOS_TEST_FOR_EXCEPTION(
307 Y.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRowsCols()),
308 std::logic_error,
"Ifpack2::TriDiContainer::applyImpl: The output " 309 "multivector Y has incompatible dimensions from those of the " 310 "inverse operator (" << Y.getLocalLength () <<
" vs. " 311 << (mode == Teuchos::NO_TRANS ? diagBlock_.numRowsCols() : diagBlock_.numRowsCols ())
312 <<
"). Please report this bug to the Ifpack2 developers.");
314 typedef Teuchos::ScalarTraits<local_scalar_type> STS;
315 const int numVecs =
static_cast<int> (X.getNumVectors ());
316 if (alpha == STS::zero ()) {
317 if (beta == STS::zero ()) {
321 Y.putScalar (STS::zero ());
328 Teuchos::LAPACK<int, local_scalar_type> lapack;
332 RCP<local_mv_type> Y_tmp;
333 if (beta == STS::zero () ){
335 Y_tmp = rcpFromRef (Y);
338 Y_tmp = rcp (
new local_mv_type (X));
340 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
341 ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
346 (mode == Teuchos::CONJ_TRANS ?
'C' : (mode == Teuchos::TRANS ?
'T' :
'N'));
347 lapack.GTTRS (trans, diagBlock_.numRowsCols(),numVecs,
352 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
353 TEUCHOS_TEST_FOR_EXCEPTION(
354 INFO != 0, std::runtime_error,
"Ifpack2::TriDiContainer::applyImpl: " 355 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) " 356 "failed with INFO = " << INFO <<
" != 0.");
358 if (beta != STS::zero ()) {
359 Y.update (alpha, *Y_tmp, beta);
361 else if (! Y.isConstantStride ()) {
362 Tpetra::deep_copy (Y, *Y_tmp);
367 template<
class MatrixType,
class LocalScalarType>
371 Teuchos::ETransp mode,
372 LocalScalarType alpha,
373 LocalScalarType beta)
const 377 template<
class MatrixType,
class LocalScalarType>
379 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
380 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
381 Teuchos::ETransp mode,
383 scalar_type beta)
const 385 using Teuchos::ArrayView;
400 const size_t numVecs = X.getNumVectors ();
402 TEUCHOS_TEST_FOR_EXCEPTION(
403 ! IsComputed_, std::runtime_error,
"Ifpack2::TriDiContainer::apply: " 404 "You must have called the compute() method before you may call apply(). " 405 "You may call the apply() method as many times as you want after calling " 406 "compute() once, but you must have called compute() at least once.");
407 TEUCHOS_TEST_FOR_EXCEPTION(
408 numVecs != Y.getNumVectors (), std::runtime_error,
409 "Ifpack2::TriDiContainer::apply: X and Y have different numbers of " 410 "vectors. X has " << X.getNumVectors ()
411 <<
", but Y has " << X.getNumVectors () <<
".");
443 X_ = rcp (
new local_mv_type (localMap_, numVecs));
445 RCP<local_mv_type> X_local = X_;
446 TEUCHOS_TEST_FOR_EXCEPTION(
447 X_local->getLocalLength () != numRows_, std::logic_error,
448 "Ifpack2::TriDiContainer::apply: " 449 "X_local has length " << X_local->getLocalLength () <<
", which does " 450 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 451 "the Ifpack2 developers.");
452 ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
453 mvgs.gather (*X_local, X, localRows);
461 Y_ = rcp (
new local_mv_type (localMap_, numVecs));
463 RCP<local_mv_type> Y_local = Y_;
464 TEUCHOS_TEST_FOR_EXCEPTION(
465 Y_local->getLocalLength () != numRows_, std::logic_error,
466 "Ifpack2::TriDiContainer::apply: " 467 "Y_local has length " << X_local->getLocalLength () <<
", which does " 468 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 469 "the Ifpack2 developers.");
470 mvgs.gather (*Y_local, Y, localRows);
474 this->applyImpl (*X_local, *Y_local, mode, as<local_scalar_type> (alpha),
475 as<local_scalar_type> (beta));
479 mvgs.scatter (Y, *Y_local, localRows);
482 template<
class MatrixType,
class LocalScalarType>
484 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
485 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
492 template<
class MatrixType,
class LocalScalarType>
494 weightedApply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
495 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
496 const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
497 Teuchos::ETransp mode,
499 scalar_type beta)
const 501 using Teuchos::ArrayRCP;
502 using Teuchos::ArrayView;
503 using Teuchos::Range1D;
506 using Teuchos::rcp_const_cast;
509 typedef Teuchos::ScalarTraits<scalar_type> STS;
525 global_ordinal_type,
node_type> local_vec_type;
528 const size_t numVecs = X.getNumVectors ();
530 TEUCHOS_TEST_FOR_EXCEPTION(
531 ! IsComputed_, std::runtime_error,
"Ifpack2::TriDiContainer::" 532 "weightedApply: You must have called the compute() method before you may " 533 "call apply(). You may call the apply() method as many times as you want " 534 "after calling compute() once, but you must have called compute() at least " 536 TEUCHOS_TEST_FOR_EXCEPTION(
537 numVecs != Y.getNumVectors (), std::runtime_error,
538 "Ifpack2::TriDiContainer::weightedApply: X and Y have different numbers " 539 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 540 << X.getNumVectors () <<
".");
571 X_ = rcp (
new local_mv_type (localMap_, numVecs));
573 RCP<local_mv_type> X_local = X_;
574 TEUCHOS_TEST_FOR_EXCEPTION(
575 X_local->getLocalLength () != numRows_, std::logic_error,
576 "Ifpack2::TriDiContainer::weightedApply: " 577 "X_local has length " << X_local->getLocalLength () <<
", which does " 578 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 579 "the Ifpack2 developers.");
580 ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
581 mvgs.gather (*X_local, X, localRows);
589 Y_ = rcp (
new local_mv_type (localMap_, numVecs));
591 RCP<local_mv_type> Y_local = Y_;
592 TEUCHOS_TEST_FOR_EXCEPTION(
593 Y_local->getLocalLength () != numRows_, std::logic_error,
594 "Ifpack2::TriDiContainer::weightedApply: " 595 "Y_local has length " << X_local->getLocalLength () <<
", which does " 596 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 597 "the Ifpack2 developers.");
598 mvgs.gather (*Y_local, Y, localRows);
610 RCP<local_vec_type> D_local = rcp (
new local_vec_type (localMap_));
611 TEUCHOS_TEST_FOR_EXCEPTION(
612 D_local->getLocalLength () != numRows_, std::logic_error,
613 "Ifpack2::TriDiContainer::weightedApply: " 614 "D_local has length " << X_local->getLocalLength () <<
", which does " 615 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 616 "the Ifpack2 developers.");
617 mvgs.gather (*D_local, D, localRows);
618 RCP<local_mv_type> X_scaled = rcp (
new local_mv_type (localMap_, numVecs));
619 X_scaled->elementWiseMultiply (STS::one (), *D_local, *X_local, STS::zero ());
626 RCP<local_mv_type> Y_temp;
627 if (beta == STS::zero ()) {
630 Y_temp = rcp (
new local_mv_type (localMap_, numVecs));
634 applyImpl (*X_scaled, *Y_temp, mode, STS::one (), STS::zero ());
640 Y_local->elementWiseMultiply (alpha, *D_local, *Y_temp, beta);
644 mvgs.scatter (Y, *Y_local, localRows);
647 template<
class MatrixType,
class LocalScalarType>
649 weightedApply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
650 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
651 const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
658 template<
class MatrixType,
class LocalScalarType>
661 Teuchos::FancyOStream fos(Teuchos::rcp(&os,
false));
662 fos.setOutputToRootOnly(0);
667 template<
class MatrixType,
class LocalScalarType>
673 template<
class MatrixType,
class LocalScalarType>
676 std::ostringstream oss;
677 oss << Teuchos::Describable::description();
678 if (isInitialized()) {
680 oss <<
"{status = initialized, computed";
683 oss <<
"{status = initialized, not computed";
687 oss <<
"{status = not initialized, not computed";
694 template<
class MatrixType,
class LocalScalarType>
700 template<
class MatrixType,
class LocalScalarType>
704 const Teuchos::EVerbosityLevel verbLevel)
const 707 if(verbLevel==Teuchos::VERB_NONE)
return;
708 os <<
"================================================================================" << endl;
709 os <<
"Ifpack2::TriDiContainer" << endl;
710 os <<
"Number of rows = " << numRows_ << endl;
711 os <<
"isInitialized() = " << IsInitialized_ << endl;
712 os <<
"isComputed() = " << IsComputed_ << endl;
713 os <<
"================================================================================" << endl;
717 template<
class MatrixType,
class LocalScalarType>
721 const Teuchos::EVerbosityLevel )
const 725 template<
class MatrixType,
class LocalScalarType>
728 extract (
const Teuchos::RCP<const row_matrix_type>& globalMatrix)
730 using Teuchos::Array;
731 using Teuchos::ArrayView;
732 using Teuchos::toString;
735 typedef Tpetra::Map<LO, GO, node_type> map_type;
736 const size_t inputMatrixNumRows = globalMatrix->getNodeNumRows ();
740 const int myRank = globalMatrix->getRowMap ()->getComm ()->getRank ();
741 const int numProcs = globalMatrix->getRowMap ()->getComm ()->getSize ();
745 ArrayView<const LO> localRows = this->getLocalRows ();
746 for (
size_t j = 0; j < numRows_; ++j) {
747 TEUCHOS_TEST_FOR_EXCEPTION(
749 static_cast<size_t> (localRows[j]) >= inputMatrixNumRows,
750 std::runtime_error,
"Ifpack2::TriDiContainer::extract: On process " <<
751 myRank <<
" of " << numProcs <<
", localRows[j=" << j <<
"] = " <<
752 localRows[j] <<
", which is out of the valid range of local row indices " 753 "indices [0, " << (inputMatrixNumRows - 1) <<
"] for the input matrix.");
766 const map_type& globalRowMap = * (globalMatrix->getRowMap ());
767 const map_type& globalColMap = * (globalMatrix->getColMap ());
768 const map_type& globalDomMap = * (globalMatrix->getDomainMap ());
770 bool rowIndsValid =
true;
771 bool colIndsValid =
true;
772 Array<LO> localCols (numRows_);
775 Array<LO> invalidLocalRowInds;
776 Array<GO> invalidGlobalColInds;
777 for (
size_t i = 0; i < numRows_; ++i) {
779 const LO ii_local = localRows[i];
784 const GO jj_global = globalRowMap.getGlobalElement (ii_local);
785 if (jj_global == Teuchos::OrdinalTraits<GO>::invalid ()) {
791 rowIndsValid =
false;
792 invalidLocalRowInds.push_back (ii_local);
797 if (globalDomMap.isNodeGlobalElement (jj_global)) {
806 const LO jj_local = globalColMap.getLocalElement (jj_global);
807 if (jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
808 colIndsValid =
false;
809 invalidGlobalColInds.push_back (jj_global);
812 localCols[i] = jj_local;
815 TEUCHOS_TEST_FOR_EXCEPTION(
816 ! rowIndsValid, std::logic_error,
"Ifpack2::TriDiContainer::extract: " 817 "On process " << myRank <<
", at least one row index in the set of local " 818 "row indices given to the constructor is not a valid local row index in " 819 "the input matrix's row Map on this process. This should be impossible " 820 "because the constructor checks for this case. Here is the complete set " 821 "of invalid local row indices: " <<
toString (invalidLocalRowInds) <<
". " 822 "Please report this bug to the Ifpack2 developers.");
823 TEUCHOS_TEST_FOR_EXCEPTION(
824 ! colIndsValid, std::runtime_error,
"Ifpack2::TriDiContainer::extract: " 825 "On process " << myRank <<
", " 826 "At least one row index in the set of row indices given to the constructor " 827 "does not have a corresponding column index in the input matrix's column " 828 "Map. This probably means that the column(s) in question is/are empty on " 829 "this process, which would make the submatrix to extract structurally " 830 "singular. Here is the compete set of invalid global column indices: " 831 <<
toString (invalidGlobalColInds) <<
".");
833 diagBlock_.putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
835 const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries ();
836 Array<scalar_type> val (maxNumEntriesInRow);
837 Array<local_ordinal_type> ind (maxNumEntriesInRow);
840 Teuchos::OrdinalTraits<local_ordinal_type>::invalid ();
841 for (
size_t i = 0; i < numRows_; ++i) {
844 globalMatrix->getLocalRowCopy (localRow, ind (), val (), numEntries);
846 for (
size_t k = 0; k < numEntries; ++k) {
858 if (localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows) {
862 for (
size_t kk = 0; kk < numRows_; ++kk) {
863 if (localRows[kk] == localCol) {
868 diagBlock_ (i, jj) += val[k];
875 template<
class MatrixType,
class LocalScalarType>
878 extract (
const Teuchos::RCP<const row_matrix_type>& )
884 #define IFPACK2_TRIDICONTAINER_INSTANT(S,LO,GO,N) \ 885 template class Ifpack2::TriDiContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >; 887 #endif // IFPACK2_TRIDICONTAINER_HPP MatrixType::node_type node_type
The Node type of the input (global) matrix.
Definition: Ifpack2_TriDiContainer_decl.hpp:137
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_TriDiContainer_decl.hpp:317
Store and solve a local TriDi linear problem.
Definition: Ifpack2_TriDiContainer_decl.hpp:109
Ifpack2::TriDiContainer class declaration.
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_TriDiContainer_decl.hpp:322
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_TriDiContainer_decl.hpp:324
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_TriDiContainer_decl.hpp:128
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_TriDiContainer_decl.hpp:133
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:103
BELOS_DEPRECATED const char * toString(const StatusType status)
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
TriDiContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::ArrayView< const local_ordinal_type > &localRows)
Constructor.
Definition: Ifpack2_TriDiContainer_def.hpp:61
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_TriDiContainer_decl.hpp:135
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85