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" 49 #include "Tpetra_Map.hpp" 53 # include "Teuchos_DefaultMpiComm.hpp" 55 # include "Teuchos_DefaultSerialComm.hpp" 66 template<
class MatrixType>
70 initializeTime_ (0.0),
76 isInitialized_ (false),
81 template<
class MatrixType>
82 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::map_type>
84 TEUCHOS_TEST_FOR_EXCEPTION(
85 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 86 "getDomainMap: The input matrix A is null. Please call setMatrix() with a " 87 "nonnull input matrix before calling this method.");
90 return A_->getRangeMap ();
94 template<
class MatrixType>
95 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::map_type>
97 TEUCHOS_TEST_FOR_EXCEPTION(
98 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 99 "getRangeMap: The input matrix A is null. Please call setMatrix() with a " 100 "nonnull input matrix before calling this method.");
103 return A_->getDomainMap ();
107 template<
class MatrixType>
115 template<
class MatrixType>
118 return isInitialized_;
122 template<
class MatrixType>
129 template<
class MatrixType>
132 return numInitialize_;
136 template<
class MatrixType>
143 template<
class MatrixType>
150 template<
class MatrixType>
153 return initializeTime_;
157 template<
class MatrixType>
164 template<
class MatrixType>
171 template<
class MatrixType>
172 Teuchos::RCP<const typename DenseSolver<MatrixType, false>::row_matrix_type>
178 template<
class MatrixType>
182 isInitialized_ =
false;
184 A_local_ = Teuchos::null;
185 A_local_dense_.reshape (0, 0);
190 template<
class MatrixType>
192 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
196 if (! A_.is_null ()) {
197 const global_size_t numRows = A->getRangeMap ()->getGlobalNumElements ();
198 const global_size_t numCols = A->getDomainMap ()->getGlobalNumElements ();
199 TEUCHOS_TEST_FOR_EXCEPTION(
200 numRows != numCols, std::invalid_argument,
"Ifpack2::Details::DenseSolver::" 201 "setMatrix: Input matrix must be (globally) square. " 202 "The matrix you provided is " << numRows <<
" by " << numCols <<
".");
212 template<
class MatrixType>
220 using Teuchos::TimeMonitor;
221 const std::string timerName (
"Ifpack2::Details::DenseSolver::initialize");
223 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
224 if (timer.is_null ()) {
225 timer = TimeMonitor::getNewCounter (timerName);
229 Teuchos::TimeMonitor timeMon (*timer);
231 TEUCHOS_TEST_FOR_EXCEPTION(
232 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 233 "initialize: The input matrix A is null. Please call setMatrix() " 234 "with a nonnull input before calling this method.");
236 TEUCHOS_TEST_FOR_EXCEPTION(
237 ! A_->hasColMap (), std::invalid_argument,
"Ifpack2::Details::DenseSolver: " 238 "The constructor's input matrix must have a column Map, " 239 "so that it has local indices.");
245 if (A_->getComm ()->getSize () > 1) {
251 TEUCHOS_TEST_FOR_EXCEPTION(
252 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::DenseSolver::" 253 "initialize: A_local_ is null after it was supposed to have been " 254 "initialized. Please report this bug to the Ifpack2 developers.");
257 const size_t numRows = A_local_->getNodeNumRows ();
258 const size_t numCols = A_local_->getNodeNumCols ();
259 TEUCHOS_TEST_FOR_EXCEPTION(
260 numRows != numCols, std::logic_error,
"Ifpack2::Details::DenseSolver::" 261 "initialize: Local filter matrix is not square. This should never happen. " 262 "Please report this bug to the Ifpack2 developers.");
263 A_local_dense_.reshape (numRows, numCols);
264 ipiv_.resize (std::min (numRows, numCols));
265 std::fill (ipiv_.begin (), ipiv_.end (), 0);
267 isInitialized_ =
true;
273 initializeTime_ = timer->totalElapsedTime ();
277 template<
class MatrixType>
282 template<
class MatrixType>
286 const std::string timerName (
"Ifpack2::Details::DenseSolver::compute");
288 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
289 if (timer.is_null ()) {
290 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
295 Teuchos::TimeMonitor timeMon (*timer);
296 TEUCHOS_TEST_FOR_EXCEPTION(
297 A_.is_null (), std::runtime_error,
"Ifpack2::Details::DenseSolver::" 298 "compute: The input matrix A is null. Please call setMatrix() with a " 299 "nonnull input, then call initialize(), before calling this method.");
301 TEUCHOS_TEST_FOR_EXCEPTION(
302 A_local_.is_null (), std::logic_error,
"Ifpack2::Details::DenseSolver::" 303 "compute: A_local_ is null. Please report this bug to the Ifpack2 " 310 extract (A_local_dense_, *A_local_);
311 factor (A_local_dense_, ipiv_ ());
318 computeTime_ = timer->totalElapsedTime ();
321 template<
class MatrixType>
323 factor (Teuchos::SerialDenseMatrix<int, scalar_type>& A,
324 const Teuchos::ArrayView<int>& ipiv)
327 std::fill (ipiv.begin (), ipiv.end (), 0);
329 Teuchos::LAPACK<int, scalar_type> lapack;
331 lapack.GETRF (A.numRows (), A.numCols (), A.values (), A.stride (),
332 ipiv.getRawPtr (), &INFO);
334 TEUCHOS_TEST_FOR_EXCEPTION(
335 INFO < 0, std::logic_error,
"Ifpack2::Details::DenseSolver::factor: " 336 "LAPACK's _GETRF (LU factorization with partial pivoting) was called " 337 "incorrectly. INFO = " << INFO <<
" < 0. " 338 "Please report this bug to the Ifpack2 developers.");
342 TEUCHOS_TEST_FOR_EXCEPTION(
343 INFO > 0, std::runtime_error,
"Ifpack2::Details::DenseSolver::factor: " 344 "LAPACK's _GETRF (LU factorization with partial pivoting) reports that the " 345 "computed U factor is exactly singular. U(" << INFO <<
"," << INFO <<
") " 346 "(one-based index i) is exactly zero. This probably means that the input " 347 "matrix has a singular diagonal block.");
351 template<
class MatrixType>
352 void DenseSolver<MatrixType, false>::
353 applyImpl (
const MV& X,
355 const Teuchos::ETransp mode,
356 const scalar_type alpha,
357 const scalar_type beta)
const 359 using Teuchos::ArrayRCP;
362 using Teuchos::rcpFromRef;
363 using Teuchos::CONJ_TRANS;
364 using Teuchos::TRANS;
366 const int numVecs =
static_cast<int> (X.getNumVectors ());
367 if (alpha == STS::zero ()) {
368 if (beta == STS::zero ()) {
372 Y.putScalar (STS::zero ());
375 Y.scale (STS::zero ());
379 Teuchos::LAPACK<int, scalar_type> lapack;
385 if (beta == STS::zero () && Y.isConstantStride () && alpha == STS::one ()) {
387 Y_tmp = rcpFromRef (Y);
390 Y_tmp = rcp (
new MV (X, Teuchos::Copy));
391 if (alpha != STS::one ()) {
392 Y_tmp->scale (alpha);
395 const int Y_stride =
static_cast<int> (Y_tmp->getStride ());
396 ArrayRCP<scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
397 scalar_type*
const Y_ptr = Y_view.getRawPtr ();
400 (mode == CONJ_TRANS ?
'C' : (mode ==
TRANS ?
'T' :
'N'));
401 lapack.GETRS (trans, A_local_dense_.numRows (), numVecs,
402 A_local_dense_.values (), A_local_dense_.stride (),
403 ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
404 TEUCHOS_TEST_FOR_EXCEPTION(
405 INFO != 0, std::runtime_error,
"Ifpack2::Details::DenseSolver::" 406 "applyImpl: LAPACK's _GETRS (solve using LU factorization with " 407 "partial pivoting) failed with INFO = " << INFO <<
" != 0.");
409 if (beta != STS::zero ()) {
410 Y.update (alpha, *Y_tmp, beta);
412 else if (! Y.isConstantStride ()) {
413 deep_copy (Y, *Y_tmp);
419 template<
class MatrixType>
421 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
422 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
423 Teuchos::ETransp mode,
427 using Teuchos::ArrayView;
431 using Teuchos::rcpFromRef;
433 const std::string timerName (
"Ifpack2::Details::DenseSolver::apply");
434 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
435 if (timer.is_null ()) {
436 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
441 Teuchos::TimeMonitor timeMon (*timer);
443 TEUCHOS_TEST_FOR_EXCEPTION(
444 ! isComputed_, std::runtime_error,
"Ifpack2::Details::DenseSolver::apply: " 445 "You must have called the compute() method before you may call apply(). " 446 "You may call the apply() method as many times as you want after calling " 447 "compute() once, but you must have called compute() at least once.");
449 const size_t numVecs = X.getNumVectors ();
451 TEUCHOS_TEST_FOR_EXCEPTION(
452 numVecs != Y.getNumVectors (), std::runtime_error,
453 "Ifpack2::Details::DenseSolver::apply: X and Y have different numbers " 454 "of vectors. X has " << X.getNumVectors () <<
", but Y has " 455 << X.getNumVectors () <<
".");
462 RCP<const MV> X_local;
464 const bool multipleProcs = (A_->getRowMap ()->getComm ()->getSize () >= 1);
469 X_local = X.offsetView (A_local_->getDomainMap (), 0);
470 Y_local = Y.offsetViewNonConst (A_local_->getRangeMap (), 0);
474 X_local = rcpFromRef (X);
475 Y_local = rcpFromRef (Y);
480 this->applyImpl (*X_local, *Y_local, mode, alpha, beta);
487 applyTime_ = timer->totalElapsedTime ();
491 template<
class MatrixType>
495 std::ostringstream out;
500 out <<
"\"Ifpack2::Details::DenseSolver\": ";
502 if (this->getObjectLabel () !=
"") {
503 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
505 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", " 506 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
509 out <<
"Matrix: null";
512 out <<
"Matrix: not null" 513 <<
", Global matrix dimensions: [" 514 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
522 template<
class MatrixType>
525 const Teuchos::EVerbosityLevel verbLevel)
const 527 using Teuchos::FancyOStream;
528 using Teuchos::OSTab;
530 using Teuchos::rcpFromRef;
533 if (verbLevel == Teuchos::VERB_NONE) {
537 RCP<FancyOStream> ptrOut = rcpFromRef (out);
539 if (this->getObjectLabel () !=
"") {
540 out <<
"label: " << this->getObjectLabel () << endl;
542 out <<
"initialized: " << (isInitialized_ ?
"true" :
"false") << endl
543 <<
"computed: " << (isComputed_ ?
"true" :
"false") << endl
544 <<
"number of initialize calls: " << numInitialize_ << endl
545 <<
"number of compute calls: " << numCompute_ << endl
546 <<
"number of apply calls: " << numApply_ << endl
547 <<
"total time in seconds in initialize: " << initializeTime_ << endl
548 <<
"total time in seconds in compute: " << computeTime_ << endl
549 <<
"total time in seconds in apply: " << applyTime_ << endl;
550 if (verbLevel >= Teuchos::VERB_EXTREME) {
551 out <<
"A_local_dense_:" << endl;
555 for (
int i = 0; i < A_local_dense_.numRows (); ++i) {
556 for (
int j = 0; j < A_local_dense_.numCols (); ++j) {
557 out << A_local_dense_(i,j);
558 if (j + 1 < A_local_dense_.numCols ()) {
562 if (i + 1 < A_local_dense_.numRows ()) {
568 out <<
"ipiv_: " << Teuchos::toString (ipiv_) << endl;
574 template<
class MatrixType>
577 const Teuchos::EVerbosityLevel verbLevel)
const 579 using Teuchos::FancyOStream;
580 using Teuchos::OSTab;
582 using Teuchos::rcpFromRef;
585 RCP<FancyOStream> ptrOut = rcpFromRef (out);
592 if (verbLevel > Teuchos::VERB_NONE) {
593 out <<
"Ifpack2::Details::DenseSolver:" << endl;
595 describeLocal (out, verbLevel);
600 const Teuchos::Comm<int>& comm = * (A_->getRowMap ()->getComm ());
601 const int myRank = comm.getRank ();
602 const int numProcs = comm.getSize ();
603 if (verbLevel > Teuchos::VERB_NONE && myRank == 0) {
604 out <<
"Ifpack2::Details::DenseSolver:" << endl;
607 for (
int p = 0; p < numProcs; ++p) {
609 out <<
"Process " << myRank <<
":" << endl;
610 describeLocal (out, verbLevel);
620 template<
class MatrixType>
622 extract (Teuchos::SerialDenseMatrix<int, scalar_type>& A_local_dense,
623 const row_matrix_type& A_local)
625 using Teuchos::Array;
626 using Teuchos::ArrayView;
627 typedef local_ordinal_type LO;
628 typedef typename Teuchos::ArrayView<LO>::size_type size_type;
631 A_local_dense.putScalar (STS::zero ());
639 const map_type& rowMap = * (A_local.getRowMap ());
643 const size_type maxNumRowEntries =
644 static_cast<size_type
> (A_local.getNodeMaxNumRowEntries ());
645 Array<LO> localIndices (maxNumRowEntries);
646 Array<scalar_type> values (maxNumRowEntries);
648 const LO numLocalRows =
static_cast<LO
> (rowMap.getNodeNumElements ());
649 const LO minLocalRow = rowMap.getMinLocalIndex ();
652 const LO maxLocalRow = minLocalRow + numLocalRows;
653 for (LO localRow = minLocalRow; localRow < maxLocalRow; ++localRow) {
660 const size_type numEntriesInRow =
661 static_cast<size_type
> (A_local.getNumEntriesInLocalRow (localRow));
662 size_t numEntriesOut = 0;
663 A_local.getLocalRowCopy (localRow,
664 localIndices (0, numEntriesInRow),
665 values (0, numEntriesInRow),
667 for (LO k = 0; k < numEntriesInRow; ++k) {
668 const LO localCol = localIndices[k];
669 const scalar_type val = values[k];
672 A_local_dense(localRow, localCol) += val;
681 template<
class MatrixType>
682 DenseSolver<MatrixType, true>::
683 DenseSolver (
const Teuchos::RCP<const row_matrix_type>& A) {
684 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
688 template<
class MatrixType>
689 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::map_type>
691 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
695 template<
class MatrixType>
696 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::map_type>
698 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
702 template<
class MatrixType>
706 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
710 template<
class MatrixType>
713 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
717 template<
class MatrixType>
720 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
724 template<
class MatrixType>
727 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
731 template<
class MatrixType>
734 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
738 template<
class MatrixType>
741 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
745 template<
class MatrixType>
748 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
752 template<
class MatrixType>
755 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
759 template<
class MatrixType>
762 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
766 template<
class MatrixType>
767 Teuchos::RCP<const typename DenseSolver<MatrixType, true>::row_matrix_type>
769 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
773 template<
class MatrixType>
775 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
777 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
781 template<
class MatrixType>
784 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
788 template<
class MatrixType>
789 DenseSolver<MatrixType, true>::~DenseSolver ()
795 template<
class MatrixType>
798 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
802 template<
class MatrixType>
804 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
805 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
806 Teuchos::ETransp mode,
808 scalar_type beta)
const 810 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
814 template<
class MatrixType>
816 DenseSolver<MatrixType, true>::description ()
const 818 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
822 template<
class MatrixType>
823 void DenseSolver<MatrixType, true>::
824 describe (Teuchos::FancyOStream& out,
825 const Teuchos::EVerbosityLevel verbLevel)
const 827 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented");
834 #define IFPACK2_DETAILS_DENSESOLVER_INSTANT(S,LO,GO,N) \ 835 template class Ifpack2::Details::DenseSolver< Tpetra::RowMatrix<S, LO, GO, N> >; 837 #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:576
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:68
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:75
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: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.
MatrixType::scalar_type scalar_type
The type of entries in the input (global) matrix.
Definition: Ifpack2_Details_DenseSolver_decl.hpp:108
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.