43 #ifndef IFPACK2_DETAILS_DENSESOLVER_DEF_HPP 44 #define IFPACK2_DETAILS_DENSESOLVER_DEF_HPP 46 #include "Ifpack2_LocalFilter.hpp" 47 #include "Teuchos_LAPACK.hpp" 48 #include "Ifpack2_Details_DenseSolver.hpp" 52 # include "Teuchos_DefaultMpiComm.hpp" 54 # include "Teuchos_DefaultSerialComm.hpp" 65 template<
class MatrixType>
69 initializeTime_ (0.0),
75 isInitialized_ (false),
80 template<
class MatrixType>
81 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::map_type>
83 TEUCHOS_TEST_FOR_EXCEPTION(
84 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 85 "getDomainMap: The input matrix A is null. Please call setMatrix() with a " 86 "nonnull input matrix before calling this method.");
89 return A_->getRangeMap ();
93 template<
class MatrixType>
94 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::map_type>
96 TEUCHOS_TEST_FOR_EXCEPTION(
97 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 98 "getRangeMap: The input matrix A is null. Please call setMatrix() with a " 99 "nonnull input matrix before calling this method.");
102 return A_->getDomainMap ();
106 template<
class MatrixType>
114 template<
class MatrixType>
117 return isInitialized_;
121 template<
class MatrixType>
128 template<
class MatrixType>
131 return numInitialize_;
135 template<
class MatrixType>
142 template<
class MatrixType>
149 template<
class MatrixType>
152 return initializeTime_;
156 template<
class MatrixType>
163 template<
class MatrixType>
170 template<
class MatrixType>
171 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::row_matrix_type>
177 template<
class MatrixType>
181 isInitialized_ =
false;
183 A_local_ = Teuchos::null;
184 A_local_dense_.reshape (0, 0);
189 template<
class MatrixType>
191 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
195 if (! A_.is_null ()) {
196 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
197 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
198 TEUCHOS_TEST_FOR_EXCEPTION(
199 numRows != numCols, std::invalid_argument,
"Ifpack2::Details::DenseSolver::" 200 "setMatrix: Input matrix must be (globally) square. " 201 "The matrix you provided is " << numRows <<
" by " << numCols <<
".");
211 template<
class MatrixType>
219 using Teuchos::TimeMonitor;
220 const std::string timerName (
"Ifpack2::Details::DenseSolver::initialize");
222 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
223 if (timer.is_null ()) {
224 timer = TimeMonitor::getNewCounter (timerName);
228 Teuchos::TimeMonitor timeMon (*timer);
230 TEUCHOS_TEST_FOR_EXCEPTION(
231 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 232 "initialize: The input matrix A is null. Please call setMatrix() " 233 "with a nonnull input before calling this method.");
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 ! A_->hasColMap (), std::invalid_argument,
"Ifpack2::Details::DenseSolver: " 237 "The constructor's input matrix must have a column Map, " 238 "so that it has local indices.");
244 if (A_->getComm ()->getSize () > 1) {
250 TEUCHOS_TEST_FOR_EXCEPTION(
251 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::DenseSolver::" 252 "initialize: A_local_ is null after it was supposed to have been " 253 "initialized. Please report this bug to the Ifpack2 developers.");
256 const size_t numRows = A_local_->getNodeNumRows ();
257 const size_t numCols = A_local_->getNodeNumCols ();
258 TEUCHOS_TEST_FOR_EXCEPTION(
259 numRows != numCols, std::logic_error,
"Ifpack2::Details::DenseSolver::" 260 "initialize: Local filter matrix is not square. This should never happen. " 261 "Please report this bug to the Ifpack2 developers.");
262 A_local_dense_.reshape (numRows, numCols);
263 ipiv_.resize (std::min (numRows, numCols));
264 std::fill (ipiv_.begin (), ipiv_.end (), 0);
266 isInitialized_ =
true;
272 initializeTime_ = timer->totalElapsedTime ();
276 template<
class MatrixType>
281 template<
class MatrixType>
285 const std::string timerName (
"Ifpack2::Details::DenseSolver::compute");
287 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
288 if (timer.is_null ()) {
289 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
294 Teuchos::TimeMonitor timeMon (*timer);
295 TEUCHOS_TEST_FOR_EXCEPTION(
296 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 297 "compute: The input matrix A is null. Please call setMatrix() with a " 298 "nonnull input, then call initialize(), before calling this method.");
300 TEUCHOS_TEST_FOR_EXCEPTION(
301 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::DenseSolver::" 302 "compute: A_local_ is null. Please report this bug to the Ifpack2 " 309 extract (A_local_dense_, *A_local_);
310 factor (A_local_dense_, ipiv_ ());
317 computeTime_ = timer->totalElapsedTime ();
320 template<
class MatrixType>
322 factor (Teuchos::SerialDenseMatrix<int, scalar_type>& A,
323 const Teuchos::ArrayView<int>& ipiv)
326 std::fill (ipiv.begin (), ipiv.end (), 0);
328 Teuchos::LAPACK<int, scalar_type> lapack;
330 lapack.GETRF (A.numRows (), A.numCols (), A.values (), A.stride (),
331 ipiv.getRawPtr (), &INFO);
333 TEUCHOS_TEST_FOR_EXCEPTION(
334 INFO < 0, std::logic_error,
"Ifpack2::Details::DenseSolver::factor: " 335 "LAPACK's _GETRF (LU factorization with partial pivoting) was called " 336 "incorrectly. INFO = " << INFO <<
" < 0. " 337 "Please report this bug to the Ifpack2 developers.");
341 TEUCHOS_TEST_FOR_EXCEPTION(
342 INFO > 0, std::runtime_error,
"Ifpack2::Details::DenseSolver::factor: " 343 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the " 344 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") " 345 "(one-based index i) is exactly zero. This probably means that the input " 346 "matrix has a singular diagonal block.");
350 template<
class MatrixType>
351 void DenseSolver<MatrixType, false>::
352 applyImpl (
const MV& X,
354 const Teuchos::ETransp mode,
355 const scalar_type alpha,
356 const scalar_type beta)
const 358 using Teuchos::ArrayRCP;
361 using Teuchos::rcpFromRef;
362 using Teuchos::CONJ_TRANS;
363 using Teuchos::TRANS;
365 const int numVecs =
static_cast<int> (X.getNumVectors ());
366 if (alpha == STS::zero ()) {
367 if (beta == STS::zero ()) {
371 Y.putScalar (STS::zero ());
374 Y.scale (STS::zero ());
378 Teuchos::LAPACK<int, scalar_type> lapack;
384 if (beta == STS::zero () && Y.isConstantStride () && alpha == STS::one ()) {
386 Y_tmp = rcpFromRef (Y);
389 Y_tmp = rcp (
new MV (X, Teuchos::Copy));
390 if (alpha != STS::one ()) {
391 Y_tmp->scale (alpha);
394 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
395 ArrayRCP<scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
396 scalar_type*
const Y_ptr = Y_view.getRawPtr ();
399 (mode == CONJ_TRANS ?
'C' : (mode == TRANS ?
'T' :
'N'));
400 lapack.GETRS (trans, A_local_dense_.numRows (), numVecs,
401 A_local_dense_.values (), A_local_dense_.stride (),
402 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
403 TEUCHOS_TEST_FOR_EXCEPTION(
404 INFO != 0, std::runtime_error,
"Ifpack2::Details::DenseSolver::" 405 "applyImpl: LAPACK's _GETRS (solve using LU factorization with " 406 "partial pivoting) failed with INFO = " << INFO <<
" != 0.");
408 if (beta != STS::zero ()) {
409 Y.update (alpha, *Y_tmp, beta);
411 else if (! Y.isConstantStride ()) {
412 deep_copy (Y, *Y_tmp);
418 template<
class MatrixType>
420 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
421 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
422 Teuchos::ETransp mode,
426 using Teuchos::ArrayView;
430 using Teuchos::rcpFromRef;
432 const std::string timerName (
"Ifpack2::Details::DenseSolver::apply");
433 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
434 if (timer.is_null ()) {
435 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
440 Teuchos::TimeMonitor timeMon (*timer);
442 TEUCHOS_TEST_FOR_EXCEPTION(
443 ! isComputed_, std::runtime_error,
"Ifpack2::Details::DenseSolver::apply: " 444 "You must have called the compute() method before you may call apply(). " 445 "You may call the apply() method as many times as you want after calling " 446 "compute() once, but you must have called compute() at least once.");
448 const size_t numVecs = X.getNumVectors ();
450 TEUCHOS_TEST_FOR_EXCEPTION(
451 numVecs != Y.getNumVectors (), std::runtime_error,
452 "Ifpack2::Details::DenseSolver::apply: X and Y have different numbers " 453 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 454 << X.getNumVectors () <<
".");
461 RCP<const MV> X_local;
463 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
468 X_local = X.offsetView (A_local_->getDomainMap (), 0);
469 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
473 X_local = rcpFromRef (X);
474 Y_local = rcpFromRef (Y);
479 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
486 applyTime_ = timer->totalElapsedTime ();
490 template<
class MatrixType>
494 std::ostringstream out;
499 out <<
"\"Ifpack2::Details::DenseSolver\": ";
501 if (this->getObjectLabel () !=
"") {
502 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
504 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", " 505 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
508 out <<
"Matrix: null";
511 out <<
"Matrix: not null" 512 <<
", Global matrix dimensions: [" 513 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
521 template<
class MatrixType>
524 const Teuchos::EVerbosityLevel verbLevel)
const 526 using Teuchos::FancyOStream;
527 using Teuchos::OSTab;
529 using Teuchos::rcpFromRef;
532 if (verbLevel == Teuchos::VERB_NONE) {
536 RCP<FancyOStream> ptrOut = rcpFromRef (out);
538 if (this->getObjectLabel () !=
"") {
539 out <<
"label: " << this->getObjectLabel () << endl;
541 out <<
"initialized: " << (isInitialized_ ?
"true" :
"false") << endl
542 <<
"computed: " << (isComputed_ ?
"true" :
"false") << endl
543 <<
"number of initialize calls: " << numInitialize_ << endl
544 <<
"number of compute calls: " << numCompute_ << endl
545 <<
"number of apply calls: " << numApply_ << endl
546 <<
"total time in seconds in initialize: " << initializeTime_ << endl
547 <<
"total time in seconds in compute: " << computeTime_ << endl
548 <<
"total time in seconds in apply: " << applyTime_ << endl;
549 if (verbLevel >= Teuchos::VERB_EXTREME) {
550 out <<
"A_local_dense_:" << endl;
554 for (
int i = 0; i < A_local_dense_.numRows (); ++i) {
555 for (
int j = 0; j < A_local_dense_.numCols (); ++j) {
556 out << A_local_dense_(i,j);
557 if (j + 1 < A_local_dense_.numCols ()) {
561 if (i + 1 < A_local_dense_.numRows ()) {
567 out <<
"ipiv_: " << Teuchos::toString (ipiv_) << endl;
573 template<
class MatrixType>
576 const Teuchos::EVerbosityLevel verbLevel)
const 578 using Teuchos::FancyOStream;
579 using Teuchos::OSTab;
581 using Teuchos::rcpFromRef;
584 RCP<FancyOStream> ptrOut = rcpFromRef (out);
591 if (verbLevel > Teuchos::VERB_NONE) {
592 out <<
"Ifpack2::Details::DenseSolver:" << endl;
594 describeLocal (out, verbLevel);
599 const Teuchos::Comm<int>& comm = * (A_->getRowMap ()->getComm ());
600 const int myRank = comm.getRank ();
601 const int numProcs = comm.getSize ();
602 if (verbLevel > Teuchos::VERB_NONE && myRank == 0) {
603 out <<
"Ifpack2::Details::DenseSolver:" << endl;
606 for (
int p = 0; p < numProcs; ++p) {
608 out <<
"Process " << myRank <<
":" << endl;
609 describeLocal (out, verbLevel);
619 template<
class MatrixType>
621 extract (Teuchos::SerialDenseMatrix<int, scalar_type>& A_local_dense,
622 const row_matrix_type& A_local)
624 using Teuchos::Array;
625 using Teuchos::ArrayView;
626 typedef local_ordinal_type LO;
627 typedef typename Teuchos::ArrayView<LO>::size_type size_type;
630 A_local_dense.putScalar (STS::zero ());
638 const map_type& rowMap = * (A_local.getRowMap ());
642 const size_type maxNumRowEntries =
643 static_cast<size_type
> (A_local.getNodeMaxNumRowEntries ());
644 Array<LO> localIndices (maxNumRowEntries);
645 Array<scalar_type> values (maxNumRowEntries);
647 const LO numLocalRows =
static_cast<LO
> (rowMap.getNodeNumElements ());
648 const LO minLocalRow = rowMap.getMinLocalIndex ();
651 const LO maxLocalRow = minLocalRow + numLocalRows;
652 for (LO localRow = minLocalRow; localRow < maxLocalRow; ++localRow) {
659 const size_type numEntriesInRow =
660 static_cast<size_type
> (A_local.getNumEntriesInLocalRow (localRow));
661 size_t numEntriesOut = 0;
662 A_local.getLocalRowCopy (localRow,
663 localIndices (0, numEntriesInRow),
664 values (0, numEntriesInRow),
666 for (LO k = 0; k < numEntriesInRow; ++k) {
667 const LO localCol = localIndices[k];
668 const scalar_type val = values[k];
671 A_local_dense(localRow, localCol) += val;
680 template<
class MatrixType>
681 DenseSolver<MatrixType, true>::
682 DenseSolver (
const Teuchos::RCP<const row_matrix_type>& A) {
683 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
687 template<
class MatrixType>
688 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::map_type>
690 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
694 template<
class MatrixType>
695 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::map_type>
697 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
701 template<
class MatrixType>
705 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
709 template<
class MatrixType>
712 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
716 template<
class MatrixType>
719 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
723 template<
class MatrixType>
726 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
730 template<
class MatrixType>
733 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
737 template<
class MatrixType>
740 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
744 template<
class MatrixType>
747 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
751 template<
class MatrixType>
754 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
758 template<
class MatrixType>
761 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
765 template<
class MatrixType>
766 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::row_matrix_type>
768 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
772 template<
class MatrixType>
774 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
776 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
780 template<
class MatrixType>
783 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
787 template<
class MatrixType>
788 DenseSolver<MatrixType, true>::~DenseSolver ()
794 template<
class MatrixType>
797 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
801 template<
class MatrixType>
803 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
804 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
805 Teuchos::ETransp mode,
807 scalar_type beta)
const 809 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
813 template<
class MatrixType>
815 DenseSolver<MatrixType, true>::description ()
const 817 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
821 template<
class MatrixType>
822 void DenseSolver<MatrixType, true>::
823 describe (Teuchos::FancyOStream& out,
824 const Teuchos::EVerbosityLevel verbLevel)
const 826 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
833 #define IFPACK2_DETAILS_DENSESOLVER_INSTANT(S,LO,GO,N) \ 834 template class Ifpack2::Details::DenseSolver< Tpetra::RowMatrix<S, LO, GO, N> >; 836 #endif // IFPACK2_DETAILS_DENSESOLVER_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.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to the given FancyOStream.
Definition: Ifpack2_Details_DenseSolver_def.hpp:575
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.
DenseSolver(const Teuchos::RCP< const row_matrix_type > &matrix)
Constructor.
Definition: Ifpack2_Details_DenseSolver_def.hpp:67
virtual void compute()=0
Set up the numerical values in this preconditioner.
"Preconditioner" that uses LAPACK's dense LU.
Definition: Ifpack2_Details_DenseSolver_decl.hpp:73
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().
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.
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_Details_DenseSolver_decl.hpp:106
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.