43 #ifndef IFPACK2_DENSECONTAINER_DEF_HPP 44 #define IFPACK2_DENSECONTAINER_DEF_HPP 46 #include "Tpetra_CrsMatrix.hpp" 47 #include "Teuchos_LAPACK.hpp" 48 #include "Tpetra_Experimental_BlockMultiVector.hpp" 52 # include "Teuchos_DefaultMpiComm.hpp" 54 # include "Teuchos_DefaultSerialComm.hpp" 60 template<
class MatrixType,
class LocalScalarType>
63 const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
64 Container<MatrixType> (matrix, localRows),
65 hasBlockCrsMatrix_(false),
66 numRows_ (localRows.size ()),
67 diagBlock_ (numRows_, numRows_),
71 using Teuchos::ArrayView;
74 using Teuchos::toString;
75 typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
76 typedef typename ArrayView<const local_ordinal_type>::size_type size_type;
77 TEUCHOS_TEST_FOR_EXCEPTION(
78 ! matrix->hasColMap (), std::invalid_argument,
"Ifpack2::DenseContainer: " 79 "The constructor's input matrix must have a column Map.");
82 Teuchos::RCP<const block_crs_matrix_type> global_bcrs =
83 Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type> (matrix);
84 if (!global_bcrs.is_null ()) {
85 hasBlockCrsMatrix_ =
true;
88 blockSize_ = global_bcrs->getBlockSize();
91 diagBlock_.shape(numRows_*blockSize_,numRows_*blockSize_);
92 ipiv_.resize(numRows_*blockSize_);
96 const map_type& rowMap = * (matrix->getRowMap ());
97 const size_type numRows = localRows.size ();
98 bool rowIndicesValid =
true;
99 Array<local_ordinal_type> invalidLocalRowIndices;
100 for (size_type i = 0; i < numRows; ++i) {
101 if (! rowMap.isNodeLocalElement (localRows[i])) {
102 rowIndicesValid =
false;
103 invalidLocalRowIndices.push_back (localRows[i]);
107 TEUCHOS_TEST_FOR_EXCEPTION(
108 ! rowIndicesValid, std::invalid_argument,
"Ifpack2::DenseContainer: " 109 "On process " << rowMap.getComm ()->getRank () <<
" of " 110 << rowMap.getComm ()->getSize () <<
", in the given set of local row " 111 "indices localRows = " <<
toString (localRows) <<
", the following " 112 "entries are not valid local row indices on the calling process: " 113 <<
toString (invalidLocalRowIndices) <<
".");
116 RCP<const Teuchos::Comm<int> > localComm =
117 rcp (
new Teuchos::MpiComm<int> (MPI_COMM_SELF));
119 RCP<const Teuchos::Comm<int> > localComm =
120 rcp (
new Teuchos::SerialComm<int> ());
126 localMap_ = rcp (
new map_type (numRows_, indexBase, localComm));
129 template<
class MatrixType,
class LocalScalarType>
132 const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
133 Container<MatrixType> (matrix, localRows)
135 TEUCHOS_TEST_FOR_EXCEPTION
136 (
true, std::logic_error,
"Ifpack2::DenseContainer: Not implemented for " 137 "LocalScalarType = " << Teuchos::TypeNameTraits<LocalScalarType>::name ()
141 template<
class MatrixType,
class LocalScalarType>
145 template<
class MatrixType,
class LocalScalarType>
150 template<
class MatrixType,
class LocalScalarType>
157 template<
class MatrixType,
class LocalScalarType>
165 template<
class MatrixType,
class LocalScalarType>
175 IsInitialized_ =
false;
179 diagBlock_.putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
180 std::fill (ipiv_.begin (), ipiv_.end (), 0);
182 IsInitialized_ =
true;
185 template<
class MatrixType,
class LocalScalarType>
191 template<
class MatrixType,
class LocalScalarType>
203 if (! this->isInitialized ()) {
208 extract (this->getMatrix ());
214 template<
class MatrixType,
class LocalScalarType>
220 template<
class MatrixType,
class LocalScalarType>
225 Teuchos::LAPACK<int, local_scalar_type> lapack;
227 lapack.GETRF (diagBlock_.numRows (), diagBlock_.numCols (),
228 diagBlock_.values (), diagBlock_.stride (),
229 ipiv_.getRawPtr (), &INFO);
231 TEUCHOS_TEST_FOR_EXCEPTION(
232 INFO < 0, std::logic_error,
"Ifpack2::DenseContainer::factor: " 233 "LAPACK's _GETRF (LU factorization with partial pivoting) was called " 234 "incorrectly. INFO = " << INFO <<
" < 0. " 235 "Please report this bug to the Ifpack2 developers.");
239 TEUCHOS_TEST_FOR_EXCEPTION(
240 INFO > 0, std::runtime_error,
"Ifpack2::DenseContainer::factor: " 241 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the " 242 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") " 243 "(one-based index i) is exactly zero. This probably means that the input " 244 "matrix has a singular diagonal block.");
247 template<
class MatrixType,
class LocalScalarType>
253 template<
class MatrixType,
class LocalScalarType>
258 Teuchos::ETransp mode,
259 LocalScalarType alpha,
260 LocalScalarType beta)
const 262 using Teuchos::ArrayRCP;
265 using Teuchos::rcpFromRef;
267 TEUCHOS_TEST_FOR_EXCEPTION(
268 X.getLocalLength () != Y.getLocalLength (),
269 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: X and Y have " 270 "incompatible dimensions (" << X.getLocalLength () <<
" resp. " 271 << Y.getLocalLength () <<
"). Please report this bug to " 272 "the Ifpack2 developers.");
273 TEUCHOS_TEST_FOR_EXCEPTION(
274 localMap_->getNodeNumElements()*blockSize_ != X.getLocalLength (),
275 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The inverse " 276 "operator and X have incompatible dimensions (" <<
277 localMap_->getNodeNumElements()*blockSize_ <<
" resp. " 278 << X.getLocalLength () <<
"). Please report this bug to " 279 "the Ifpack2 developers.");
280 TEUCHOS_TEST_FOR_EXCEPTION(
281 localMap_->getNodeNumElements()*blockSize_ != Y.getLocalLength (),
282 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The inverse " 283 "operator and Y have incompatible dimensions (" <<
284 localMap_->getNodeNumElements()*blockSize_ <<
" resp. " 285 << Y.getLocalLength () <<
"). Please report this bug to " 286 "the Ifpack2 developers.");
287 TEUCHOS_TEST_FOR_EXCEPTION(
288 X.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRows ()),
289 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The input " 290 "multivector X has incompatible dimensions from those of the " 291 "inverse operator (" << X.getLocalLength () <<
" vs. " 292 << (mode == Teuchos::NO_TRANS ? diagBlock_.numCols () : diagBlock_.numRows ())
293 <<
"). Please report this bug to the Ifpack2 developers.");
294 TEUCHOS_TEST_FOR_EXCEPTION(
295 X.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRows ()),
296 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The output " 297 "multivector Y has incompatible dimensions from those of the " 298 "inverse operator (" << Y.getLocalLength () <<
" vs. " 299 << (mode == Teuchos::NO_TRANS ? diagBlock_.numRows () : diagBlock_.numCols ())
300 <<
"). Please report this bug to the Ifpack2 developers.");
302 typedef Teuchos::ScalarTraits<local_scalar_type> STS;
303 const int numVecs =
static_cast<int> (X.getNumVectors ());
304 if (alpha == STS::zero ()) {
305 if (beta == STS::zero ()) {
309 Y.putScalar (STS::zero ());
316 Teuchos::LAPACK<int, local_scalar_type> lapack;
320 RCP<local_mv_type> Y_tmp;
321 if (beta == STS::zero () ){
322 Tpetra::deep_copy (Y, X);
323 Y_tmp = rcpFromRef (Y);
326 Y_tmp = rcp (
new local_mv_type (X, Teuchos::Copy));
328 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
329 ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
334 (mode == Teuchos::CONJ_TRANS ?
'C' : (mode == Teuchos::TRANS ?
'T' :
'N'));
335 lapack.GETRS (trans, diagBlock_.numRows (), numVecs,
336 diagBlock_.values (), diagBlock_.stride (),
337 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
338 TEUCHOS_TEST_FOR_EXCEPTION(
339 INFO != 0, std::runtime_error,
"Ifpack2::DenseContainer::applyImpl: " 340 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) " 341 "failed with INFO = " << INFO <<
" != 0.");
343 if (beta != STS::zero ()) {
344 Y.update (alpha, *Y_tmp, beta);
346 else if (! Y.isConstantStride ()) {
347 Tpetra::deep_copy (Y, *Y_tmp);
352 template<
class MatrixType,
class LocalScalarType>
357 Teuchos::ETransp mode,
358 LocalScalarType alpha,
359 LocalScalarType beta)
const 361 using Teuchos::ArrayRCP;
364 using Teuchos::rcpFromRef;
366 if(hasBlockCrsMatrix_) {
367 applyImplBlockCrs(X,Y,mode,alpha,beta);
371 TEUCHOS_TEST_FOR_EXCEPTION(
372 X.getLocalLength () != Y.getLocalLength (),
373 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: X and Y have " 374 "incompatible dimensions (" << X.getLocalLength () <<
" resp. " 375 << Y.getLocalLength () <<
"). Please report this bug to " 376 "the Ifpack2 developers.");
377 TEUCHOS_TEST_FOR_EXCEPTION(
378 localMap_->getNodeNumElements () != X.getLocalLength (),
379 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The inverse " 380 "operator and X have incompatible dimensions (" <<
381 localMap_->getNodeNumElements () <<
" resp. " 382 << X.getLocalLength () <<
"). Please report this bug to " 383 "the Ifpack2 developers.");
384 TEUCHOS_TEST_FOR_EXCEPTION(
385 localMap_->getNodeNumElements () != Y.getLocalLength (),
386 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The inverse " 387 "operator and Y have incompatible dimensions (" <<
388 localMap_->getNodeNumElements () <<
" resp. " 389 << Y.getLocalLength () <<
"). Please report this bug to " 390 "the Ifpack2 developers.");
391 TEUCHOS_TEST_FOR_EXCEPTION(
392 X.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRows ()),
393 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The input " 394 "multivector X has incompatible dimensions from those of the " 395 "inverse operator (" << X.getLocalLength () <<
" vs. " 396 << (mode == Teuchos::NO_TRANS ? diagBlock_.numCols () : diagBlock_.numRows ())
397 <<
"). Please report this bug to the Ifpack2 developers.");
398 TEUCHOS_TEST_FOR_EXCEPTION(
399 X.getLocalLength () !=
static_cast<size_t> (diagBlock_.numRows ()),
400 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: The output " 401 "multivector Y has incompatible dimensions from those of the " 402 "inverse operator (" << Y.getLocalLength () <<
" vs. " 403 << (mode == Teuchos::NO_TRANS ? diagBlock_.numRows () : diagBlock_.numCols ())
404 <<
"). Please report this bug to the Ifpack2 developers.");
406 typedef Teuchos::ScalarTraits<local_scalar_type> STS;
407 const int numVecs =
static_cast<int> (X.getNumVectors ());
408 if (alpha == STS::zero ()) {
409 if (beta == STS::zero ()) {
413 Y.putScalar (STS::zero ());
420 Teuchos::LAPACK<int, local_scalar_type> lapack;
424 RCP<local_mv_type> Y_tmp;
425 if (beta == STS::zero () ){
426 Tpetra::deep_copy (Y, X);
427 Y_tmp = rcpFromRef (Y);
430 Y_tmp = rcp (
new local_mv_type (X, Teuchos::Copy));
432 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
433 ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
438 (mode == Teuchos::CONJ_TRANS ?
'C' : (mode == Teuchos::TRANS ?
'T' :
'N'));
439 lapack.GETRS (trans, diagBlock_.numRows (), numVecs,
440 diagBlock_.values (), diagBlock_.stride (),
441 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
442 TEUCHOS_TEST_FOR_EXCEPTION(
443 INFO != 0, std::runtime_error,
"Ifpack2::DenseContainer::applyImpl: " 444 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) " 445 "failed with INFO = " << INFO <<
" != 0.");
447 if (beta != STS::zero ()) {
448 Y.update (alpha, *Y_tmp, beta);
450 else if (! Y.isConstantStride ()) {
451 Tpetra::deep_copy (Y, *Y_tmp);
456 template<
class MatrixType,
class LocalScalarType>
463 LocalScalarType )
const 466 template<
class MatrixType,
class LocalScalarType>
473 LocalScalarType )
const 477 template<
class MatrixType,
class LocalScalarType>
480 applyBlockCrs (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
481 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
482 Teuchos::ETransp mode,
484 scalar_type beta)
const 486 using Teuchos::ArrayView;
487 using Teuchos::ArrayRCP;
498 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
499 TEUCHOS_TEST_FOR_EXCEPTION(
500 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 501 "compute() method before you may call this method. You may call " 502 "apply() as many times as you want after calling compute() once, " 503 "but you must have called compute() at least once first.");
504 const size_t numVecs = X.getNumVectors ();
505 TEUCHOS_TEST_FOR_EXCEPTION(
506 numVecs != Y.getNumVectors (), std::runtime_error,
507 prefix <<
"X and Y have different numbers of vectors (columns). X has " 508 << X.getNumVectors () <<
", but Y has " << X.getNumVectors () <<
".");
539 if (blockX_.is_null ()) {
540 blockX_ = rcp (
new local_block_mv_type (*localMap_, blockSize_, numVecs));
542 local_mv_type X_local = blockX_->getMultiVectorView();
543 TEUCHOS_TEST_FOR_EXCEPTION(
544 X_local.getLocalLength () != numRows_*blockSize_, std::logic_error,
545 "Ifpack2::DenseContainer::apply: " 546 "X_local has length " << X_local.getLocalLength () <<
", which does " 547 "not match numRows_*blockSize_ = " << numRows_*blockSize_ <<
". Please report this bug to " 548 "the Ifpack2 developers.");
549 ArrayView<const local_ordinal_type> localRows = this->getLocalRows();
552 for (
size_t j = 0; j < numVecs; ++j) {
553 ArrayRCP<const scalar_type> X_in_j = X.getData (j);
554 ArrayRCP<local_scalar_type> X_out_j = X_local.getDataNonConst (j);
556 for (
size_t i = 0; i < numRows_; ++i) {
557 const size_t i_perm = localRows[i];
558 for (
int k = 0; k < blockSize_; k++)
559 X_out_j[i*blockSize_+k] = X_in_j[i_perm*blockSize_+k];
568 if (blockY_.is_null ()) {
569 blockY_ = rcp (
new local_block_mv_type (*localMap_, blockSize_, numVecs));
571 local_mv_type Y_local = blockY_->getMultiVectorView();
572 TEUCHOS_TEST_FOR_EXCEPTION(
573 Y_local.getLocalLength () != numRows_*blockSize_, std::logic_error,
574 "Ifpack2::DenseContainer::apply: " 575 "Y_local has length " << X_local.getLocalLength () <<
", which does " 576 "not match numRows_ = " << numRows_*blockSize_ <<
". Please report this bug to " 577 "the Ifpack2 developers.");
580 for (
size_t j = 0; j < numVecs; ++j) {
581 ArrayRCP<const scalar_type> Y_in_j = Y.getData (j);
582 ArrayRCP<local_scalar_type> Y_out_j = Y_local.getDataNonConst (j);
584 for (
size_t i = 0; i < numRows_; ++i) {
585 const size_t i_perm = localRows[i];
586 for (
int k = 0; k < blockSize_; k++)
587 Y_out_j[i*blockSize_+k] = Y_in_j[i_perm*blockSize_+k];
593 this->applyImpl (X_local, Y_local, mode, as<local_scalar_type> (alpha),
594 as<local_scalar_type> (beta));
598 for (
size_t j = 0; j < numVecs; ++j) {
599 ArrayRCP<local_scalar_type> Y_in_j = Y.getDataNonConst (j);
600 ArrayRCP<const scalar_type> Y_out_j = Y_local.getData (j);
602 for (
size_t i = 0; i < numRows_; ++i) {
603 const size_t i_perm = localRows[i];
604 for (
int k = 0; k < blockSize_; k++)
605 Y_in_j[i_perm*blockSize_+k] = Y_out_j[i*blockSize_+k];
611 template<
class MatrixType,
class LocalScalarType>
614 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
615 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
616 Teuchos::ETransp mode,
618 scalar_type beta)
const 620 using Teuchos::ArrayView;
626 if(hasBlockCrsMatrix_) {
627 applyBlockCrs(X,Y,mode,alpha,beta);
638 typedef Tpetra::MultiVector<scalar_type,
641 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
642 TEUCHOS_TEST_FOR_EXCEPTION(
643 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 644 "compute() method before you may call this method. You may call " 645 "apply() as many times as you want after calling compute() once, " 646 "but you must have called compute() at least once first.");
647 const size_t numVecs = X.getNumVectors ();
648 TEUCHOS_TEST_FOR_EXCEPTION(
649 numVecs != Y.getNumVectors (), std::runtime_error,
650 prefix <<
"X and Y have different numbers of vectors (columns). X has " 651 << X.getNumVectors () <<
", but Y has " << X.getNumVectors () <<
".");
683 X_ = rcp (
new local_mv_type (localMap_, numVecs));
685 RCP<local_mv_type> X_local = X_;
686 TEUCHOS_TEST_FOR_EXCEPTION(
687 X_local->getLocalLength () != numRows_, std::logic_error,
688 "Ifpack2::DenseContainer::apply: " 689 "X_local has length " << X_local->getLocalLength () <<
", which does " 690 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 691 "the Ifpack2 developers.");
692 ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
695 mvgs.gather (*X_local, X, localRows);
703 Y_ = rcp (
new local_mv_type (localMap_, numVecs));
705 RCP<local_mv_type> Y_local = Y_;
706 TEUCHOS_TEST_FOR_EXCEPTION(
707 Y_local->getLocalLength () != numRows_, std::logic_error,
708 "Ifpack2::DenseContainer::apply: " 709 "Y_local has length " << X_local->getLocalLength () <<
", which does " 710 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 711 "the Ifpack2 developers.");
712 mvgs.gather (*Y_local, Y, localRows);
716 this->applyImpl (*X_local, *Y_local, mode, as<local_scalar_type> (alpha),
717 as<local_scalar_type> (beta));
721 mvgs.scatter (Y, *Y_local, localRows);
724 template<
class MatrixType,
class LocalScalarType>
727 applyBlockCrs (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
728 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
734 template<
class MatrixType,
class LocalScalarType>
737 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
738 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
744 template<
class MatrixType,
class LocalScalarType>
746 weightedApply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
747 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
748 const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
749 Teuchos::ETransp mode,
751 scalar_type beta)
const 753 using Teuchos::ArrayRCP;
754 using Teuchos::ArrayView;
755 using Teuchos::Range1D;
758 using Teuchos::rcp_const_cast;
760 typedef Teuchos::ScalarTraits<scalar_type> STS;
770 typedef Tpetra::MultiVector<scalar_type,
776 global_ordinal_type,
node_type> local_vec_type;
778 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
779 TEUCHOS_TEST_FOR_EXCEPTION(
780 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 781 "compute() method before you may call this method. You may call " 782 "weightedApply() as many times as you want after calling compute() once, " 783 "but you must have called compute() at least once first.");
784 const size_t numVecs = X.getNumVectors ();
785 TEUCHOS_TEST_FOR_EXCEPTION(
786 numVecs != Y.getNumVectors (), std::runtime_error,
787 prefix <<
"X and Y have different numbers of vectors (columns). X has " 788 << X.getNumVectors () <<
", but Y has " << X.getNumVectors () <<
".");
819 X_ = rcp (
new local_mv_type (localMap_, numVecs));
821 RCP<local_mv_type> X_local = X_;
822 TEUCHOS_TEST_FOR_EXCEPTION(
823 X_local->getLocalLength () != numRows_, std::logic_error,
824 "Ifpack2::DenseContainer::weightedApply: " 825 "X_local has length " << X_local->getLocalLength () <<
", which does " 826 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 827 "the Ifpack2 developers.");
828 ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
831 mvgs.gather (*X_local, X, localRows);
839 Y_ = rcp (
new local_mv_type (localMap_, numVecs));
841 RCP<local_mv_type> Y_local = Y_;
842 TEUCHOS_TEST_FOR_EXCEPTION(
843 Y_local->getLocalLength () != numRows_, std::logic_error,
844 "Ifpack2::DenseContainer::weightedApply: " 845 "Y_local has length " << X_local->getLocalLength () <<
", which does " 846 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 847 "the Ifpack2 developers.");
848 mvgs.gather (*Y_local, Y, localRows);
860 local_vec_type D_local (localMap_);
861 TEUCHOS_TEST_FOR_EXCEPTION(
862 D_local.getLocalLength () != numRows_, std::logic_error,
863 "Ifpack2::DenseContainer::weightedApply: " 864 "D_local has length " << D_local.getLocalLength () <<
", which does " 865 "not match numRows_ = " << numRows_ <<
". Please report this bug to " 866 "the Ifpack2 developers.");
867 mvgs.gather (D_local, D, localRows);
868 local_mv_type X_scaled (localMap_, numVecs);
869 X_scaled.elementWiseMultiply (STS::one (), D_local, *X_local, STS::zero ());
876 RCP<local_mv_type> Y_temp;
877 if (beta == STS::zero ()) {
880 Y_temp = rcp (
new local_mv_type (localMap_, numVecs));
884 this->applyImpl (X_scaled, *Y_temp, mode, STS::one (), STS::zero ());
890 Y_local->elementWiseMultiply (alpha, D_local, *Y_temp, beta);
894 mvgs.scatter (Y, *Y_local, localRows);
897 template<
class MatrixType,
class LocalScalarType>
900 weightedApply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
901 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
902 const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& ,
909 template<
class MatrixType,
class LocalScalarType>
914 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
915 fos.setOutputToRootOnly (0);
916 this->describe (fos);
920 template<
class MatrixType,
class LocalScalarType>
928 template<
class MatrixType,
class LocalScalarType>
933 std::ostringstream oss;
934 oss << Teuchos::Describable::description();
935 if (isInitialized()) {
937 oss <<
"{status = initialized, computed";
940 oss <<
"{status = initialized, not computed";
944 oss <<
"{status = not initialized, not computed";
951 template<
class MatrixType,
class LocalScalarType>
959 template<
class MatrixType,
class LocalScalarType>
963 const Teuchos::EVerbosityLevel verbLevel)
const 966 if(verbLevel==Teuchos::VERB_NONE)
return;
967 os <<
"================================================================================" << endl;
968 os <<
"Ifpack2::DenseContainer" << endl;
969 os <<
"Number of rows = " << numRows_ << endl;
970 os <<
"isInitialized() = " << IsInitialized_ << endl;
971 os <<
"isComputed() = " << IsComputed_ << endl;
972 os <<
"================================================================================" << endl;
976 template<
class MatrixType,
class LocalScalarType>
980 const Teuchos::EVerbosityLevel verbLevel)
const 984 template<
class MatrixType,
class LocalScalarType>
987 extractBlockCrs (
const Teuchos::RCP<const block_crs_matrix_type>& globalMatrix)
989 using Teuchos::Array;
990 using Teuchos::ArrayView;
991 using Teuchos::toString;
994 typedef Tpetra::Map<LO, GO, node_type> map_type;
995 const size_t inputMatrixNumRows = globalMatrix->getNodeNumRows ();
999 const int myRank = globalMatrix->getRowMap ()->getComm ()->getRank ();
1000 const int numProcs = globalMatrix->getRowMap ()->getComm ()->getSize ();
1004 ArrayView<const LO> localRows = this->getLocalRows ();
1005 for (
size_t j = 0; j < numRows_; ++j) {
1006 TEUCHOS_TEST_FOR_EXCEPTION(
1008 static_cast<size_t> (localRows[j]) >= inputMatrixNumRows,
1009 std::runtime_error,
"Ifpack2::DenseContainer::extract: On process " <<
1010 myRank <<
" of " << numProcs <<
", localRows[j=" << j <<
"] = " <<
1011 localRows[j] <<
", which is out of the valid range of local row indices " 1012 "indices [0, " << (inputMatrixNumRows - 1) <<
"] for the input matrix.");
1025 Teuchos::RCP<const map_type> globalRowMap = globalMatrix->getRowMap ();
1026 Teuchos::RCP<const map_type> globalColMap = globalMatrix->getColMap ();
1027 Teuchos::RCP<const map_type> globalDomMap = globalMatrix->getDomainMap ();
1029 bool rowIndsValid =
true;
1030 bool colIndsValid =
true;
1031 Array<LO> localCols (numRows_);
1034 Array<LO> invalidLocalRowInds;
1035 Array<GO> invalidGlobalColInds;
1036 for (
size_t i = 0; i < numRows_; ++i) {
1038 const LO ii_local = localRows[i];
1043 const GO jj_global = globalRowMap->getGlobalElement (ii_local);
1044 if (jj_global == Teuchos::OrdinalTraits<GO>::invalid ()) {
1050 rowIndsValid =
false;
1051 invalidLocalRowInds.push_back (ii_local);
1056 if (globalDomMap->isNodeGlobalElement (jj_global)) {
1065 const LO jj_local = globalColMap->getLocalElement (jj_global);
1066 if (jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
1067 colIndsValid =
false;
1068 invalidGlobalColInds.push_back (jj_global);
1071 localCols[i] = jj_local;
1074 TEUCHOS_TEST_FOR_EXCEPTION(
1075 ! rowIndsValid, std::logic_error,
"Ifpack2::DenseContainer::extract: " 1076 "On process " << myRank <<
", at least one row index in the set of local " 1077 "row indices given to the constructor is not a valid local row index in " 1078 "the input matrix's row Map on this process. This should be impossible " 1079 "because the constructor checks for this case. Here is the complete set " 1080 "of invalid local row indices: " <<
toString (invalidLocalRowInds) <<
". " 1081 "Please report this bug to the Ifpack2 developers.");
1082 TEUCHOS_TEST_FOR_EXCEPTION(
1083 ! colIndsValid, std::runtime_error,
"Ifpack2::DenseContainer::extract: " 1084 "On process " << myRank <<
", " 1085 "At least one row index in the set of row indices given to the constructor " 1086 "does not have a corresponding column index in the input matrix's column " 1087 "Map. This probably means that the column(s) in question is/are empty on " 1088 "this process, which would make the submatrix to extract structurally " 1089 "singular. Here is the compete set of invalid global column indices: " 1090 <<
toString (invalidGlobalColInds) <<
".");
1092 diagBlock_.putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
1094 const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries ();
1095 Array<local_ordinal_type> ind (maxNumEntriesInRow);
1098 Teuchos::OrdinalTraits<local_ordinal_type>::invalid ();
1100 Array<scalar_type> val (maxNumEntriesInRow*blockSize_*blockSize_);
1101 for (
size_t i = 0; i < numRows_; ++i) {
1104 globalMatrix->getLocalRowCopy (localRow, ind (), val (), numEntries);
1106 for (
size_t k = 0; k < numEntries; ++k) {
1118 if (localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows) {
1122 for (
size_t kk = 0; kk < numRows_; ++kk) {
1123 if (localRows[kk] == localCol) {
1127 if (jj != INVALID) {
1134 diagBlock_(blockSize_*i+r,blockSize_*jj+c) = val[k*blockSize_*blockSize_+r+c*blockSize_];
1145 template<
class MatrixType,
class LocalScalarType>
1148 extract (
const Teuchos::RCP<const row_matrix_type>& globalMatrix)
1150 using Teuchos::Array;
1151 using Teuchos::ArrayView;
1152 using Teuchos::toString;
1155 typedef Tpetra::Map<LO, GO, node_type> map_type;
1156 const size_t inputMatrixNumRows = globalMatrix->getNodeNumRows ();
1160 const int myRank = globalMatrix->getRowMap ()->getComm ()->getRank ();
1161 const int numProcs = globalMatrix->getRowMap ()->getComm ()->getSize ();
1164 if(hasBlockCrsMatrix_) {
1165 extractBlockCrs(Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (globalMatrix));
1171 ArrayView<const LO> localRows = this->getLocalRows ();
1172 for (
size_t j = 0; j < numRows_; ++j) {
1173 TEUCHOS_TEST_FOR_EXCEPTION(
1175 static_cast<size_t> (localRows[j]) >= inputMatrixNumRows,
1176 std::runtime_error,
"Ifpack2::DenseContainer::extract: On process " <<
1177 myRank <<
" of " << numProcs <<
", localRows[j=" << j <<
"] = " <<
1178 localRows[j] <<
", which is out of the valid range of local row indices " 1179 "indices [0, " << (inputMatrixNumRows - 1) <<
"] for the input matrix.");
1192 const map_type& globalRowMap = * (globalMatrix->getRowMap ());
1193 const map_type& globalColMap = * (globalMatrix->getColMap ());
1194 const map_type& globalDomMap = * (globalMatrix->getDomainMap ());
1196 bool rowIndsValid =
true;
1197 bool colIndsValid =
true;
1198 Array<LO> localCols (numRows_);
1201 Array<LO> invalidLocalRowInds;
1202 Array<GO> invalidGlobalColInds;
1203 for (
size_t i = 0; i < numRows_; ++i) {
1205 const LO ii_local = localRows[i];
1210 const GO jj_global = globalRowMap.getGlobalElement (ii_local);
1211 if (jj_global == Teuchos::OrdinalTraits<GO>::invalid ()) {
1217 rowIndsValid =
false;
1218 invalidLocalRowInds.push_back (ii_local);
1223 if (globalDomMap.isNodeGlobalElement (jj_global)) {
1232 const LO jj_local = globalColMap.getLocalElement (jj_global);
1233 if (jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid ()) {
1234 colIndsValid =
false;
1235 invalidGlobalColInds.push_back (jj_global);
1238 localCols[i] = jj_local;
1241 TEUCHOS_TEST_FOR_EXCEPTION(
1242 ! rowIndsValid, std::logic_error,
"Ifpack2::DenseContainer::extract: " 1243 "On process " << myRank <<
", at least one row index in the set of local " 1244 "row indices given to the constructor is not a valid local row index in " 1245 "the input matrix's row Map on this process. This should be impossible " 1246 "because the constructor checks for this case. Here is the complete set " 1247 "of invalid local row indices: " <<
toString (invalidLocalRowInds) <<
". " 1248 "Please report this bug to the Ifpack2 developers.");
1249 TEUCHOS_TEST_FOR_EXCEPTION(
1250 ! colIndsValid, std::runtime_error,
"Ifpack2::DenseContainer::extract: " 1251 "On process " << myRank <<
", " 1252 "At least one row index in the set of row indices given to the constructor " 1253 "does not have a corresponding column index in the input matrix's column " 1254 "Map. This probably means that the column(s) in question is/are empty on " 1255 "this process, which would make the submatrix to extract structurally " 1256 "singular. Here is the compete set of invalid global column indices: " 1257 <<
toString (invalidGlobalColInds) <<
".");
1259 diagBlock_.putScalar (Teuchos::ScalarTraits<local_scalar_type>::zero ());
1261 const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries ();
1262 Array<local_ordinal_type> ind (maxNumEntriesInRow);
1265 Teuchos::OrdinalTraits<local_ordinal_type>::invalid ();
1267 Array<scalar_type> val (maxNumEntriesInRow);
1268 for (
size_t i = 0; i < numRows_; ++i) {
1271 globalMatrix->getLocalRowCopy (localRow, ind (), val (), numEntries);
1273 for (
size_t k = 0; k < numEntries; ++k) {
1285 if (localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows) {
1289 for (
size_t kk = 0; kk < numRows_; ++kk) {
1290 if (localRows[kk] == localCol) {
1294 if (jj != INVALID) {
1295 diagBlock_ (i, jj) += val[k];
1302 template<
class MatrixType,
class LocalScalarType>
1308 template<
class MatrixType,
class LocalScalarType>
1311 extract (
const Teuchos::RCP<const row_matrix_type>& )
1320 #define IFPACK2_DENSECONTAINER_INSTANT(S,LO,GO,N) \ 1321 template class Ifpack2::DenseContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >; 1323 #endif // IFPACK2_SPARSECONTAINER_HPP DenseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::ArrayView< const local_ordinal_type > &localRows)
Constructor.
Definition: Ifpack2_DenseContainer_def.hpp:62
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_DenseContainer_decl.hpp:365
Store and solve a local dense linear problem.
Definition: Ifpack2_DenseContainer_decl.hpp:109
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_DenseContainer_decl.hpp:358
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input (global) matrix.
Definition: Ifpack2_DenseContainer_decl.hpp:135
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_DenseContainer_decl.hpp:363
LocalScalarType local_scalar_type
The second template parameter of this class.
Definition: Ifpack2_DenseContainer_decl.hpp:128
MatrixType::node_type node_type
The Node type of the input (global) matrix.
Definition: Ifpack2_DenseContainer_decl.hpp:137
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input (global) matrix.
Definition: Ifpack2_DenseContainer_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
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85