43 #ifndef IFPACK2_EXPERIMENTAL_CRSRBILUK_DEF_HPP 44 #define IFPACK2_EXPERIMENTAL_CRSRBILUK_DEF_HPP 46 #include "Tpetra_BlockMultiVector.hpp" 47 #include "Tpetra_BlockView.hpp" 48 #include "Ifpack2_OverlappingRowMatrix.hpp" 49 #include "Ifpack2_LocalFilter.hpp" 51 #include "Ifpack2_RILUK.hpp" 54 #define IFPACK2_RBILUK_INITIAL_NOKK 56 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 57 #include "KokkosBatched_Gemm_Decl.hpp" 58 #include "KokkosBatched_Gemm_Serial_Impl.hpp" 59 #include "KokkosBatched_Util.hpp" 66 template<
class MatrixType>
70 A_block_(
Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(Matrix_in))
73 template<
class MatrixType>
80 template<
class MatrixType>
84 template<
class MatrixType>
94 if (A.getRawPtr () != A_block_.getRawPtr ())
96 this->isAllocated_ =
false;
97 this->isInitialized_ =
false;
98 this->isComputed_ =
false;
99 this->Graph_ = Teuchos::null;
100 L_block_ = Teuchos::null;
101 U_block_ = Teuchos::null;
102 D_block_ = Teuchos::null;
109 template<
class MatrixType>
110 const typename RBILUK<MatrixType>::block_crs_matrix_type&
113 TEUCHOS_TEST_FOR_EXCEPTION(
114 L_block_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor " 115 "is null. Please call initialize() (and preferably also compute()) " 116 "before calling this method. If the input matrix has not yet been set, " 117 "you must first call setMatrix() with a nonnull input matrix before you " 118 "may call initialize() or compute().");
123 template<
class MatrixType>
124 const typename RBILUK<MatrixType>::block_crs_matrix_type&
127 TEUCHOS_TEST_FOR_EXCEPTION(
128 D_block_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor " 129 "(of diagonal entries) is null. Please call initialize() (and " 130 "preferably also compute()) before calling this method. If the input " 131 "matrix has not yet been set, you must first call setMatrix() with a " 132 "nonnull input matrix before you may call initialize() or compute().");
137 template<
class MatrixType>
138 const typename RBILUK<MatrixType>::block_crs_matrix_type&
141 TEUCHOS_TEST_FOR_EXCEPTION(
142 U_block_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U factor " 143 "is null. Please call initialize() (and preferably also compute()) " 144 "before calling this method. If the input matrix has not yet been set, " 145 "you must first call setMatrix() with a nonnull input matrix before you " 146 "may call initialize() or compute().");
150 template<
class MatrixType>
156 if (! this->isAllocated_) {
168 L_block_ = rcp(
new block_crs_matrix_type(*this->Graph_->getL_Graph (), blockSize_) );
169 U_block_ = rcp(
new block_crs_matrix_type(*this->Graph_->getU_Graph (), blockSize_) );
170 D_block_ = rcp(
new block_crs_matrix_type(*(Ifpack2::Details::computeDiagonalGraph(*(this->Graph_->getOverlapGraph()))),
172 L_block_->setAllToScalar (STM::zero ());
173 U_block_->setAllToScalar (STM::zero ());
174 D_block_->setAllToScalar (STM::zero ());
177 this->isAllocated_ =
true;
180 template<
class MatrixType>
181 Teuchos::RCP<const typename RBILUK<MatrixType>::block_crs_matrix_type>
186 template<
class MatrixType>
191 using Teuchos::rcp_dynamic_cast;
192 const char prefix[] =
"Ifpack2::Experimental::RBILUK::initialize: ";
203 if (A_block_.is_null ()) {
208 RCP<const LocalFilter<row_matrix_type> > filteredA =
210 TEUCHOS_TEST_FOR_EXCEPTION
211 (filteredA.is_null (), std::runtime_error, prefix <<
212 "Cannot cast to filtered matrix.");
213 RCP<const OverlappingRowMatrix<row_matrix_type> > overlappedA =
215 if (! overlappedA.is_null ()) {
216 A_block_ = rcp_dynamic_cast<
const block_crs_matrix_type>(overlappedA->getUnderlyingMatrix());
219 A_block_ = rcp_dynamic_cast<
const block_crs_matrix_type>(filteredA->getUnderlyingMatrix());
223 TEUCHOS_TEST_FOR_EXCEPTION
224 (A_block_.is_null (), std::runtime_error, prefix <<
"The matrix (A_block_, " 225 "the BlockCrsMatrix) is null. Please call setMatrix() with a nonnull " 226 "input before calling this method.");
227 TEUCHOS_TEST_FOR_EXCEPTION
228 (! A_block_->isFillComplete (), std::runtime_error, prefix <<
"The matrix " 229 "(A_block_, the BlockCrsMatrix) is not fill complete. You may not invoke " 230 "initialize() or compute() with this matrix until the matrix is fill " 231 "complete. Note: BlockCrsMatrix is fill complete if and only if its " 232 "underlying graph is fill complete.");
234 blockSize_ = A_block_->getBlockSize();
236 Teuchos::Time timer (
"RBILUK::initialize");
237 double startTime = timer.wallTime();
239 Teuchos::TimeMonitor timeMon (timer);
248 this->isInitialized_ =
false;
249 this->isAllocated_ =
false;
250 this->isComputed_ =
false;
251 this->Graph_ = Teuchos::null;
257 RCP<const crs_graph_type> matrixCrsGraph = Teuchos::rcpFromRef(A_block_->getCrsGraph() );
259 this->LevelOfFill_, 0));
261 this->Graph_->initialize ();
262 allocate_L_and_U_blocks ();
263 #ifdef IFPACK2_RBILUK_INITIAL 264 initAllValues (*A_block_);
268 this->isInitialized_ =
true;
269 this->numInitialize_ += 1;
270 this->initializeTime_ += (timer.wallTime() - startTime);
274 template<
class MatrixType>
280 typedef Tpetra::Map<LO,GO,node_type> map_type;
282 LO NumIn = 0, NumL = 0, NumU = 0;
283 bool DiagFound =
false;
284 size_t NumNonzeroDiags = 0;
285 size_t MaxNumEntries = A.getNodeMaxNumRowEntries();
286 LO blockMatSize = blockSize_*blockSize_;
291 Teuchos::ArrayView<const GO> rowGIDs = A.getRowMap()->getNodeElementList();
292 Teuchos::ArrayView<const GO> colGIDs = A.getColMap()->getNodeElementList();
293 bool gidsAreConsistentlyOrdered=
true;
294 GO indexOfInconsistentGID=0;
295 for (GO i=0; i<rowGIDs.size(); ++i) {
296 if (rowGIDs[i] != colGIDs[i]) {
297 gidsAreConsistentlyOrdered=
false;
298 indexOfInconsistentGID=i;
302 TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==
false, std::runtime_error,
303 "The ordering of the local GIDs in the row and column maps is not the same" 304 << std::endl <<
"at index " << indexOfInconsistentGID
305 <<
". Consistency is required, as all calculations are done with" 306 << std::endl <<
"local indexing.");
310 Teuchos::Array<LO> LI(MaxNumEntries);
311 Teuchos::Array<LO> UI(MaxNumEntries);
312 Teuchos::Array<scalar_type> LV(MaxNumEntries*blockMatSize);
313 Teuchos::Array<scalar_type> UV(MaxNumEntries*blockMatSize);
315 Teuchos::Array<scalar_type> diagValues(blockMatSize);
317 L_block_->setAllToScalar (STM::zero ());
318 U_block_->setAllToScalar (STM::zero ());
319 D_block_->setAllToScalar (STM::zero ());
336 RCP<const map_type> rowMap = L_block_->getRowMap ();
348 using inds_type =
typename row_matrix_type::local_inds_host_view_type;
349 using vals_type =
typename row_matrix_type::values_host_view_type;
350 for (
size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
351 LO local_row = myRow;
357 A.getLocalRowView(local_row, InI, InV);
358 NumIn = (LO)InI.size();
366 for (LO j = 0; j < NumIn; ++j) {
368 const LO blockOffset = blockMatSize*j;
370 if (k == local_row) {
373 for (LO jj = 0; jj < blockMatSize; ++jj)
374 diagValues[jj] = this->Rthresh_ * InV[blockOffset+jj] + IFPACK2_SGN(InV[blockOffset+jj]) * this->Athresh_;
375 D_block_->replaceLocalValues(local_row, &InI[j], diagValues.getRawPtr(), 1);
378 TEUCHOS_TEST_FOR_EXCEPTION(
379 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current " 380 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is " 381 "probably an artifact of the undocumented assumptions of the " 382 "original implementation (likely copied and pasted from Ifpack). " 383 "Nevertheless, the code I found here insisted on this being an error " 384 "state, so I will throw an exception here.");
386 else if (k < local_row) {
388 const LO LBlockOffset = NumL*blockMatSize;
389 for (LO jj = 0; jj < blockMatSize; ++jj)
390 LV[LBlockOffset+jj] = InV[blockOffset+jj];
393 else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
395 const LO UBlockOffset = NumU*blockMatSize;
396 for (LO jj = 0; jj < blockMatSize; ++jj)
397 UV[UBlockOffset+jj] = InV[blockOffset+jj];
408 for (LO jj = 0; jj < blockSize_; ++jj)
409 diagValues[jj*(blockSize_+1)] = this->Athresh_;
410 D_block_->replaceLocalValues(local_row, &local_row, diagValues.getRawPtr(), 1);
414 L_block_->replaceLocalValues(local_row, &LI[0], &LV[0], NumL);
418 U_block_->replaceLocalValues(local_row, &UI[0], &UV[0], NumU);
433 this->isInitialized_ =
true;
442 template<
class LittleBlockType>
443 struct GetManagedView {
444 static_assert (Kokkos::Impl::is_view<LittleBlockType>::value,
445 "LittleBlockType must be a Kokkos::View.");
446 typedef Kokkos::View<
typename LittleBlockType::non_const_data_type,
447 typename LittleBlockType::array_layout,
448 typename LittleBlockType::device_type> managed_non_const_type;
449 static_assert (static_cast<int> (managed_non_const_type::rank) ==
450 static_cast<int> (LittleBlockType::rank),
451 "managed_non_const_type::rank != LittleBlock::rank. " 452 "Please report this bug to the Ifpack2 developers.");
457 template<
class MatrixType>
460 typedef impl_scalar_type IST;
461 const char prefix[] =
"Ifpack2::Experimental::RBILUK::compute: ";
466 TEUCHOS_TEST_FOR_EXCEPTION
467 (A_block_.is_null (), std::runtime_error, prefix <<
"The matrix (A_block_, " 468 "the BlockCrsMatrix) is null. Please call setMatrix() with a nonnull " 469 "input before calling this method.");
470 TEUCHOS_TEST_FOR_EXCEPTION
471 (! A_block_->isFillComplete (), std::runtime_error, prefix <<
"The matrix " 472 "(A_block_, the BlockCrsMatrix) is not fill complete. You may not invoke " 473 "initialize() or compute() with this matrix until the matrix is fill " 474 "complete. Note: BlockCrsMatrix is fill complete if and only if its " 475 "underlying graph is fill complete.");
477 if (! this->isInitialized ()) {
484 if (! A_block_.is_null ()) {
485 Teuchos::RCP<block_crs_matrix_type> A_nc =
486 Teuchos::rcp_const_cast<block_crs_matrix_type> (A_block_);
499 Teuchos::Time timer (
"RBILUK::compute");
500 double startTime = timer.wallTime();
502 Teuchos::TimeMonitor timeMon (timer);
503 this->isComputed_ =
false;
510 initAllValues (*A_block_);
513 LO NumL, NumU, NumURead;
516 const size_t MaxNumEntries =
517 L_block_->getNodeMaxNumRowEntries () + U_block_->getNodeMaxNumRowEntries () + 1;
519 const LO blockMatSize = blockSize_*blockSize_;
524 const LO rowStride = blockSize_;
526 Teuchos::Array<int> ipiv_teuchos(blockSize_);
527 Kokkos::View<
int*, Kokkos::HostSpace,
528 Kokkos::MemoryUnmanaged> ipiv (ipiv_teuchos.getRawPtr (), blockSize_);
529 Teuchos::Array<IST> work_teuchos(blockSize_);
530 Kokkos::View<IST*, Kokkos::HostSpace,
531 Kokkos::MemoryUnmanaged> work (work_teuchos.getRawPtr (), blockSize_);
533 size_t num_cols = U_block_->getColMap()->getNodeNumElements();
534 Teuchos::Array<int> colflag(num_cols);
536 typename GetManagedView<little_block_host_type>::managed_non_const_type diagModBlock (
"diagModBlock", blockSize_, blockSize_);
537 typename GetManagedView<little_block_host_type>::managed_non_const_type matTmp (
"matTmp", blockSize_, blockSize_);
538 typename GetManagedView<little_block_host_type>::managed_non_const_type multiplier (
"multiplier", blockSize_, blockSize_);
546 for (
size_t j = 0; j < num_cols; ++j) {
549 Teuchos::Array<LO> InI(MaxNumEntries, 0);
550 Teuchos::Array<scalar_type> InV(MaxNumEntries*blockMatSize,STM::zero());
552 const LO numLocalRows = L_block_->getNodeNumRows ();
553 for (LO local_row = 0; local_row < numLocalRows; ++local_row) {
557 NumIn = MaxNumEntries;
558 local_inds_host_view_type colValsL;
559 values_host_view_type valsL;
561 L_block_->getLocalRowView(local_row, colValsL, valsL);
562 NumL = (LO) colValsL.size();
563 for (LO j = 0; j < NumL; ++j)
565 const LO matOffset = blockMatSize*j;
566 little_block_host_type lmat((
typename little_block_host_type::value_type*) &valsL[matOffset],blockSize_,rowStride);
567 little_block_host_type lmatV((
typename little_block_host_type::value_type*) &InV[matOffset],blockSize_,rowStride);
569 Tpetra::COPY (lmat, lmatV);
570 InI[j] = colValsL[j];
573 little_block_host_type dmat = D_block_->getLocalBlockHostNonConst(local_row, local_row);
574 little_block_host_type dmatV((
typename little_block_host_type::value_type*) &InV[NumL*blockMatSize], blockSize_, rowStride);
576 Tpetra::COPY (dmat, dmatV);
577 InI[NumL] = local_row;
579 local_inds_host_view_type colValsU;
580 values_host_view_type valsU;
581 U_block_->getLocalRowView(local_row, colValsU, valsU);
582 NumURead = (LO) colValsU.size();
584 for (LO j = 0; j < NumURead; ++j)
586 if (!(colValsU[j] < numLocalRows))
continue;
587 InI[NumL+1+j] = colValsU[j];
588 const LO matOffset = blockMatSize*(NumL+1+j);
589 little_block_host_type umat((
typename little_block_host_type::value_type*) &valsU[blockMatSize*j], blockSize_, rowStride);
590 little_block_host_type umatV((
typename little_block_host_type::value_type*) &InV[matOffset], blockSize_, rowStride);
592 Tpetra::COPY (umat, umatV);
598 for (
size_t j = 0; j < NumIn; ++j) {
602 #ifndef IFPACK2_RBILUK_INITIAL 603 for (LO i = 0; i < blockSize_; ++i)
604 for (LO j = 0; j < blockSize_; ++j){
606 diagModBlock(i,j) = 0;
611 Kokkos::deep_copy (diagModBlock, diagmod);
614 for (LO jj = 0; jj < NumL; ++jj) {
616 little_block_host_type currentVal((
typename little_block_host_type::value_type*) &InV[jj*blockMatSize], blockSize_, rowStride);
618 Tpetra::COPY (currentVal, multiplier);
620 const little_block_host_type dmatInverse = D_block_->getLocalBlockHostNonConst(j,j);
622 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 623 KokkosBatched::Experimental::SerialGemm
624 <KokkosBatched::Experimental::Trans::NoTranspose,
625 KokkosBatched::Experimental::Trans::NoTranspose,
626 KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
627 (STS::one (), currentVal, dmatInverse, STS::zero (), matTmp);
629 Tpetra::GEMM (
"N",
"N", STS::one (), currentVal, dmatInverse,
630 STS::zero (), matTmp);
634 Tpetra::COPY (matTmp, currentVal);
635 local_inds_host_view_type UUI;
636 values_host_view_type UUV;
638 U_block_->getLocalRowView(j, UUI, UUV);
639 NumUU = (LO) UUI.size();
641 if (this->RelaxValue_ == STM::zero ()) {
642 for (LO k = 0; k < NumUU; ++k) {
643 if (!(UUI[k] < numLocalRows))
continue;
644 const int kk = colflag[UUI[k]];
646 little_block_host_type kkval((
typename little_block_host_type::value_type*) &InV[kk*blockMatSize], blockSize_, rowStride);
647 little_block_host_type uumat((
typename little_block_host_type::value_type*) &UUV[k*blockMatSize], blockSize_, rowStride);
648 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 649 KokkosBatched::Experimental::SerialGemm
650 <KokkosBatched::Experimental::Trans::NoTranspose,
651 KokkosBatched::Experimental::Trans::NoTranspose,
652 KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
653 (
magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), kkval);
655 Tpetra::GEMM (
"N",
"N",
magnitude_type(-STM::one ()), multiplier, uumat,
663 for (LO k = 0; k < NumUU; ++k) {
664 if (!(UUI[k] < numLocalRows))
continue;
665 const int kk = colflag[UUI[k]];
666 little_block_host_type uumat((
typename little_block_host_type::value_type*) &UUV[k*blockMatSize], blockSize_, rowStride);
668 little_block_host_type kkval((
typename little_block_host_type::value_type*) &InV[kk*blockMatSize], blockSize_, rowStride);
669 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 670 KokkosBatched::Experimental::SerialGemm
671 <KokkosBatched::Experimental::Trans::NoTranspose,
672 KokkosBatched::Experimental::Trans::NoTranspose,
673 KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
674 (
magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), kkval);
676 Tpetra::GEMM (
"N",
"N",
magnitude_type(-STM::one ()), multiplier, uumat,
682 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 683 KokkosBatched::Experimental::SerialGemm
684 <KokkosBatched::Experimental::Trans::NoTranspose,
685 KokkosBatched::Experimental::Trans::NoTranspose,
686 KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
687 (
magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), diagModBlock);
689 Tpetra::GEMM (
"N",
"N",
magnitude_type(-STM::one ()), multiplier, uumat,
690 STM::one (), diagModBlock);
699 L_block_->replaceLocalValues (local_row, InI.getRawPtr (), InV.getRawPtr (), NumL);
703 Tpetra::COPY (dmatV, dmat);
705 if (this->RelaxValue_ != STM::zero ()) {
707 Tpetra::AXPY (this->RelaxValue_, diagModBlock, dmat);
721 for (
int k = 0; k < blockSize_; ++k) {
725 Tpetra::GETF2 (dmat, ipiv, lapackInfo);
727 TEUCHOS_TEST_FOR_EXCEPTION(
728 lapackInfo != 0, std::runtime_error,
"Ifpack2::Experimental::RBILUK::compute: " 729 "lapackInfo = " << lapackInfo <<
" which indicates an error in the factorization GETRF.");
731 Tpetra::GETRI (dmat, ipiv, work, lapackInfo);
733 TEUCHOS_TEST_FOR_EXCEPTION(
734 lapackInfo != 0, std::runtime_error,
"Ifpack2::Experimental::RBILUK::compute: " 735 "lapackInfo = " << lapackInfo <<
" which indicates an error in the matrix inverse GETRI.");
738 for (LO j = 0; j < NumU; ++j) {
739 little_block_host_type currentVal((
typename little_block_host_type::value_type*) &InV[(NumL+1+j)*blockMatSize], blockSize_, rowStride);
741 #ifndef IFPACK2_RBILUK_INITIAL_NOKK 742 KokkosBatched::Experimental::SerialGemm
743 <KokkosBatched::Experimental::Trans::NoTranspose,
744 KokkosBatched::Experimental::Trans::NoTranspose,
745 KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
746 (STS::one (), dmat, currentVal, STS::zero (), matTmp);
748 Tpetra::GEMM (
"N",
"N", STS::one (), dmat, currentVal,
749 STS::zero (), matTmp);
753 Tpetra::COPY (matTmp, currentVal);
758 U_block_->replaceLocalValues (local_row, &InI[NumL+1], &InV[blockMatSize*(NumL+1)], NumU);
761 #ifndef IFPACK2_RBILUK_INITIAL 763 for (
size_t j = 0; j < NumIn; ++j) {
764 colflag[InI[j]] = -1;
768 for (
size_t j = 0; j < num_cols; ++j) {
791 this->isComputed_ =
true;
792 this->numCompute_ += 1;
793 this->computeTime_ += (timer.wallTime() - startTime);
797 template<
class MatrixType>
800 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
801 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
802 Teuchos::ETransp mode,
812 TEUCHOS_TEST_FOR_EXCEPTION(
813 A_block_.is_null (), std::runtime_error,
"Ifpack2::Experimental::RBILUK::apply: The matrix is " 814 "null. Please call setMatrix() with a nonnull input, then initialize() " 815 "and compute(), before calling this method.");
816 TEUCHOS_TEST_FOR_EXCEPTION(
817 ! this->isComputed (), std::runtime_error,
818 "Ifpack2::Experimental::RBILUK::apply: If you have not yet called compute(), " 819 "you must call compute() before calling this method.");
820 TEUCHOS_TEST_FOR_EXCEPTION(
821 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
822 "Ifpack2::Experimental::RBILUK::apply: X and Y do not have the same number of columns. " 823 "X.getNumVectors() = " << X.getNumVectors ()
824 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
825 TEUCHOS_TEST_FOR_EXCEPTION(
826 STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
827 "Ifpack2::Experimental::RBILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for " 828 "complex Scalar type. Please talk to the Ifpack2 developers to get this " 829 "fixed. There is a FIXME in this file about this very issue.");
831 const LO blockMatSize = blockSize_*blockSize_;
833 const LO rowStride = blockSize_;
835 BMV yBlock (Y, * (A_block_->getGraph ()->getDomainMap ()), blockSize_);
836 const BMV xBlock (X, * (A_block_->getColMap ()), blockSize_);
838 Teuchos::Array<scalar_type> lclarray(blockSize_);
839 little_host_vec_type lclvec((
typename little_host_vec_type::value_type*)&lclarray[0], blockSize_);
843 Teuchos::Time timer (
"RBILUK::apply");
844 double startTime = timer.wallTime();
846 Teuchos::TimeMonitor timeMon (timer);
847 if (alpha == one && beta == zero) {
848 if (mode == Teuchos::NO_TRANS) {
855 const LO numVectors = xBlock.getNumVectors();
856 BMV cBlock (* (A_block_->getGraph ()->getDomainMap ()), blockSize_, numVectors);
857 BMV rBlock (* (A_block_->getGraph ()->getDomainMap ()), blockSize_, numVectors);
858 for (LO imv = 0; imv < numVectors; ++imv)
860 for (
size_t i = 0; i < D_block_->getNodeNumRows(); ++i)
863 const_host_little_vec_type xval =
864 xBlock.getLocalBlockHost(local_row, imv, Tpetra::Access::ReadOnly);
865 little_host_vec_type cval =
866 cBlock.getLocalBlockHost(local_row, imv, Tpetra::Access::OverwriteAll);
868 Tpetra::COPY (xval, cval);
870 local_inds_host_view_type colValsL;
871 values_host_view_type valsL;
872 L_block_->getLocalRowView(local_row, colValsL, valsL);
873 LO NumL = (LO) colValsL.size();
875 for (LO j = 0; j < NumL; ++j)
877 LO col = colValsL[j];
878 const_host_little_vec_type prevVal =
879 cBlock.getLocalBlockHost(col, imv, Tpetra::Access::ReadOnly);
881 const LO matOffset = blockMatSize*j;
882 little_block_host_type lij((
typename little_block_host_type::value_type*) &valsL[matOffset],blockSize_,rowStride);
885 Tpetra::GEMV (-one, lij, prevVal, cval);
891 D_block_->applyBlock(cBlock, rBlock);
894 for (LO imv = 0; imv < numVectors; ++imv)
896 const LO numRows = D_block_->getNodeNumRows();
897 for (LO i = 0; i < numRows; ++i)
899 LO local_row = (numRows-1)-i;
900 const_host_little_vec_type rval =
901 rBlock.getLocalBlockHost(local_row, imv, Tpetra::Access::ReadOnly);
902 little_host_vec_type yval =
903 yBlock.getLocalBlockHost(local_row, imv, Tpetra::Access::OverwriteAll);
905 Tpetra::COPY (rval, yval);
907 local_inds_host_view_type colValsU;
908 values_host_view_type valsU;
909 U_block_->getLocalRowView(local_row, colValsU, valsU);
910 LO NumU = (LO) colValsU.size();
912 for (LO j = 0; j < NumU; ++j)
914 LO col = colValsU[NumU-1-j];
915 const_host_little_vec_type prevVal =
916 yBlock.getLocalBlockHost(col, imv, Tpetra::Access::ReadOnly);
918 const LO matOffset = blockMatSize*(NumU-1-j);
919 little_block_host_type uij((
typename little_block_host_type::value_type*) &valsU[matOffset], blockSize_, rowStride);
922 Tpetra::GEMV (-one, uij, prevVal, yval);
928 TEUCHOS_TEST_FOR_EXCEPTION(
929 true, std::runtime_error,
930 "Ifpack2::Experimental::RBILUK::apply: transpose apply is not implemented for the block algorithm. ");
941 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
942 apply (X, Y_tmp, mode);
943 Y.update (alpha, Y_tmp, beta);
948 this->numApply_ += 1;
949 this->applyTime_ += (timer.wallTime() - startTime);
953 template<
class MatrixType>
956 std::ostringstream os;
961 os <<
"\"Ifpack2::Experimental::RBILUK\": {";
962 os <<
"Initialized: " << (this->isInitialized () ?
"true" :
"false") <<
", " 963 <<
"Computed: " << (this->isComputed () ?
"true" :
"false") <<
", ";
965 os <<
"Level-of-fill: " << this->getLevelOfFill() <<
", ";
967 if (A_block_.is_null ()) {
968 os <<
"Matrix: null";
971 os <<
"Global matrix dimensions: [" 972 << A_block_->getGlobalNumRows () <<
", " << A_block_->getGlobalNumCols () <<
"]" 973 <<
", Global nnz: " << A_block_->getGlobalNumEntries();
988 #define IFPACK2_EXPERIMENTAL_RBILUK_INSTANT(S,LO,GO,N) \ 989 template class Ifpack2::Experimental::RBILUK< Tpetra::BlockCrsMatrix<S, LO, GO, N> >; \ 990 template class Ifpack2::Experimental::RBILUK< Tpetra::RowMatrix<S, LO, GO, N> >; const block_crs_matrix_type & getDBlock() const
Return the diagonal entries of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:125
Ifpack2 features that are experimental. Use at your own risk.
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:187
Teuchos::RCP< const block_crs_matrix_type > getBlockMatrix() const
Get the input matrix.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:182
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:146
virtual ~RBILUK()
Destructor (declared virtual for memory safety).
Definition: Ifpack2_Experimental_RBILUK_def.hpp:81
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:136
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:153
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:245
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:150
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 (inverse of the) incomplete factorization to X, resulting in Y.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:800
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:142
void compute()
Compute the (numeric) incomplete factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:458
File for utility functions.
const block_crs_matrix_type & getUBlock() const
Return the U factor of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:139
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:100
void setMatrix(const Teuchos::RCP< const block_crs_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:86
Definition: Ifpack2_Container_decl.hpp:576
std::string description() const
A one-line description of this object.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:954
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:58
const block_crs_matrix_type & getLBlock() const
Return the L factor of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:111
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:159
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
ILU(k) factorization of a given Tpetra::BlockCrsMatrix.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:128