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" 59 template<
class MatrixType,
class LocalScalarType>
60 DenseContainer<MatrixType, LocalScalarType, true>::
61 DenseContainer (
const Teuchos::RCP<const row_matrix_type>& matrix,
62 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
63 const Teuchos::RCP<const import_type>& importer,
65 scalar_type DampingFactor) :
66 Container<MatrixType> (matrix, partitions, importer, OverlapLevel,
69 scalarOffsets_ (this->numBlocks_)
72 using Teuchos::ArrayView;
76 using Teuchos::toString;
77 typedef typename ArrayView<const local_ordinal_type>::size_type size_type;
78 TEUCHOS_TEST_FOR_EXCEPTION(
79 !matrix->hasColMap(), std::invalid_argument,
"Ifpack2::DenseContainer: " 80 "The constructor's input matrix must have a column Map.");
83 global_ordinal_type totalScalars = 0;
84 for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
86 scalarOffsets_[i] = totalScalars;
87 totalScalars += this->blockRows_[i] * this->blockRows_[i]
88 * this->bcrsBlockSize_ * this->bcrsBlockSize_;
90 scalars_ =
new local_scalar_type[totalScalars];
91 for(
int i = 0; i < this->numBlocks_; i++)
93 int nnodes = this->blockRows_[i];
94 int denseRows = nnodes * this->bcrsBlockSize_;
96 diagBlocks_.emplace_back(Teuchos::View, scalars_ + scalarOffsets_[i], denseRows, denseRows, denseRows);
97 diagBlocks_[i].putScalar(0);
100 ipiv_.resize(this->partitions_.size() * this->bcrsBlockSize_);
102 for(
int i = 0; i < this->numBlocks_; i++)
104 Teuchos::ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
106 const map_type& rowMap = * (matrix->getRowMap ());
107 const size_type numRows = localRows.size ();
108 bool rowIndicesValid =
true;
109 Array<local_ordinal_type> invalidLocalRowIndices;
110 for(size_type j = 0; j < numRows; j++) {
111 if(!rowMap.isNodeLocalElement(localRows[j])) {
112 rowIndicesValid =
false;
113 invalidLocalRowIndices.push_back(localRows[j]);
117 TEUCHOS_TEST_FOR_EXCEPTION(
118 !rowIndicesValid, std::invalid_argument,
"Ifpack2::DenseContainer: " 119 "On process " << rowMap.getComm()->getRank() <<
" of " 120 << rowMap.getComm()->getSize() <<
", in the given set of local row " 121 "indices localRows = " << toString(localRows) <<
", the following " 122 "entries are not valid local row indices on the calling process: " 123 << toString(invalidLocalRowIndices) <<
".");
125 IsInitialized_ =
false;
129 template<
class MatrixType,
class LocalScalarType>
130 DenseContainer<MatrixType, LocalScalarType, true>::
131 DenseContainer (
const Teuchos::RCP<const row_matrix_type>& matrix,
132 const Teuchos::Array<local_ordinal_type>& localRows) :
133 Container<MatrixType>(matrix, localRows),
136 using Teuchos::Array;
137 using Teuchos::ArrayView;
140 using Teuchos::toString;
141 typedef typename ArrayView<const local_ordinal_type>::size_type size_type;
142 TEUCHOS_TEST_FOR_EXCEPTION(
143 !matrix->hasColMap(), std::invalid_argument,
"Ifpack2::DenseContainer: " 144 "The constructor's input matrix must have a column Map.");
145 diagBlocks_.emplace_back(this->blockRows_[0] * this->bcrsBlockSize_,
146 this->blockRows_[0] * this->bcrsBlockSize_);
147 diagBlocks_[0].putScalar(0);
149 ipiv_.resize(this->partitions_.size() * this->bcrsBlockSize_);
151 for(
int i = 0; i < this->numBlocks_; i++)
153 ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
155 const map_type& rowMap = *(matrix->getRowMap());
156 const size_type numRows = localRows.size ();
157 bool rowIndicesValid =
true;
158 Array<local_ordinal_type> invalidLocalRowIndices;
159 for(size_type j = 0; j < numRows; j++)
161 if(!rowMap.isNodeLocalElement(localRows[j]))
163 rowIndicesValid =
false;
164 invalidLocalRowIndices.push_back(localRows[j]);
168 TEUCHOS_TEST_FOR_EXCEPTION(
169 !rowIndicesValid, std::invalid_argument,
"Ifpack2::DenseContainer: " 170 "On process " << rowMap.getComm()->getRank() <<
" of " 171 << rowMap.getComm()->getSize() <<
", in the given set of local row " 172 "indices localRows = " << toString (localRows) <<
", the following " 173 "entries are not valid local row indices on the calling process: " 174 << toString(invalidLocalRowIndices) <<
".");
178 IsInitialized_ =
false;
182 template<
class MatrixType,
class LocalScalarType>
183 DenseContainer<MatrixType, LocalScalarType, true>::~DenseContainer()
189 template<
class MatrixType,
class LocalScalarType>
190 void DenseContainer<MatrixType, LocalScalarType, true>::
191 setParameters (
const Teuchos::ParameterList& )
196 template<
class MatrixType,
class LocalScalarType>
198 DenseContainer<MatrixType, LocalScalarType, true>::
206 IsInitialized_ =
false;
210 for(
int i = 0; i < this->numBlocks_; i++)
211 diagBlocks_[i].putScalar(Teuchos::ScalarTraits<local_scalar_type>::zero());
212 std::fill (ipiv_.begin (), ipiv_.end (), 0);
213 IsInitialized_ =
true;
216 template<
class MatrixType,
class LocalScalarType>
218 DenseContainer<MatrixType, LocalScalarType, true>::
228 if (! this->isInitialized ()) {
239 template<
class MatrixType,
class LocalScalarType>
241 DenseContainer<MatrixType, LocalScalarType, true>::
244 Teuchos::LAPACK<int, local_scalar_type> lapack;
245 for(
int i = 0; i < this->numBlocks_; i++)
248 int* blockIpiv = ipiv_.getRawPtr() + this->partitionIndices_[i] * this->bcrsBlockSize_;
249 lapack.GETRF(diagBlocks_[i].numRows(),
250 diagBlocks_[i].numCols(),
251 diagBlocks_[i].values(),
252 diagBlocks_[i].stride(),
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 INFO < 0, std::logic_error,
"Ifpack2::DenseContainer::factor: " 257 "LAPACK's _GETRF (LU factorization with partial pivoting) was called " 258 "incorrectly. INFO = " << INFO <<
" < 0. " 259 "Please report this bug to the Ifpack2 developers.");
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 INFO > 0, std::runtime_error,
"Ifpack2::DenseContainer::factor: " 265 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the " 266 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") " 267 "(one-based index i) is exactly zero. This probably means that the input " 268 "matrix has a singular diagonal block.");
272 template<
class MatrixType,
class LocalScalarType>
274 DenseContainer<MatrixType, LocalScalarType, true>::
275 applyImplBlockCrs (HostViewLocal& X,
279 Teuchos::ETransp mode,
280 local_scalar_type alpha,
281 local_scalar_type beta)
const 283 using Teuchos::ArrayRCP;
288 using Teuchos::rcpFromRef;
290 typedef Teuchos::ScalarTraits<local_scalar_type> STS;
291 const size_t numRows = X.dimension_0();
292 const size_t numVecs = X.dimension_1();
294 TEUCHOS_TEST_FOR_EXCEPTION(
295 static_cast<size_t> (X.dimension_0 ()) != static_cast<size_t> (diagBlocks_[blockIndex].numRows ()),
296 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: X and Y have " 297 "different number of rows than block matrix (" << X.dimension_0() <<
" resp. " 298 << diagBlocks_[blockIndex].numRows() <<
"). Please report this bug to " 299 "the Ifpack2 developers.");
301 if (alpha == STS::zero ()) {
302 if (beta == STS::zero ()) {
306 for(
size_t i = 0; i < numRows; i++)
307 for(
size_t j = 0; j < numVecs; j++)
308 Y(i, j) = STS::zero();
311 for(
size_t i = 0; i < numRows; i++)
312 for(
size_t j = 0; j < numVecs; j++)
317 Teuchos::LAPACK<int, local_scalar_type> lapack;
321 Ptr<HostViewLocal> Y_tmp;
322 bool deleteYT =
false;
323 if (beta == STS::zero () ){
324 Kokkos::deep_copy(Y, X);
328 Y_tmp = ptr (
new HostViewLocal (
"", X.dimension_0(), X.dimension_1()));
329 Kokkos::deep_copy(*Y_tmp, X);
332 local_scalar_type*
const Y_ptr = (local_scalar_type*) Y_tmp->ptr_on_device();
335 (mode == Teuchos::CONJ_TRANS ?
'C' : (mode == Teuchos::TRANS ?
'T' :
'N'));
336 int* blockIpiv = (
int*) ipiv_.getRawPtr()
337 + this->partitionIndices_[blockIndex] * this->bcrsBlockSize_;
339 diagBlocks_[blockIndex].numRows (),
341 diagBlocks_[blockIndex].values (),
342 diagBlocks_[blockIndex].stride (),
346 TEUCHOS_TEST_FOR_EXCEPTION(
347 INFO != 0, std::runtime_error,
"Ifpack2::DenseContainer::applyImpl: " 348 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) " 349 "failed with INFO = " << INFO <<
" != 0.");
351 if (beta != STS::zero ()) {
352 for(
size_t i = 0; i < Y.dimension_0(); i++)
354 for(
size_t j = 0; j < Y.dimension_1(); j++)
357 Y(i, j) += alpha * (*Y_tmp)(i, j);
366 template<
class MatrixType,
class LocalScalarType>
368 DenseContainer<MatrixType, LocalScalarType, true>::
369 applyImpl (HostViewLocal& X,
373 Teuchos::ETransp mode,
374 local_scalar_type alpha,
375 local_scalar_type beta)
const 377 using Teuchos::ArrayRCP;
382 using Teuchos::rcpFromRef;
384 TEUCHOS_TEST_FOR_EXCEPTION(
385 X.dimension_0 () != Y.dimension_0 (),
386 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: X and Y have " 387 "incompatible dimensions (" << X.dimension_0 () <<
" resp. " 388 << Y.dimension_0 () <<
"). Please report this bug to " 389 "the Ifpack2 developers.");
391 TEUCHOS_TEST_FOR_EXCEPTION(
392 X.dimension_1 () != Y.dimension_1(),
393 std::logic_error,
"Ifpack2::DenseContainer::applyImpl: X and Y have " 394 "incompatible numbers of vectors (" << X.dimension_1 () <<
" resp. " 395 << Y.dimension_1 () <<
"). Please report this bug to " 396 "the Ifpack2 developers.");
398 if(this->hasBlockCrs_) {
399 applyImplBlockCrs(X,Y,blockIndex,stride,mode,alpha,beta);
403 typedef Teuchos::ScalarTraits<local_scalar_type> STS;
404 size_t numVecs = X.dimension_1();
405 if(alpha == STS::zero()) {
406 if(beta == STS::zero()) {
410 for(
size_t i = 0; i < Y.dimension_0(); i++)
412 for(
size_t j = 0; j < Y.dimension_1(); j++)
413 Y(i, j) = STS::zero();
417 for(
size_t i = 0; i < Y.dimension_0(); i++)
419 for(
size_t j = 0; j < Y.dimension_1(); j++)
424 Teuchos::LAPACK<int, local_scalar_type> lapack;
428 Ptr<HostViewLocal> Y_tmp;
429 bool deleteYT =
false;
430 if (beta == STS::zero () ){
431 Kokkos::deep_copy (Y, X);
435 Y_tmp = ptr (
new HostViewLocal (
"", Y.dimension_0(), Y.dimension_1()));
438 local_scalar_type* Y_ptr = (local_scalar_type*) Y_tmp->ptr_on_device();
440 int* blockIpiv = (
int*) ipiv_.getRawPtr() + this->partitionIndices_[blockIndex] * this->bcrsBlockSize_;
442 (mode == Teuchos::CONJ_TRANS ?
'C' : (mode == Teuchos::TRANS ?
'T' :
'N'));
444 diagBlocks_[blockIndex].numRows (),
446 diagBlocks_[blockIndex].values (),
447 diagBlocks_[blockIndex].stride (),
451 TEUCHOS_TEST_FOR_EXCEPTION(
452 INFO != 0, std::runtime_error,
"Ifpack2::DenseContainer::applyImpl: " 453 "LAPACK's _GETRS (solve using LU factorization with partial pivoting) " 454 "failed with INFO = " << INFO <<
" != 0.");
456 if (beta != STS::zero ()) {
457 for(
size_t i = 0; i < Y.dimension_0(); i++)
459 for(
size_t j = 0; j < Y.dimension_1(); j++)
460 Y(i, j) = Y(i, j) * (local_impl_scalar_type) beta + (local_impl_scalar_type) alpha * (*Y_tmp)(i, j);
468 template<
class MatrixType,
class LocalScalarType>
470 DenseContainer<MatrixType, LocalScalarType, true>::
471 applyBlockCrs (HostView& XIn,
475 Teuchos::ETransp mode,
477 scalar_type beta)
const 479 using Teuchos::ArrayView;
480 using Teuchos::ArrayRCP;
485 const size_t numRows = this->blockRows_[blockIndex];
493 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
494 TEUCHOS_TEST_FOR_EXCEPTION(
495 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 496 "compute() method before you may call this method. You may call " 497 "apply() as many times as you want after calling compute() once, " 498 "but you must have called compute() at least once first.");
499 const size_t numVecs = XIn.dimension_1 ();
500 TEUCHOS_TEST_FOR_EXCEPTION(
501 numVecs != YIn.dimension_1 (), std::runtime_error,
502 prefix <<
"X and Y have different numbers of vectors (columns). X has " 503 << XIn.dimension_1 () <<
", but Y has " << YIn.dimension_1 () <<
".");
534 if(X_local.size() == 0)
538 for(
int i = 0; i < this->numBlocks_; i++)
540 X_local.emplace_back(
"", this->blockRows_[i] * this->bcrsBlockSize_, numVecs);
542 for(
int i = 0; i < this->numBlocks_; i++)
544 Y_local.emplace_back(
"", this->blockRows_[i] * this->bcrsBlockSize_, numVecs);
547 HostViewLocal& XOut = X_local[blockIndex];
548 HostViewLocal& YOut = Y_local[blockIndex];
550 ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
552 for (
size_t j = 0; j < numVecs; ++j) {
553 for (
size_t i = 0; i < numRows; ++i) {
554 const size_t i_perm = localRows[i];
555 for (
int k = 0; k < this->bcrsBlockSize_; ++k)
556 XOut(i*this->bcrsBlockSize_+k, j) = XIn(i_perm*this->bcrsBlockSize_+k, j);
566 for (
size_t j = 0; j < numVecs; ++j) {
567 for (
size_t i = 0; i < numRows; ++i) {
568 const size_t i_perm = localRows[i];
569 for (
int k = 0; k < this->bcrsBlockSize_; ++k)
570 YOut(i*this->bcrsBlockSize_+k, j) = YIn(i_perm*this->bcrsBlockSize_+k, j);
576 this->applyImpl (XOut, YOut, blockIndex, stride, mode, as<local_scalar_type>(alpha),
577 as<local_scalar_type>(beta));
581 for(
size_t j = 0; j < numVecs; ++j) {
582 for(
size_t i = 0; i < numRows; ++i) {
583 const size_t i_perm = localRows[i];
584 for(
int k = 0; k < this->bcrsBlockSize_; ++k)
585 YIn(i_perm*this->bcrsBlockSize_+k, j) = YOut(i*this->bcrsBlockSize_+k, j);
590 template<
class MatrixType,
class LocalScalarType>
592 DenseContainer<MatrixType, LocalScalarType, true>::
597 Teuchos::ETransp mode,
599 scalar_type beta)
const 601 using Teuchos::ArrayView;
607 if(this->hasBlockCrs_) {
608 applyBlockCrs(X,Y,blockIndex,stride,mode,alpha,beta);
611 const size_t numVecs = X.dimension_1();
619 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
620 TEUCHOS_TEST_FOR_EXCEPTION(
621 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 622 "compute() method before you may call this method. You may call " 623 "apply() as many times as you want after calling compute() once, " 624 "but you must have called compute() at least once first.");
625 TEUCHOS_TEST_FOR_EXCEPTION(
626 X.dimension_1 () != Y.dimension_1 (), std::runtime_error,
627 prefix <<
"X and Y have different numbers of vectors (columns). X has " 628 << X.dimension_1 () <<
", but Y has " << Y.dimension_1 () <<
".");
659 if(X_local.size() == 0)
663 for(
int i = 0; i < this->numBlocks_; i++)
665 X_local.emplace_back(
"", this->blockRows_[i], numVecs);
667 for(
int i = 0; i < this->numBlocks_; i++)
669 Y_local.emplace_back(
"", this->blockRows_[i], numVecs);
673 const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
675 Details::MultiVectorLocalGatherScatter<mv_type, local_mv_type> mvgs;
676 mvgs.gatherViewToView (X_local[blockIndex], X, localRows);
683 mvgs.gatherViewToView (Y_local[blockIndex], Y, localRows);
687 this->applyImpl (X_local[blockIndex], Y_local[blockIndex], blockIndex, stride, mode,
688 as<local_scalar_type>(alpha), as<local_scalar_type>(beta));
692 mvgs.scatterViewToView (Y, Y_local[blockIndex], localRows);
695 template<
class MatrixType,
class LocalScalarType>
696 void DenseContainer<MatrixType, LocalScalarType, true>::
697 weightedApply (HostView& X,
702 Teuchos::ETransp mode,
704 scalar_type beta)
const 706 using Teuchos::ArrayRCP;
707 using Teuchos::ArrayView;
708 using Teuchos::Range1D;
713 using Teuchos::rcp_const_cast;
715 typedef Teuchos::ScalarTraits<scalar_type> STS;
724 const char prefix[] =
"Ifpack2::DenseContainer::weightedApply: ";
725 TEUCHOS_TEST_FOR_EXCEPTION(
726 ! IsComputed_, std::runtime_error, prefix <<
"You must have called the " 727 "compute() method before you may call this method. You may call " 728 "weightedApply() as many times as you want after calling compute() once, " 729 "but you must have called compute() at least once first.");
731 const size_t numVecs = X.dimension_1();
733 TEUCHOS_TEST_FOR_EXCEPTION(
734 X.dimension_1() != Y.dimension_1(), std::runtime_error,
735 prefix <<
"X and Y have different numbers of vectors (columns). X has " 736 << X.dimension_1() <<
", but Y has " << Y.dimension_1() <<
".");
742 const size_t numRows = this->blockRows_[blockIndex];
768 if(X_local.size() == 0)
772 for(
int i = 0; i < this->numBlocks_; i++)
774 X_local.emplace_back(
"", this->blockRows_[i], numVecs);
776 for(
int i = 0; i < this->numBlocks_; i++)
778 Y_local.emplace_back(
"", this->blockRows_[i], numVecs);
782 ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
784 Details::MultiVectorLocalGatherScatter<mv_type, local_mv_type> mvgs;
785 mvgs.gatherViewToView (X_local[blockIndex], X, localRows);
791 mvgs.gatherViewToView (Y_local[blockIndex], Y, localRows);
803 HostViewLocal D_local(
"", numRows, 1);
804 mvgs.gatherViewToView (D_local, D, localRows);
805 HostViewLocal X_scaled(
"", numRows, numVecs);
806 for(
size_t j = 0; j < numVecs; j++)
807 for(
size_t i = 0; i < numRows; i++)
808 X_scaled(i, j) = X_local[blockIndex](i, j) * D_local(i, 0);
815 Ptr<HostViewLocal> Y_temp;
816 bool deleteYT =
false;
817 if(beta == STS::zero())
819 Y_temp = ptr(&Y_local[blockIndex]);
821 Y_temp = ptr(
new HostViewLocal(
"", numRows, numVecs));
826 this->applyImpl (X_scaled, *Y_temp, blockIndex, stride, mode, STS::one(), STS::zero());
832 for(
size_t j = 0; j < numVecs; j++)
833 for(
size_t i = 0; i < numRows; i++)
834 Y_local[blockIndex](i, j) = Y_local[blockIndex](i, j) * (local_impl_scalar_type) beta + (local_impl_scalar_type) alpha * (*Y_temp)(i, j) * D_local(i, 0);
841 mvgs.scatterViewToView (Y, Y_local[blockIndex], localRows);
844 template<
class MatrixType,
class LocalScalarType>
846 DenseContainer<MatrixType, LocalScalarType, true>::
847 print (std::ostream& os)
const 849 Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
850 fos.setOutputToRootOnly (0);
851 this->describe (fos);
855 template<
class MatrixType,
class LocalScalarType>
857 DenseContainer<MatrixType, LocalScalarType, true>::
860 std::ostringstream oss;
861 oss <<
"Ifpack::DenseContainer: ";
862 if (isInitialized()) {
864 oss <<
"{status = initialized, computed";
867 oss <<
"{status = initialized, not computed";
871 oss <<
"{status = not initialized, not computed";
878 template<
class MatrixType,
class LocalScalarType>
880 DenseContainer<MatrixType, LocalScalarType, true>::
881 describe (Teuchos::FancyOStream& os,
882 const Teuchos::EVerbosityLevel verbLevel)
const 885 if(verbLevel==Teuchos::VERB_NONE)
return;
886 os <<
"================================================================================" << endl;
887 os <<
"Ifpack2::DenseContainer" << endl;
888 for(
int i = 0; i < this->numBlocks_; i++)
890 os <<
"Block " << i <<
" number of rows = " << this->blockRows_[i] << endl;
892 os <<
"isInitialized() = " << IsInitialized_ << endl;
893 os <<
"isComputed() = " << IsComputed_ << endl;
894 os <<
"================================================================================" << endl;
899 template<
class MatrixType,
class LocalScalarType>
901 DenseContainer<MatrixType, LocalScalarType, true>::
904 using Teuchos::Array;
905 using Teuchos::ArrayView;
906 using Teuchos::toString;
907 auto& A = this->inputMatrix_;
908 const size_t inputMatrixNumRows = A->getNodeNumRows();
912 const int myRank = A->getRowMap ()->getComm ()->getRank ();
913 const int numProcs = A->getRowMap ()->getComm ()->getSize ();
917 for(
int i = 0; i < this->numBlocks_; ++i) {
918 ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
919 for(local_ordinal_type j = 0; j < this->blockRows_[i]; ++j) {
920 TEUCHOS_TEST_FOR_EXCEPTION(
922 static_cast<size_t>(localRows[j]) >= inputMatrixNumRows,
923 std::runtime_error,
"Ifpack2::DenseContainer::extract: On process " <<
924 myRank <<
" of " << numProcs <<
", localRows[j=" << j <<
"] = " <<
925 localRows[j] <<
", which is out of the valid range of local row indices " 926 "indices [0, " << (inputMatrixNumRows - 1) <<
"] for the input matrix.");
940 auto globalRowMap = A->getRowMap ();
941 auto globalColMap = A->getColMap ();
942 auto globalDomMap = A->getDomainMap ();
944 for(
int blockIndex = 0; blockIndex < this->numBlocks_; blockIndex++)
946 const local_ordinal_type numRows_ = this->blockRows_[blockIndex];
947 Teuchos::ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
948 bool rowIndsValid =
true;
949 bool colIndsValid =
true;
950 Array<local_ordinal_type> localCols(numRows_);
953 Array<local_ordinal_type> invalidLocalRowInds;
954 Array<global_ordinal_type> invalidGlobalColInds;
955 for (local_ordinal_type i = 0; i < numRows_; i++)
958 const local_ordinal_type ii_local = localRows[i];
963 const global_ordinal_type jj_global = globalRowMap->getGlobalElement(ii_local);
964 if(jj_global == Teuchos::OrdinalTraits<global_ordinal_type>::invalid())
971 rowIndsValid =
false;
972 invalidLocalRowInds.push_back(ii_local);
977 if(globalDomMap->isNodeGlobalElement(jj_global))
987 const local_ordinal_type jj_local = globalColMap->getLocalElement(jj_global);
988 if(jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid())
990 colIndsValid =
false;
991 invalidGlobalColInds.push_back(jj_global);
994 localCols[i] = jj_local;
997 TEUCHOS_TEST_FOR_EXCEPTION(
998 !rowIndsValid, std::logic_error,
"Ifpack2::DenseContainer::extract: " 999 "On process " << myRank <<
", at least one row index in the set of local " 1000 "row indices given to the constructor is not a valid local row index in " 1001 "the input matrix's row Map on this process. This should be impossible " 1002 "because the constructor checks for this case. Here is the complete set " 1003 "of invalid local row indices: " << toString(invalidLocalRowInds) <<
". " 1004 "Please report this bug to the Ifpack2 developers.");
1005 TEUCHOS_TEST_FOR_EXCEPTION(
1006 !colIndsValid, std::runtime_error,
"Ifpack2::DenseContainer::extract: " 1007 "On process " << myRank <<
", " 1008 "At least one row index in the set of row indices given to the constructor " 1009 "does not have a corresponding column index in the input matrix's column " 1010 "Map. This probably means that the column(s) in question is/are empty on " 1011 "this process, which would make the submatrix to extract structurally " 1012 "singular. Here is the compete set of invalid global column indices: " 1013 << toString(invalidGlobalColInds) <<
".");
1015 diagBlocks_[blockIndex].putScalar(Teuchos::ScalarTraits<local_scalar_type>::zero());
1017 const size_t maxNumEntriesInRow = A->getNodeMaxNumRowEntries();
1018 Array<local_ordinal_type> ind(maxNumEntriesInRow);
1020 const local_ordinal_type INVALID = Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
1022 Array<scalar_type> val(maxNumEntriesInRow * this->bcrsBlockSize_ * this->bcrsBlockSize_);
1023 for(local_ordinal_type i = 0; i < numRows_; i++)
1025 const local_ordinal_type localRow = localRows[i];
1027 A->getLocalRowCopy(localRow, ind(), val(), numEntries);
1029 for(
size_t k = 0; k < numEntries; k++)
1031 const local_ordinal_type localCol = ind[k];
1041 if(localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows)
1045 local_ordinal_type jj = INVALID;
1046 for(local_ordinal_type kk = 0; kk < numRows_; kk++)
1048 if(localRows[kk] == localCol)
1054 for(local_ordinal_type c = 0; c < this->bcrsBlockSize_; c++)
1056 for(local_ordinal_type r = 0; r < this->bcrsBlockSize_; r++)
1057 diagBlocks_[blockIndex](this->bcrsBlockSize_ * i + r,
1058 this->bcrsBlockSize_ * jj + c)
1059 = val[k * (this->bcrsBlockSize_ * this->bcrsBlockSize_)
1060 + (r + this->bcrsBlockSize_ * c)];
1070 template<
class MatrixType,
class LocalScalarType>
1072 DenseContainer<MatrixType, LocalScalarType, true>::
1075 using Teuchos::Array;
1076 using Teuchos::ArrayView;
1077 using Teuchos::toString;
1078 auto& A = *this->inputMatrix_;
1079 const size_t inputMatrixNumRows = A.getNodeNumRows();
1083 const int myRank = A.getRowMap ()->getComm ()->getRank ();
1084 const int numProcs = A.getRowMap ()->getComm ()->getSize ();
1086 for(
int blockIndex = 0; blockIndex < this->numBlocks_; blockIndex++)
1088 local_ordinal_type numRows_ = this->blockRows_[blockIndex];
1090 if(this->hasBlockCrs_)
1098 ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
1099 for(local_ordinal_type j = 0; j < numRows_; j++)
1101 TEUCHOS_TEST_FOR_EXCEPTION(
1103 static_cast<size_t> (localRows[j]) >= inputMatrixNumRows,
1104 std::runtime_error,
"Ifpack2::DenseContainer::extract: On process " <<
1105 myRank <<
" of " << numProcs <<
", localRows[j=" << j <<
"] = " <<
1106 localRows[j] <<
", which is out of the valid range of local row indices " 1107 "indices [0, " << (inputMatrixNumRows - 1) <<
"] for the input matrix.");
1120 const map_type& globalRowMap = * (A.getRowMap ());
1121 const map_type& globalColMap = * (A.getColMap ());
1122 const map_type& globalDomMap = * (A.getDomainMap ());
1124 bool rowIndsValid =
true;
1125 bool colIndsValid =
true;
1126 Array<local_ordinal_type> localCols(numRows_);
1129 Array<local_ordinal_type> invalidLocalRowInds;
1130 Array<global_ordinal_type> invalidGlobalColInds;
1131 for(local_ordinal_type i = 0; i < numRows_; i++)
1134 const local_ordinal_type ii_local = localRows[i];
1139 const global_ordinal_type jj_global = globalRowMap.getGlobalElement(ii_local);
1140 if(jj_global == Teuchos::OrdinalTraits<global_ordinal_type>::invalid())
1147 rowIndsValid =
false;
1148 invalidLocalRowInds.push_back(ii_local);
1153 if(globalDomMap.isNodeGlobalElement(jj_global))
1163 const local_ordinal_type jj_local = globalColMap.getLocalElement(jj_global);
1164 if(jj_local == Teuchos::OrdinalTraits<local_ordinal_type>::invalid())
1166 colIndsValid =
false;
1167 invalidGlobalColInds.push_back(jj_global);
1170 localCols[i] = jj_local;
1173 TEUCHOS_TEST_FOR_EXCEPTION(
1174 !rowIndsValid, std::logic_error,
"Ifpack2::DenseContainer::extract: " 1175 "On process " << myRank <<
", at least one row index in the set of local " 1176 "row indices given to the constructor is not a valid local row index in " 1177 "the input matrix's row Map on this process. This should be impossible " 1178 "because the constructor checks for this case. Here is the complete set " 1179 "of invalid local row indices: " << toString(invalidLocalRowInds) <<
". " 1180 "Please report this bug to the Ifpack2 developers.");
1181 TEUCHOS_TEST_FOR_EXCEPTION(
1182 !colIndsValid, std::runtime_error,
"Ifpack2::DenseContainer::extract: " 1183 "On process " << myRank <<
", " 1184 "At least one row index in the set of row indices given to the constructor " 1185 "does not have a corresponding column index in the input matrix's column " 1186 "Map. This probably means that the column(s) in question is/are empty on " 1187 "this process, which would make the submatrix to extract structurally " 1188 "singular. Here is the compete set of invalid global column indices: " 1189 << toString(invalidGlobalColInds) <<
".");
1191 diagBlocks_[blockIndex].putScalar(Teuchos::ScalarTraits<local_scalar_type>::zero());
1193 const size_t maxNumEntriesInRow = A.getNodeMaxNumRowEntries();
1194 Array<local_ordinal_type> ind(maxNumEntriesInRow);
1196 const local_ordinal_type INVALID = Teuchos::OrdinalTraits<local_ordinal_type>::invalid();
1198 Array<scalar_type> val(maxNumEntriesInRow);
1199 for (local_ordinal_type i = 0; i < numRows_; i++)
1201 const local_ordinal_type localRow = localRows[i];
1203 A.getLocalRowCopy(localRow, ind(), val(), numEntries);
1204 for (
size_t k = 0; k < numEntries; ++k)
1206 const local_ordinal_type localCol = ind[k];
1216 if(localCol >= 0 && static_cast<size_t> (localCol) < inputMatrixNumRows)
1220 local_ordinal_type jj = INVALID;
1221 for(local_ordinal_type kk = 0; kk < numRows_; kk++)
1223 if(localRows[kk] == localCol)
1227 diagBlocks_[blockIndex](i, jj) += val[k];
1234 template<
class MatrixType,
class LocalScalarType>
1235 void DenseContainer<MatrixType, LocalScalarType, true>::clearBlocks()
1237 std::vector<Teuchos::SerialDenseMatrix<int, local_scalar_type>> empty1;
1238 std::swap(diagBlocks_, empty1);
1239 Teuchos::Array<int> empty2;
1240 Teuchos::swap(ipiv_, empty2);
1241 std::vector<HostViewLocal> empty3;
1242 std::swap(X_local, empty3);
1243 std::vector<HostViewLocal> empty4;
1244 std::swap(Y_local, empty4);
1245 Container<MatrixType>::clearBlocks();
1248 template<
class MatrixType,
class LocalScalarType>
1249 std::string DenseContainer<MatrixType, LocalScalarType, true>::getName()
1254 template<
class MatrixType,
class LocalScalarType>
1255 DenseContainer<MatrixType, LocalScalarType, false>::
1256 DenseContainer (
const Teuchos::RCP<const row_matrix_type>& matrix,
1257 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
1258 const Teuchos::RCP<const import_type>& importer,
1260 scalar_type DampingFactor) :
1261 Container<MatrixType> (matrix, partitions, importer, OverlapLevel,
1264 TEUCHOS_TEST_FOR_EXCEPTION
1265 (
true, std::logic_error,
"Ifpack2::DenseContainer: Not implemented for " 1266 "LocalScalarType = " << Teuchos::TypeNameTraits<LocalScalarType>::name ()
1270 template<
class MatrixType,
class LocalScalarType>
1271 DenseContainer<MatrixType, LocalScalarType, false>::
1272 DenseContainer (
const Teuchos::RCP<const row_matrix_type>& matrix,
1273 const Teuchos::Array<local_ordinal_type>& localRows) :
1274 Container<MatrixType>(matrix, localRows)
1276 TEUCHOS_TEST_FOR_EXCEPTION
1277 (
true, std::logic_error,
"Ifpack2::DenseContainer: Not implemented for " 1278 "LocalScalarType = " << Teuchos::TypeNameTraits<LocalScalarType>::name ()
1282 template<
class MatrixType,
class LocalScalarType>
1283 DenseContainer<MatrixType, LocalScalarType, false>::~DenseContainer() {}
1285 template<
class MatrixType,
class LocalScalarType>
1286 void DenseContainer<MatrixType, LocalScalarType, false>::
1287 setParameters (
const Teuchos::ParameterList& ) {}
1289 template<
class MatrixType,
class LocalScalarType>
1290 void DenseContainer<MatrixType, LocalScalarType, false>::initialize() {}
1292 template<
class MatrixType,
class LocalScalarType>
1293 void DenseContainer<MatrixType, LocalScalarType, false>::compute() {}
1295 template<
class MatrixType,
class LocalScalarType>
1296 void DenseContainer<MatrixType, LocalScalarType, false>::factor() {}
1298 template<
class MatrixType,
class LocalScalarType>
1299 void DenseContainer<MatrixType, LocalScalarType, false>::
1300 applyImplBlockCrs (HostViewLocal& X,
1304 Teuchos::ETransp mode,
1305 local_scalar_type alpha,
1306 local_scalar_type beta)
const {}
1308 template<
class MatrixType,
class LocalScalarType>
1309 void DenseContainer<MatrixType, LocalScalarType, false>::
1310 applyImpl (HostViewLocal& X,
1314 Teuchos::ETransp mode,
1315 local_scalar_type alpha,
1316 local_scalar_type beta)
const {}
1318 template<
class MatrixType,
class LocalScalarType>
1319 void DenseContainer<MatrixType, LocalScalarType, false>::
1320 applyBlockCrs (HostView& XIn,
1324 Teuchos::ETransp mode,
1326 scalar_type beta)
const {}
1328 template<
class MatrixType,
class LocalScalarType>
1329 void DenseContainer<MatrixType, LocalScalarType, false>::
1334 Teuchos::ETransp mode,
1336 scalar_type beta)
const {}
1338 template<
class MatrixType,
class LocalScalarType>
1339 void DenseContainer<MatrixType, LocalScalarType, false>::
1340 weightedApply (HostView& X,
1345 Teuchos::ETransp mode,
1347 scalar_type beta)
const {}
1349 template<
class MatrixType,
class LocalScalarType>
1350 std::ostream& DenseContainer<MatrixType, LocalScalarType, false>::
1351 print (std::ostream& os)
const 1356 template<
class MatrixType,
class LocalScalarType>
1357 std::string DenseContainer<MatrixType, LocalScalarType, false>::
1358 description ()
const 1363 template<
class MatrixType,
class LocalScalarType>
1364 void DenseContainer<MatrixType, LocalScalarType, false>::
1365 describe (Teuchos::FancyOStream& os,
1366 const Teuchos::EVerbosityLevel verbLevel)
const {}
1368 template<
class MatrixType,
class LocalScalarType>
1369 void DenseContainer<MatrixType, LocalScalarType, false>::
1370 extractBlockCrs () {}
1372 template<
class MatrixType,
class LocalScalarType>
1373 void DenseContainer<MatrixType, LocalScalarType, false>::
1376 template<
class MatrixType,
class LocalScalarType>
1377 void DenseContainer<MatrixType, LocalScalarType, false>::clearBlocks() {}
1379 template<
class MatrixType,
class LocalScalarType>
1380 std::string DenseContainer<MatrixType, LocalScalarType, false>::getName()
1391 #define IFPACK2_DENSECONTAINER_INSTANT(S,LO,GO,N) \ 1392 template class Ifpack2::DenseContainer< Tpetra::RowMatrix<S, LO, GO, N>, S >; 1394 #endif // IFPACK2_DENSECONTAINER_DEF_HPP Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72