43 #ifndef BELOS_LINEAR_PROBLEM_HPP 44 #define BELOS_LINEAR_PROBLEM_HPP 51 #include "Teuchos_ParameterList.hpp" 52 #include "Teuchos_TimeMonitor.hpp" 81 template <
class ScalarType,
class MV,
class OP>
104 const Teuchos::RCP<MV> &X,
105 const Teuchos::RCP<const MV> &B);
134 void setLHS (
const Teuchos::RCP<MV> &X) {
143 void setRHS (
const Teuchos::RCP<const MV> &B) {
177 void setLSIndex (
const std::vector<int>& index);
197 void setLabel (
const std::string& label);
236 updateSolution (
const Teuchos::RCP<MV>& update = Teuchos::null,
237 bool updateLP =
false,
238 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
258 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() )
const 292 setProblem (
const Teuchos::RCP<MV> &newX = Teuchos::null,
293 const Teuchos::RCP<const MV> &newB = Teuchos::null);
304 Teuchos::RCP<MV>
getLHS()
const {
return(X_); }
307 Teuchos::RCP<const MV>
getRHS()
const {
return(B_); }
332 Teuchos::RCP<MV> getCurrLHSVec();
348 Teuchos::RCP<const MV> getCurrRHSVec();
377 const std::vector<int>
getLSIndex()
const {
return(rhsIndex_); }
393 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
394 return Teuchos::tuple(timerOp_,timerPrec_);
441 void apply(
const MV& x, MV& y )
const;
450 void applyOp(
const MV& x, MV& y )
const;
458 void applyLeftPrec(
const MV& x, MV& y )
const;
466 void applyRightPrec(
const MV& x, MV& y )
const;
473 void computeCurrResVec( MV* R ,
const MV* X = 0,
const MV* B = 0 )
const;
480 void computeCurrPrecResVec( MV* R,
const MV* X = 0,
const MV* B = 0 )
const;
487 Teuchos::RCP<const OP> A_;
493 Teuchos::RCP<MV> curX_;
496 Teuchos::RCP<const MV> B_;
499 Teuchos::RCP<const MV> curB_;
502 Teuchos::RCP<MV> R0_;
505 Teuchos::RCP<MV> PR0_;
508 Teuchos::RCP<const OP> LP_;
511 Teuchos::RCP<const OP> RP_;
514 mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
523 std::vector<int> rhsIndex_;
545 bool solutionUpdated_;
560 template <
class ScalarType,
class MV,
class OP>
570 solutionUpdated_(false),
575 template <
class ScalarType,
class MV,
class OP>
577 const Teuchos::RCP<MV> &X,
578 const Teuchos::RCP<const MV> &B
591 solutionUpdated_(false),
596 template <
class ScalarType,
class MV,
class OP>
600 curX_(Problem.curX_),
602 curB_(Problem.curB_),
607 timerOp_(Problem.timerOp_),
608 timerPrec_(Problem.timerPrec_),
609 blocksize_(Problem.blocksize_),
610 num2Solve_(Problem.num2Solve_),
611 rhsIndex_(Problem.rhsIndex_),
612 lsNum_(Problem.lsNum_),
613 Left_Scale_(Problem.Left_Scale_),
614 Right_Scale_(Problem.Right_Scale_),
615 isSet_(Problem.isSet_),
616 isHermitian_(Problem.isHermitian_),
617 solutionUpdated_(Problem.solutionUpdated_),
618 label_(Problem.label_)
622 template <
class ScalarType,
class MV,
class OP>
626 template <
class ScalarType,
class MV,
class OP>
634 curB_ = Teuchos::null;
635 curX_ = Teuchos::null;
638 int validIdx = 0, ivalidIdx = 0;
639 blocksize_ = rhsIndex_.size();
640 std::vector<int> vldIndex( blocksize_ );
641 std::vector<int> newIndex( blocksize_ );
642 std::vector<int> iIndex( blocksize_ );
643 for (
int i=0; i<blocksize_; ++i) {
644 if (rhsIndex_[i] > -1) {
645 vldIndex[validIdx] = rhsIndex_[i];
646 newIndex[validIdx] = i;
650 iIndex[ivalidIdx] = i;
654 vldIndex.resize(validIdx);
655 newIndex.resize(validIdx);
656 iIndex.resize(ivalidIdx);
657 num2Solve_ = validIdx;
660 if (num2Solve_ != blocksize_) {
661 newIndex.resize(num2Solve_);
662 vldIndex.resize(num2Solve_);
668 Teuchos::RCP<MV> tmpCurB =
MVT::Clone( *B_, blocksize_ );
680 solutionUpdated_ =
false;
693 template <
class ScalarType,
class MV,
class OP>
700 if (num2Solve_ < blocksize_) {
705 std::vector<int> newIndex( num2Solve_ );
706 std::vector<int> vldIndex( num2Solve_ );
707 for (
int i=0; i<blocksize_; ++i) {
708 if ( rhsIndex_[i] > -1 ) {
709 vldIndex[validIdx] = rhsIndex_[i];
710 newIndex[validIdx] = i;
721 curX_ = Teuchos::null;
722 curB_ = Teuchos::null;
727 template <
class ScalarType,
class MV,
class OP>
738 if (update.is_null())
756 RCP<MV> rightPrecUpdate =
759 #ifdef BELOS_TEUCHOS_TIME_MONITOR 760 Teuchos::TimeMonitor PrecTimer (*timerPrec_);
762 OPT::Apply( *RP_, *update, *rightPrecUpdate );
765 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
767 solutionUpdated_ =
true;
783 RCP<MV> rightPrecUpdate =
786 #ifdef BELOS_TEUCHOS_TIME_MONITOR 787 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
789 OPT::Apply( *RP_, *update, *rightPrecUpdate );
792 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
799 template <
class ScalarType,
class MV,
class OP>
802 if (label != label_) {
805 if (timerOp_ != Teuchos::null) {
806 std::string opLabel = label_ +
": Operation Op*x";
807 #ifdef BELOS_TEUCHOS_TIME_MONITOR 808 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
811 if (timerPrec_ != Teuchos::null) {
812 std::string precLabel = label_ +
": Operation Prec*x";
813 #ifdef BELOS_TEUCHOS_TIME_MONITOR 814 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
820 template <
class ScalarType,
class MV,
class OP>
824 const Teuchos::RCP<const MV> &newB)
827 if (timerOp_ == Teuchos::null) {
828 std::string opLabel = label_ +
": Operation Op*x";
829 #ifdef BELOS_TEUCHOS_TIME_MONITOR 830 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
833 if (timerPrec_ == Teuchos::null) {
834 std::string precLabel = label_ +
": Operation Prec*x";
835 #ifdef BELOS_TEUCHOS_TIME_MONITOR 836 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
841 if (newX != Teuchos::null)
843 if (newB != Teuchos::null)
848 curX_ = Teuchos::null;
849 curB_ = Teuchos::null;
853 if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
861 solutionUpdated_ =
false;
869 if (LP_!=Teuchos::null) {
874 #ifdef BELOS_TEUCHOS_TIME_MONITOR 875 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
891 template <
class ScalarType,
class MV,
class OP>
898 return Teuchos::null;
902 template <
class ScalarType,
class MV,
class OP>
909 return Teuchos::null;
913 template <
class ScalarType,
class MV,
class OP>
919 const bool leftPrec = LP_ != null;
920 const bool rightPrec = RP_ != null;
931 if (!leftPrec && !rightPrec){
932 #ifdef BELOS_TEUCHOS_TIME_MONITOR 933 Teuchos::TimeMonitor OpTimer(*timerOp_);
940 else if( leftPrec && rightPrec )
943 #ifdef BELOS_TEUCHOS_TIME_MONITOR 944 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
949 #ifdef BELOS_TEUCHOS_TIME_MONITOR 950 Teuchos::TimeMonitor OpTimer(*timerOp_);
955 #ifdef BELOS_TEUCHOS_TIME_MONITOR 956 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
967 #ifdef BELOS_TEUCHOS_TIME_MONITOR 968 Teuchos::TimeMonitor OpTimer(*timerOp_);
973 #ifdef BELOS_TEUCHOS_TIME_MONITOR 974 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
985 #ifdef BELOS_TEUCHOS_TIME_MONITOR 986 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
991 #ifdef BELOS_TEUCHOS_TIME_MONITOR 992 Teuchos::TimeMonitor OpTimer(*timerOp_);
999 template <
class ScalarType,
class MV,
class OP>
1002 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1003 Teuchos::TimeMonitor OpTimer(*timerOp_);
1008 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1009 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1013 template <
class ScalarType,
class MV,
class OP>
1015 if (LP_!=Teuchos::null) {
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1017 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1022 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1023 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1027 template <
class ScalarType,
class MV,
class OP>
1029 if (RP_!=Teuchos::null) {
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1031 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1036 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1037 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1041 template <
class ScalarType,
class MV,
class OP>
1047 if (LP_!=Teuchos::null)
1051 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1052 Teuchos::TimeMonitor OpTimer(*timerOp_);
1058 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1059 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1067 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1068 Teuchos::TimeMonitor OpTimer(*timerOp_);
1077 Teuchos::RCP<const MV> localB, localX;
1079 localB = Teuchos::rcp( B,
false );
1084 localX = Teuchos::rcp( X,
false );
1088 if (LP_!=Teuchos::null)
1092 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1093 Teuchos::TimeMonitor OpTimer(*timerOp_);
1099 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1100 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1109 Teuchos::TimeMonitor OpTimer(*timerOp_);
1120 template <
class ScalarType,
class MV,
class OP>
1127 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1128 Teuchos::TimeMonitor OpTimer(*timerOp_);
1136 Teuchos::RCP<const MV> localB, localX;
1138 localB = Teuchos::rcp( B,
false );
1143 localX = Teuchos::rcp( X,
false );
1148 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1149 Teuchos::TimeMonitor OpTimer(*timerOp_);
Teuchos::RCP< MV > getLHS() const
A pointer to the left-hand side X.
Exception thrown to signal error with the Belos::LinearProblem object.
void apply(const MV &x, MV &y) const
Apply the composite operator of this linear problem to x, returning y.
void setHermitian(bool isSym=true)
Tell the linear problem that the (preconditioned) operator is Hermitian.
bool isProblemSet() const
Whether the problem has been set.
Teuchos::RCP< const MV > getCurrRHSVec()
Get a pointer to the current right-hand side of the linear system.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Teuchos::RCP< MV > getCurrLHSVec()
Get a pointer to the current left-hand side (solution) of the linear system.
Teuchos::RCP< MV > updateSolution(const Teuchos::RCP< MV > &update=Teuchos::null, bool updateLP=false, ScalarType scale=Teuchos::ScalarTraits< ScalarType >::one())
Compute the new solution to the linear system using the given update vector.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
virtual ~LinearProblem(void)
Destructor (declared virtual for memory safety of derived classes).
bool setProblem(const Teuchos::RCP< MV > &newX=Teuchos::null, const Teuchos::RCP< const MV > &newB=Teuchos::null)
Set up the linear problem manager.
LinearProblem(void)
Default constructor.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
void setOperator(const Teuchos::RCP< const OP > &A)
Set the operator A of the linear problem .
Teuchos::RCP< MV > updateSolution(const Teuchos::RCP< MV > &update=Teuchos::null, ScalarType scale=Teuchos::ScalarTraits< ScalarType >::one()) const
Compute the new solution to the linear system using the given update vector.
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
Class which defines basic traits for the operator type.
void setLSIndex(const std::vector< int > &index)
Tell the linear problem which linear system(s) need to be solved next.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
void computeCurrPrecResVec(MV *R, const MV *X=0, const MV *B=0) const
Compute a residual R for this operator given a solution X, and right-hand side B. ...
void setRHS(const Teuchos::RCP< const MV > &B)
Set right-hand-side B of linear problem .
bool isRightPrec() const
Whether the linear system is being preconditioned on the right.
Traits class which defines basic operations on multivectors.
void applyOp(const MV &x, MV &y) const
Apply ONLY the operator to x, returning y.
Teuchos::RCP< const OP > getOperator() const
A pointer to the (unpreconditioned) operator A.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
bool isSolutionUpdated() const
Has the current approximate solution been updated?
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
void setLeftPrec(const Teuchos::RCP< const OP > &LP)
Set left preconditioner (LP) of linear problem .
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
A linear system to solve, and its associated information.
int getLSNumber() const
The number of linear systems that have been set.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
The timers for this object.
bool isLeftPrec() const
Whether the linear system is being preconditioned on the left.
void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
void setLabel(const std::string &label)
Set the label prefix used by the timers in this object.
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
void applyLeftPrec(const MV &x, MV &y) const
Apply ONLY the left preconditioner to x, returning y.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
const std::vector< int > getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
void computeCurrResVec(MV *R, const MV *X=0, const MV *B=0) const
Compute a residual R for this operator given a solution X, and right-hand side B. ...
LinearProblemError(const std::string &what_arg)
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
void setCurrLS()
Tell the linear problem that the solver is finished with the current linear system.
Class which defines basic traits for the operator type.
bool isHermitian() const
Whether the (preconditioned) operator is Hermitian.
Parent class to all Belos exceptions.
void setRightPrec(const Teuchos::RCP< const OP > &RP)
Set right preconditioner (RP) of linear problem .
void setLHS(const Teuchos::RCP< MV > &X)
Set left-hand-side X of linear problem .