43 #ifndef IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP
44 #define IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP
46 #include "Ifpack2_LocalFilter.hpp"
48 #include "Tpetra_MultiVector.hpp"
49 #include "Tpetra_Map.hpp"
50 #include "Tpetra_Import.hpp"
51 #include "Tpetra_Export.hpp"
57 # include "Teuchos_DefaultSerialComm.hpp"
68 template<
class MatrixType>
72 initializeTime_ (0.0),
78 isInitialized_ (false),
83 template<
class MatrixType>
87 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::"
88 "getDomainMap: The input matrix A is null. Please call setMatrix() with a "
89 "nonnull input matrix before calling this method.");
92 return A_->getRangeMap ();
96 template<
class MatrixType>
100 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::"
101 "getRangeMap: The input matrix A is null. Please call setMatrix() with a "
102 "nonnull input matrix before calling this method.");
105 return A_->getDomainMap ();
109 template<
class MatrixType>
117 template<
class MatrixType>
120 return isInitialized_;
124 template<
class MatrixType>
131 template<
class MatrixType>
134 return numInitialize_;
138 template<
class MatrixType>
145 template<
class MatrixType>
152 template<
class MatrixType>
155 return initializeTime_;
159 template<
class MatrixType>
166 template<
class MatrixType>
173 template<
class MatrixType>
180 template<
class MatrixType>
184 isInitialized_ =
false;
186 A_local_ = Teuchos::null;
187 A_local_tridi_.reshape (0);
192 template<
class MatrixType>
198 if (! A_.is_null ()) {
199 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
200 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
202 numRows != numCols, std::invalid_argument,
"Ifpack2::Details::TriDiSolver::"
203 "setMatrix: Input matrix must be (globally) square. "
204 "The matrix you provided is " << numRows <<
" by " << numCols <<
".");
214 template<
class MatrixType>
223 const std::string timerName (
"Ifpack2::Details::TriDiSolver::initialize");
225 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
226 if (timer.is_null ()) {
227 timer = TimeMonitor::getNewCounter (timerName);
234 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::"
235 "initialize: The input matrix A is null. Please call setMatrix() "
236 "with a nonnull input before calling this method.");
239 ! A_->hasColMap (), std::invalid_argument,
"Ifpack2::Details::TriDiSolver: "
240 "The constructor's input matrix must have a column Map, "
241 "so that it has local indices.");
247 if (A_->getComm ()->getSize () > 1) {
254 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::"
255 "initialize: A_local_ is null after it was supposed to have been "
256 "initialized. Please report this bug to the Ifpack2 developers.");
259 const size_t numRows = A_local_->getNodeNumRows ();
260 const size_t numCols = A_local_->getNodeNumCols ();
262 numRows != numCols, std::logic_error,
"Ifpack2::Details::TriDiSolver::"
263 "initialize: Local filter matrix is not square. This should never happen. "
264 "Please report this bug to the Ifpack2 developers.");
265 A_local_tridi_.reshape (numRows);
266 ipiv_.resize (numRows);
267 std::fill (ipiv_.begin (), ipiv_.end (), 0);
269 isInitialized_ =
true;
275 initializeTime_ = timer->totalElapsedTime ();
279 template<
class MatrixType>
284 template<
class MatrixType>
288 const std::string timerName (
"Ifpack2::Details::TriDiSolver::compute");
291 if (timer.is_null ()) {
299 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::"
300 "compute: The input matrix A is null. Please call setMatrix() with a "
301 "nonnull input, then call initialize(), before calling this method.");
304 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::"
305 "compute: A_local_ is null. Please report this bug to the Ifpack2 "
312 extract (A_local_tridi_, *A_local_);
314 factor (A_local_tridi_, ipiv_ ());
321 computeTime_ = timer->totalElapsedTime ();
324 template<
class MatrixType>
329 std::fill (ipiv.
begin (), ipiv.
end (), 0);
341 INFO < 0, std::logic_error,
"Ifpack2::Details::TriDiSolver::factor: "
342 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) "
343 "was called incorrectly. INFO = " << INFO <<
" < 0. "
344 "Please report this bug to the Ifpack2 developers.");
349 INFO > 0, std::runtime_error,
"Ifpack2::Details::TriDiSolver::factor: "
350 "LAPACK's _GTTRF (tridiagonal LU factorization with partial pivoting) "
351 "reports that the computed U factor is exactly singular. U(" << INFO <<
352 "," << INFO <<
") (one-based index i) is exactly zero. This probably "
353 "means that the input matrix has a singular diagonal block.");
357 template<
class MatrixType>
358 void TriDiSolver<MatrixType, false>::
359 applyImpl (
const MV& X,
362 const scalar_type alpha,
363 const scalar_type beta)
const
368 using Teuchos::rcpFromRef;
372 const int numVecs =
static_cast<int> (X.getNumVectors ());
373 if (alpha == STS::zero ()) {
374 if (beta == STS::zero ()) {
378 Y.putScalar (STS::zero ());
381 Y.scale (STS::zero ());
391 if (beta == STS::zero () && Y.isConstantStride () && alpha == STS::one ()) {
393 Y_tmp = rcpFromRef (Y);
396 Y_tmp =
rcp (
new MV (createCopy(X)));
397 if (alpha != STS::one ()) {
398 Y_tmp->scale (alpha);
401 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
402 ArrayRCP<scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
403 scalar_type*
const Y_ptr = Y_view.getRawPtr ();
407 lapack.
GTTRS (trans, A_local_tridi_.numRowsCols(), numVecs,
411 A_local_tridi_.DU2(),
412 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
414 INFO != 0, std::runtime_error,
"Ifpack2::Details::TriDiSolver::"
415 "applyImpl: LAPACK's _GTTRS (tridiagonal solve using LU factorization "
416 "with partial pivoting) failed with INFO = " << INFO <<
" != 0.");
418 if (beta != STS::zero ()) {
419 Y.update (alpha, *Y_tmp, beta);
421 else if (! Y.isConstantStride ()) {
422 deep_copy(Y, *Y_tmp);
428 template<
class MatrixType>
430 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
431 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
440 using Teuchos::rcpFromRef;
442 const std::string timerName (
"Ifpack2::Details::TriDiSolver::apply");
444 if (timer.is_null ()) {
453 ! isComputed_, std::runtime_error,
"Ifpack2::Details::TriDiSolver::apply: "
454 "You must have called the compute() method before you may call apply(). "
455 "You may call the apply() method as many times as you want after calling "
456 "compute() once, but you must have called compute() at least once.");
458 const size_t numVecs = X.getNumVectors ();
461 numVecs != Y.getNumVectors (), std::runtime_error,
462 "Ifpack2::Details::TriDiSolver::apply: X and Y have different numbers "
463 "of vectors. X has " << X.getNumVectors () <<
", but Y has "
464 << X.getNumVectors () <<
".");
471 RCP<const MV> X_local;
473 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
478 X_local = X.offsetView (A_local_->getDomainMap (), 0);
479 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
483 X_local = rcpFromRef (X);
484 Y_local = rcpFromRef (Y);
489 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
496 applyTime_ = timer->totalElapsedTime ();
500 template<
class MatrixType>
503 std::ostringstream out;
508 out <<
"\"Ifpack2::Details::TriDiSolver\": ";
510 if (this->getObjectLabel () !=
"") {
511 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
513 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", "
514 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
517 out <<
"Matrix: null";
520 out <<
"Matrix: not null"
521 <<
", Global matrix dimensions: ["
522 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
530 template<
class MatrixType>
536 using Teuchos::rcpFromRef;
543 RCP<FancyOStream> ptrOut = rcpFromRef (out);
545 if (this->getObjectLabel () !=
"") {
546 out <<
"label: " << this->getObjectLabel () << endl;
548 out <<
"initialized: " << (isInitialized_ ?
"true" :
"false") << endl
549 <<
"computed: " << (isComputed_ ?
"true" :
"false") << endl
550 <<
"number of initialize calls: " << numInitialize_ << endl
551 <<
"number of compute calls: " << numCompute_ << endl
552 <<
"number of apply calls: " << numApply_ << endl
553 <<
"total time in seconds in initialize: " << initializeTime_ << endl
554 <<
"total time in seconds in compute: " << computeTime_ << endl
555 <<
"total time in seconds in apply: " << applyTime_ << endl;
557 out <<
"A_local_tridi_:" << endl;
558 A_local_tridi_.print(out);
564 template<
class MatrixType>
571 using Teuchos::rcpFromRef;
574 RCP<FancyOStream> ptrOut = rcpFromRef (out);
582 out <<
"Ifpack2::Details::TriDiSolver:" << endl;
584 describeLocal (out, verbLevel);
590 const int myRank = comm.
getRank ();
591 const int numProcs = comm.
getSize ();
593 out <<
"Ifpack2::Details::TriDiSolver:" << endl;
596 for (
int p = 0; p < numProcs; ++p) {
598 out <<
"Process " << myRank <<
":" << endl;
599 describeLocal (out, verbLevel);
608 template<
class MatrixType>
610 const row_matrix_type& A_local)
614 typedef local_ordinal_type LO;
626 const map_type& rowMap = * (A_local.getRowMap ());
630 const size_type maxNumRowEntries =
631 static_cast<size_type
> (A_local.getNodeMaxNumRowEntries ());
632 Array<LO> localIndices (maxNumRowEntries);
633 Array<scalar_type> values (maxNumRowEntries);
635 const LO numLocalRows =
static_cast<LO
> (rowMap.getNodeNumElements ());
636 const LO minLocalRow = rowMap.getMinLocalIndex ();
639 const LO maxLocalRow = minLocalRow + numLocalRows;
640 for (LO localRow = minLocalRow; localRow < maxLocalRow; ++localRow) {
647 const size_type numEntriesInRow =
648 static_cast<size_type
> (A_local.getNumEntriesInLocalRow (localRow));
649 size_t numEntriesOut = 0;
650 A_local.getLocalRowCopy (localRow,
651 localIndices (0, numEntriesInRow),
652 values (0, numEntriesInRow),
654 for (LO k = 0; k < numEntriesInRow; ++k) {
655 const LO localCol = localIndices[k];
656 const scalar_type val = values[k];
660 if( localCol >= localRow-1 && localCol <= localRow+1 )
661 A_local_tridi(localRow, localCol) += val;
670 template<
class MatrixType>
676 template<
class MatrixType>
683 template<
class MatrixType>
690 template<
class MatrixType>
697 template<
class MatrixType>
704 template<
class MatrixType>
711 template<
class MatrixType>
718 template<
class MatrixType>
725 template<
class MatrixType>
732 template<
class MatrixType>
739 template<
class MatrixType>
746 template<
class MatrixType>
753 template<
class MatrixType>
760 template<
class MatrixType>
767 template<
class MatrixType>
774 template<
class MatrixType>
775 TriDiSolver<MatrixType, true>::~TriDiSolver ()
781 template<
class MatrixType>
788 template<
class MatrixType>
790 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
791 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
794 scalar_type beta)
const
800 template<
class MatrixType>
802 TriDiSolver<MatrixType, true>::description ()
const
808 template<
class MatrixType>
818 #define IFPACK2_DETAILS_TRIDISOLVER_INSTANT(S,LO,GO,N) \
819 template class Ifpack2::Details::TriDiSolver< Tpetra::RowMatrix<S, LO, GO, N> >;
821 #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 int getSize() const =0
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.
OrdinalType numRowsCols() const
virtual int getRank() const =0
"Preconditioner" that uses LAPACK's tridi LU.
Definition: Ifpack2_Details_TriDiSolver_decl.hpp:76
basic_OSTab< char > OSTab
void GTTRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *dl, const ScalarType *d, const ScalarType *du, const ScalarType *du2, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
void GTTRF(const OrdinalType &n, ScalarType *dl, ScalarType *d, ScalarType *du, ScalarType *du2, OrdinalType *IPIV, OrdinalType *info) const
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_Details_TriDiSolver_decl.hpp:109
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.
virtual void barrier() const =0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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:70
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.
TypeTo as(const TypeFrom &t)
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
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
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_typealpha=Teuchos::ScalarTraits< MatrixType::scalar_type >::one(), MatrixType::scalar_typebeta=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().
std::string toString(const T &t)
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.