43 #ifndef IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP 44 #define IFPACK2_DETAILS_TRIDISOLVER_DEF_HPP 46 #include "Ifpack2_LocalFilter.hpp" 47 #include "Teuchos_LAPACK.hpp" 48 #include "Tpetra_MultiVector.hpp" 49 #include "Tpetra_Map.hpp" 50 #include "Tpetra_Import.hpp" 51 #include "Tpetra_Export.hpp" 55 # include "Teuchos_DefaultMpiComm.hpp" 57 # include "Teuchos_DefaultSerialComm.hpp" 68 template<
class MatrixType>
72 initializeTime_ (0.0),
78 isInitialized_ (false),
83 template<
class MatrixType>
84 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
86 TEUCHOS_TEST_FOR_EXCEPTION(
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>
97 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::map_type>
99 TEUCHOS_TEST_FOR_EXCEPTION(
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>
174 Teuchos::RCP<const typename TriDiSolver<MatrixType, false>::row_matrix_type>
180 template<
class MatrixType>
184 isInitialized_ =
false;
186 A_local_ = Teuchos::null;
187 A_local_tridi_.reshape (0);
192 template<
class MatrixType>
194 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
198 if (! A_.is_null ()) {
199 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
200 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
201 TEUCHOS_TEST_FOR_EXCEPTION(
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>
222 using Teuchos::TimeMonitor;
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);
230 double startTime = timer->wallTime();
233 Teuchos::TimeMonitor timeMon (*timer);
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 237 "initialize: The input matrix A is null. Please call setMatrix() " 238 "with a nonnull input before calling this method.");
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 ! A_->hasColMap (), std::invalid_argument,
"Ifpack2::Details::TriDiSolver: " 242 "The constructor's input matrix must have a column Map, " 243 "so that it has local indices.");
249 if (A_->getComm ()->getSize () > 1) {
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::" 257 "initialize: A_local_ is null after it was supposed to have been " 258 "initialized. Please report this bug to the Ifpack2 developers.");
261 const size_t numRows = A_local_->getNodeNumRows ();
262 const size_t numCols = A_local_->getNodeNumCols ();
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 numRows != numCols, std::logic_error,
"Ifpack2::Details::TriDiSolver::" 265 "initialize: Local filter matrix is not square. This should never happen. " 266 "Please report this bug to the Ifpack2 developers.");
267 A_local_tridi_.reshape (numRows);
268 ipiv_.resize (numRows);
269 std::fill (ipiv_.begin (), ipiv_.end (), 0);
271 isInitialized_ =
true;
275 initializeTime_ += (timer->wallTime() - startTime);
279 template<
class MatrixType>
284 template<
class MatrixType>
288 const std::string timerName (
"Ifpack2::Details::TriDiSolver::compute");
290 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
291 if (timer.is_null ()) {
292 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
295 double startTime = timer->wallTime();
299 Teuchos::TimeMonitor timeMon (*timer);
300 TEUCHOS_TEST_FOR_EXCEPTION(
301 A_.is_null (), std::runtime_error,
"Ifpack2::Details::TriDiSolver::" 302 "compute: The input matrix A is null. Please call setMatrix() with a " 303 "nonnull input, then call initialize(), before calling this method.");
305 TEUCHOS_TEST_FOR_EXCEPTION(
306 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::TriDiSolver::" 307 "compute: A_local_ is null. Please report this bug to the Ifpack2 " 314 extract (A_local_tridi_, *A_local_);
316 factor (A_local_tridi_, ipiv_ ());
321 computeTime_ += (timer->wallTime() - startTime);
324 template<
class MatrixType>
326 const Teuchos::ArrayView<int>& ipiv)
329 std::fill (ipiv.begin (), ipiv.end (), 0);
331 Teuchos::LAPACK<int, scalar_type> lapack;
333 lapack.GTTRF (A.numRowsCols (),
338 ipiv.getRawPtr (), &INFO);
340 TEUCHOS_TEST_FOR_EXCEPTION(
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.");
348 TEUCHOS_TEST_FOR_EXCEPTION(
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,
361 const Teuchos::ETransp mode,
362 const scalar_type alpha,
363 const scalar_type beta)
const 365 using Teuchos::ArrayRCP;
368 using Teuchos::rcpFromRef;
369 using Teuchos::CONJ_TRANS;
370 using Teuchos::TRANS;
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 ());
385 Teuchos::LAPACK<int, scalar_type> lapack;
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 ();
406 (mode == CONJ_TRANS ?
'C' : (mode ==
TRANS ?
'T' :
'N'));
407 lapack.GTTRS (trans, A_local_tridi_.numRowsCols(), numVecs,
411 A_local_tridi_.DU2(),
412 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
413 TEUCHOS_TEST_FOR_EXCEPTION(
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,
432 Teuchos::ETransp mode,
436 using Teuchos::ArrayView;
440 using Teuchos::rcpFromRef;
442 const std::string timerName (
"Ifpack2::Details::TriDiSolver::apply");
443 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
444 if (timer.is_null ()) {
445 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
448 double startTime = timer->wallTime();
452 Teuchos::TimeMonitor timeMon (*timer);
454 TEUCHOS_TEST_FOR_EXCEPTION(
455 ! isComputed_, std::runtime_error,
"Ifpack2::Details::TriDiSolver::apply: " 456 "You must have called the compute() method before you may call apply(). " 457 "You may call the apply() method as many times as you want after calling " 458 "compute() once, but you must have called compute() at least once.");
460 const size_t numVecs = X.getNumVectors ();
462 TEUCHOS_TEST_FOR_EXCEPTION(
463 numVecs != Y.getNumVectors (), std::runtime_error,
464 "Ifpack2::Details::TriDiSolver::apply: X and Y have different numbers " 465 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 466 << X.getNumVectors () <<
".");
473 RCP<const MV> X_local;
475 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
480 X_local = X.offsetView (A_local_->getDomainMap (), 0);
481 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
485 X_local = rcpFromRef (X);
486 Y_local = rcpFromRef (Y);
491 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
496 applyTime_ += (timer->wallTime() - startTime);
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>
532 const Teuchos::EVerbosityLevel verbLevel)
const {
533 using Teuchos::FancyOStream;
534 using Teuchos::OSTab;
536 using Teuchos::rcpFromRef;
539 if (verbLevel == Teuchos::VERB_NONE) {
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;
556 if (verbLevel >= Teuchos::VERB_EXTREME) {
557 out <<
"A_local_tridi_:" << endl;
558 A_local_tridi_.print(out);
560 out <<
"ipiv_: " << Teuchos::toString (ipiv_) << endl;
564 template<
class MatrixType>
566 const Teuchos::EVerbosityLevel verbLevel)
const 568 using Teuchos::FancyOStream;
569 using Teuchos::OSTab;
571 using Teuchos::rcpFromRef;
574 RCP<FancyOStream> ptrOut = rcpFromRef (out);
581 if (verbLevel > Teuchos::VERB_NONE) {
582 out <<
"Ifpack2::Details::TriDiSolver:" << endl;
584 describeLocal (out, verbLevel);
589 const Teuchos::Comm<int>& comm = * (A_->getRowMap ()->getComm ());
590 const int myRank = comm.getRank ();
591 const int numProcs = comm.getSize ();
592 if (verbLevel > Teuchos::VERB_NONE && myRank == 0) {
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)
612 using Teuchos::Array;
613 using Teuchos::ArrayView;
614 typedef local_ordinal_type LO;
615 typedef typename Teuchos::ArrayView<LO>::size_type size_type;
618 A_local_tridi.putScalar (STS::zero ());
626 const map_type& rowMap = * (A_local.getRowMap ());
630 const size_type maxNumRowEntries =
631 static_cast<size_type
> (A_local.getNodeMaxNumRowEntries ());
632 nonconst_local_inds_host_view_type localIndices(
"localIndices",maxNumRowEntries);
633 nonconst_values_host_view_type values (
"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,
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>
671 TriDiSolver<MatrixType, true>::TriDiSolver (
const Teuchos::RCP<const row_matrix_type>& A) {
672 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
676 template<
class MatrixType>
677 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
679 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
683 template<
class MatrixType>
684 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::map_type>
686 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
690 template<
class MatrixType>
693 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
697 template<
class MatrixType>
700 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
704 template<
class MatrixType>
707 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
711 template<
class MatrixType>
714 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
718 template<
class MatrixType>
721 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
725 template<
class MatrixType>
728 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
732 template<
class MatrixType>
735 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
739 template<
class MatrixType>
742 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
746 template<
class MatrixType>
749 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
753 template<
class MatrixType>
754 Teuchos::RCP<const typename TriDiSolver<MatrixType, true>::row_matrix_type>
756 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
760 template<
class MatrixType>
763 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
767 template<
class MatrixType>
770 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
774 template<
class MatrixType>
775 TriDiSolver<MatrixType, true>::~TriDiSolver ()
781 template<
class MatrixType>
784 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
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,
792 Teuchos::ETransp mode,
794 scalar_type beta)
const 796 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
800 template<
class MatrixType>
802 TriDiSolver<MatrixType, true>::description ()
const 804 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
808 template<
class MatrixType>
809 void TriDiSolver<MatrixType, true>::describe(Teuchos::FancyOStream& out,
810 const Teuchos::EVerbosityLevel verbLevel)
const 812 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
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 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:76
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.
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: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.
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:73
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.