43 #ifndef IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP 44 #define IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP 46 #include "Ifpack2_LocalFilter.hpp" 47 #include "Teuchos_LAPACK.hpp" 51 # include "Teuchos_DefaultMpiComm.hpp" 53 # include "Teuchos_DefaultSerialComm.hpp" 64 template<
class MatrixType>
68 initializeTime_ (0.0),
74 isInitialized_ (false),
79 template<
class MatrixType>
80 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
82 TEUCHOS_TEST_FOR_EXCEPTION(
83 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 84 "getDomainMap: The input matrix A is null. Please call setMatrix() with a " 85 "nonnull input matrix before calling this method.");
88 return A_->getRangeMap ();
92 template<
class MatrixType>
93 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
95 TEUCHOS_TEST_FOR_EXCEPTION(
96 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 97 "getRangeMap: The input matrix A is null. Please call setMatrix() with a " 98 "nonnull input matrix before calling this method.");
101 return A_->getDomainMap ();
105 template<
class MatrixType>
113 template<
class MatrixType>
116 return isInitialized_;
120 template<
class MatrixType>
127 template<
class MatrixType>
130 return numInitialize_;
134 template<
class MatrixType>
141 template<
class MatrixType>
148 template<
class MatrixType>
151 return initializeTime_;
155 template<
class MatrixType>
162 template<
class MatrixType>
169 template<
class MatrixType>
170 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::row_matrix_type>
176 template<
class MatrixType>
180 isInitialized_ =
false;
182 A_local_ = Teuchos::null;
183 A_local_tridi_.reshape (0);
188 template<
class MatrixType>
190 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
194 if (! A_.is_null ()) {
195 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
196 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
197 TEUCHOS_TEST_FOR_EXCEPTION(
198 numRows != numCols, std::invalid_argument,
"Ifpack2::Details::TriDiSolver::" 199 "setMatrix: Input matrix must be (globally) square. " 200 "The matrix you provided is " << numRows <<
" by " << numCols <<
".");
210 template<
class MatrixType>
218 using Teuchos::TimeMonitor;
219 const std::string timerName (
"Ifpack2::Details::TriDiSolver::initialize");
221 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
222 if (timer.is_null ()) {
223 timer = TimeMonitor::getNewCounter (timerName);
227 Teuchos::TimeMonitor timeMon (*timer);
229 TEUCHOS_TEST_FOR_EXCEPTION(
230 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 231 "initialize: The input matrix A is null. Please call setMatrix() " 232 "with a nonnull input before calling this method.");
234 TEUCHOS_TEST_FOR_EXCEPTION(
235 ! A_->hasColMap (), std::invalid_argument,
"Ifpack2::Details::TriDiSolver: " 236 "The constructor's input matrix must have a column Map, " 237 "so that it has local indices.");
243 if (A_->getComm ()->getSize () > 1) {
249 TEUCHOS_TEST_FOR_EXCEPTION(
250 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::" 251 "initialize: A_local_ is null after it was supposed to have been " 252 "initialized. Please report this bug to the Ifpack2 developers.");
255 const size_t numRows = A_local_->getNodeNumRows ();
256 const size_t numCols = A_local_->getNodeNumCols ();
257 TEUCHOS_TEST_FOR_EXCEPTION(
258 numRows != numCols, std::logic_error,
"Ifpack2::Details::TriDiSolver::" 259 "initialize: Local filter matrix is not square. This should never happen. " 260 "Please report this bug to the Ifpack2 developers.");
261 A_local_tridi_.reshape (numRows);
262 ipiv_.resize (numRows);
263 std::fill (ipiv_.begin (), ipiv_.end (), 0);
265 isInitialized_ =
true;
271 initializeTime_ = timer->totalElapsedTime ();
275 template<
class MatrixType>
280 template<
class MatrixType>
284 const std::string timerName (
"Ifpack2::Details::TriDiSolver::compute");
286 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
287 if (timer.is_null ()) {
288 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
293 Teuchos::TimeMonitor timeMon (*timer);
294 TEUCHOS_TEST_FOR_EXCEPTION(
295 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 296 "compute: The input matrix A is null. Please call setMatrix() with a " 297 "nonnull input, then call initialize(), before calling this method.");
299 TEUCHOS_TEST_FOR_EXCEPTION(
300 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::" 301 "compute: A_local_ is null. Please report this bug to the Ifpack2 " 308 extract (A_local_tridi_, *A_local_);
310 factor (A_local_tridi_, ipiv_ ());
317 computeTime_ = timer->totalElapsedTime ();
320 template<
class MatrixType>
322 const Teuchos::ArrayView<int>& ipiv)
325 std::fill (ipiv.begin (), ipiv.end (), 0);
327 Teuchos::LAPACK<int, scalar_type> lapack;
329 lapack.GTTRF (A.numRowsCols (),
334 ipiv.getRawPtr (), &INFO);
336 TEUCHOS_TEST_FOR_EXCEPTION(
337 INFO < 0, std::logic_error,
"Ifpack2::Details::TriDiSolver::factor: " 338 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) " 339 "was called incorrectly. INFO = " << INFO <<
" < 0. " 340 "Please report this bug to the Ifpack2 developers.");
344 TEUCHOS_TEST_FOR_EXCEPTION(
345 INFO > 0, std::runtime_error,
"Ifpack2::Details::TriDiSolver::factor: " 346 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) " 347 "reports that the computed U factor is exactly singular. U(" << INFO <<
348 "," << INFO <<
") (one-based index i) is exactly zero. This probably " 349 "means that the input matrix has a singular diagonal block.");
353 template<
class MatrixType>
354 void TriDiSolver<MatrixType, false>::
355 applyImpl (
const MV& X,
357 const Teuchos::ETransp mode,
358 const scalar_type alpha,
359 const scalar_type beta)
const 361 using Teuchos::ArrayRCP;
364 using Teuchos::rcpFromRef;
365 using Teuchos::CONJ_TRANS;
366 using Teuchos::TRANS;
368 const int numVecs =
static_cast<int> (X.getNumVectors ());
369 if (alpha == STS::zero ()) {
370 if (beta == STS::zero ()) {
374 Y.putScalar (STS::zero ());
377 Y.scale (STS::zero ());
381 Teuchos::LAPACK<int, scalar_type> lapack;
387 if (beta == STS::zero () && Y.isConstantStride () && alpha == STS::one ()) {
389 Y_tmp = rcpFromRef (Y);
392 Y_tmp = rcp (
new MV (createCopy(X)));
393 if (alpha != STS::one ()) {
394 Y_tmp->scale (alpha);
397 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
398 ArrayRCP<scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
399 scalar_type*
const Y_ptr = Y_view.getRawPtr ();
402 (mode == CONJ_TRANS ?
'C' : (mode ==
TRANS ?
'T' :
'N'));
403 lapack.GTTRS (trans, A_local_tridi_.numRowsCols(), numVecs,
407 A_local_tridi_.DU2(),
408 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
409 TEUCHOS_TEST_FOR_EXCEPTION(
410 INFO != 0, std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 411 "applyImpl: LAPACK's _GTTRS (tridiagonal solve using LU factorization " 412 "with partial pivoting) failed with INFO = " << INFO <<
" != 0.");
414 if (beta != STS::zero ()) {
415 Y.update (alpha, *Y_tmp, beta);
417 else if (! Y.isConstantStride ()) {
418 deep_copy(Y, *Y_tmp);
424 template<
class MatrixType>
426 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
427 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
428 Teuchos::ETransp mode,
432 using Teuchos::ArrayView;
436 using Teuchos::rcpFromRef;
438 const std::string timerName (
"Ifpack2::Details::TriDiSolver::apply");
439 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
440 if (timer.is_null ()) {
441 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
446 Teuchos::TimeMonitor timeMon (*timer);
448 TEUCHOS_TEST_FOR_EXCEPTION(
449 ! isComputed_, std::runtime_error,
"Ifpack2::Details::TriDiSolver::apply: " 450 "You must have called the compute() method before you may call apply(). " 451 "You may call the apply() method as many times as you want after calling " 452 "compute() once, but you must have called compute() at least once.");
454 const size_t numVecs = X.getNumVectors ();
456 TEUCHOS_TEST_FOR_EXCEPTION(
457 numVecs != Y.getNumVectors (), std::runtime_error,
458 "Ifpack2::Details::TriDiSolver::apply: X and Y have different numbers " 459 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 460 << X.getNumVectors () <<
".");
467 RCP<const MV> X_local;
469 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
474 X_local = X.offsetView (A_local_->getDomainMap (), 0);
475 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
479 X_local = rcpFromRef (X);
480 Y_local = rcpFromRef (Y);
485 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
492 applyTime_ = timer->totalElapsedTime ();
496 template<
class MatrixType>
499 std::ostringstream out;
504 out <<
"\"Ifpack2::Details::TriDiSolver\": ";
506 if (this->getObjectLabel () !=
"") {
507 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
509 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", " 510 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
513 out <<
"Matrix: null";
516 out <<
"Matrix: not null" 517 <<
", Global matrix dimensions: [" 518 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
526 template<
class MatrixType>
528 const Teuchos::EVerbosityLevel verbLevel)
const {
529 using Teuchos::FancyOStream;
530 using Teuchos::OSTab;
532 using Teuchos::rcpFromRef;
535 if (verbLevel == Teuchos::VERB_NONE) {
539 RCP<FancyOStream> ptrOut = rcpFromRef (out);
541 if (this->getObjectLabel () !=
"") {
542 out <<
"label: " << this->getObjectLabel () << endl;
544 out <<
"initialized: " << (isInitialized_ ?
"true" :
"false") << endl
545 <<
"computed: " << (isComputed_ ?
"true" :
"false") << endl
546 <<
"number of initialize calls: " << numInitialize_ << endl
547 <<
"number of compute calls: " << numCompute_ << endl
548 <<
"number of apply calls: " << numApply_ << endl
549 <<
"total time in seconds in initialize: " << initializeTime_ << endl
550 <<
"total time in seconds in compute: " << computeTime_ << endl
551 <<
"total time in seconds in apply: " << applyTime_ << endl;
552 if (verbLevel >= Teuchos::VERB_EXTREME) {
553 out <<
"A_local_tridi_:" << endl;
554 A_local_tridi_.print(out);
556 out <<
"ipiv_: " << Teuchos::toString (ipiv_) << endl;
560 template<
class MatrixType>
562 const Teuchos::EVerbosityLevel verbLevel)
const 564 using Teuchos::FancyOStream;
565 using Teuchos::OSTab;
567 using Teuchos::rcpFromRef;
570 RCP<FancyOStream> ptrOut = rcpFromRef (out);
577 if (verbLevel > Teuchos::VERB_NONE) {
578 out <<
"Ifpack2::Details::TriDiSolver:" << endl;
580 describeLocal (out, verbLevel);
585 const Teuchos::Comm<int>& comm = * (A_->getRowMap ()->getComm ());
586 const int myRank = comm.getRank ();
587 const int numProcs = comm.getSize ();
588 if (verbLevel > Teuchos::VERB_NONE && myRank == 0) {
589 out <<
"Ifpack2::Details::TriDiSolver:" << endl;
592 for (
int p = 0; p < numProcs; ++p) {
594 out <<
"Process " << myRank <<
":" << endl;
595 describeLocal (out, verbLevel);
604 template<
class MatrixType>
606 const row_matrix_type& A_local)
608 using Teuchos::Array;
609 using Teuchos::ArrayView;
610 typedef local_ordinal_type LO;
611 typedef typename Teuchos::ArrayView<LO>::size_type size_type;
614 A_local_tridi.putScalar (STS::zero ());
622 const map_type& rowMap = * (A_local.getRowMap ());
626 const size_type maxNumRowEntries =
627 static_cast<size_type
> (A_local.getNodeMaxNumRowEntries ());
628 Array<LO> localIndices (maxNumRowEntries);
629 Array<scalar_type> values (maxNumRowEntries);
631 const LO numLocalRows =
static_cast<LO
> (rowMap.getNodeNumElements ());
632 const LO minLocalRow = rowMap.getMinLocalIndex ();
635 const LO maxLocalRow = minLocalRow + numLocalRows;
636 for (LO localRow = minLocalRow; localRow < maxLocalRow; ++localRow) {
643 const size_type numEntriesInRow =
644 static_cast<size_type
> (A_local.getNumEntriesInLocalRow (localRow));
645 size_t numEntriesOut = 0;
646 A_local.getLocalRowCopy (localRow,
647 localIndices (0, numEntriesInRow),
648 values (0, numEntriesInRow),
650 for (LO k = 0; k < numEntriesInRow; ++k) {
651 const LO localCol = localIndices[k];
652 const scalar_type val = values[k];
656 if( localCol >= localRow-1 && localCol <= localRow+1 )
657 A_local_tridi(localRow, localCol) += val;
666 template<
class MatrixType>
667 TriDiSolver<MatrixType, true>::TriDiSolver (
const Teuchos::RCP<const row_matrix_type>& A) {
668 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
672 template<
class MatrixType>
673 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
675 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
679 template<
class MatrixType>
680 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
682 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
686 template<
class MatrixType>
689 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
693 template<
class MatrixType>
696 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
700 template<
class MatrixType>
703 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
707 template<
class MatrixType>
710 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
714 template<
class MatrixType>
717 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
721 template<
class MatrixType>
724 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
728 template<
class MatrixType>
731 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
735 template<
class MatrixType>
738 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
742 template<
class MatrixType>
745 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
749 template<
class MatrixType>
750 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::row_matrix_type>
752 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
756 template<
class MatrixType>
759 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
763 template<
class MatrixType>
766 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
770 template<
class MatrixType>
771 TriDiSolver<MatrixType, true>::~TriDiSolver ()
777 template<
class MatrixType>
780 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
784 template<
class MatrixType>
786 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
787 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
788 Teuchos::ETransp mode,
790 scalar_type beta)
const 792 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
796 template<
class MatrixType>
798 TriDiSolver<MatrixType, true>::description ()
const 800 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
804 template<
class MatrixType>
805 void TriDiSolver<MatrixType, true>::describe(Teuchos::FancyOStream& out,
806 const Teuchos::EVerbosityLevel verbLevel)
const 808 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
814 #define IFPACK2_DETAILS_TRIDISOLVER_INSTANT(S,LO,GO,N) \ 815 template class Ifpack2::Details::TriDiSolver< Tpetra::RowMatrix<S, LO, GO, N> >; 817 #endif // IFPACK2_DETAILS_TRIDISOLVER_HPP virtual Teuchos::RCP< const Tpetra::Map< MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getDomainMap() const=0
The domain Map of this operator.
virtual int getNumApply() const=0
The number of calls to apply().
virtual Teuchos::RCP< const Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getMatrix() const=0
The input matrix given to the constructor.
"Preconditioner" that uses LAPACK's tridi LU.
Definition: Ifpack2_Details_TriDiSolver_decl.hpp:74
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_Details_TriDiSolver_decl.hpp:107
virtual void compute()=0
Set up the numerical values in this preconditioner.
virtual Teuchos::RCP< const Tpetra::Map< MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > getRangeMap() const=0
The range Map of this operator.
Ifpack2 implementation details.
virtual int getNumCompute() const=0
The number of calls to compute().
TriDiSolver(const Teuchos::RCP< const row_matrix_type > &matrix)
Constructor.
Definition: Ifpack2_Details_TriDiSolver_def.hpp:66
virtual double getInitializeTime() const=0
The time (in seconds) spent in initialize().
virtual void initialize()=0
Set up the graph structure of this preconditioner.
virtual bool isComputed() const=0
True if the preconditioner has been successfully computed, else false.
virtual int getNumInitialize() const=0
The number of calls to initialize().
virtual void setParameters(const Teuchos::ParameterList &List)=0
Set this preconditioner's parameters.
virtual double getComputeTime() const=0
The time (in seconds) spent in compute().
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual void apply(const Tpetra::MultiVector< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > &X, Tpetra::MultiVector< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, MatrixType::scalar_type alpha=Teuchos::ScalarTraits< MatrixType::scalar_type >::one(), MatrixType::scalar_type beta=Teuchos::ScalarTraits< MatrixType::scalar_type >::zero()) const=0
Apply the preconditioner to X, putting the result in Y.
virtual bool isInitialized() const=0
True if the preconditioner has been successfully initialized, else false.
virtual double getApplyTime() const=0
The time (in seconds) spent in apply().
virtual void setMatrix(const Teuchos::RCP< const Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > > &A)=0
Set the new matrix.