43 #ifndef BELOS_LINEAR_PROBLEM_HPP 44 #define BELOS_LINEAR_PROBLEM_HPP 81 template <
class ScalarType,
class MV,
class OP>
196 void setLSIndex (
const std::vector<int>& index);
216 void setLabel (
const std::string& label);
256 bool updateLP =
false,
396 const std::vector<int>
getLSIndex()
const {
return(rhsIndex_); }
413 return Teuchos::tuple(timerOp_,timerPrec_);
460 void apply(
const MV& x, MV& y )
const;
469 void applyOp(
const MV& x, MV& y )
const;
548 std::vector<int> rhsIndex_;
564 bool solutionUpdated_;
579 template <
class ScalarType,
class MV,
class OP>
587 solutionUpdated_(false),
592 template <
class ScalarType,
class MV,
class OP>
606 solutionUpdated_(false),
611 template <
class ScalarType,
class MV,
class OP>
625 timerPrec_(
Problem.timerPrec_),
626 blocksize_(
Problem.blocksize_),
627 num2Solve_(
Problem.num2Solve_),
631 isHermitian_(
Problem.isHermitian_),
632 solutionUpdated_(
Problem.solutionUpdated_),
637 template <
class ScalarType,
class MV,
class OP>
641 template <
class ScalarType,
class MV,
class OP>
649 curB_ = Teuchos::null;
650 curX_ = Teuchos::null;
653 int validIdx = 0, ivalidIdx = 0;
654 blocksize_ = rhsIndex_.size();
655 std::vector<int> vldIndex( blocksize_ );
656 std::vector<int> newIndex( blocksize_ );
657 std::vector<int> iIndex( blocksize_ );
658 for (
int i=0; i<blocksize_; ++i) {
659 if (rhsIndex_[i] > -1) {
660 vldIndex[validIdx] = rhsIndex_[i];
661 newIndex[validIdx] = i;
665 iIndex[ivalidIdx] = i;
669 vldIndex.resize(validIdx);
670 newIndex.resize(validIdx);
671 iIndex.resize(ivalidIdx);
672 num2Solve_ = validIdx;
675 if (num2Solve_ != blocksize_) {
676 newIndex.resize(num2Solve_);
677 vldIndex.resize(num2Solve_);
681 curX_ = MVT::Clone( *X_, blocksize_ );
684 MVT::MvRandom(*tmpCurB);
688 MVT::SetBlock( *tptr, newIndex, *tmpCurB );
692 tptr = MVT::CloneView( *X_, vldIndex );
693 MVT::SetBlock( *tptr, newIndex, *curX_ );
695 solutionUpdated_ =
false;
698 curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
699 curB_ = MVT::CloneView( *B_, rhsIndex_ );
708 template <
class ScalarType,
class MV,
class OP>
715 if (num2Solve_ < blocksize_) {
720 std::vector<int> newIndex( num2Solve_ );
721 std::vector<int> vldIndex( num2Solve_ );
722 for (
int i=0; i<blocksize_; ++i) {
723 if ( rhsIndex_[i] > -1 ) {
724 vldIndex[validIdx] = rhsIndex_[i];
725 newIndex[validIdx] = i;
730 MVT::SetBlock( *tptr, vldIndex, *X_ );
736 curX_ = Teuchos::null;
737 curB_ = Teuchos::null;
742 template <
class ScalarType,
class MV,
class OP>
766 MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
771 RCP<MV> rightPrecUpdate =
772 MVT::Clone (*update, MVT::GetNumberVecs (*update));
774 #ifdef BELOS_TEUCHOS_TIME_MONITOR 777 OPT::Apply( *RP_, *update, *rightPrecUpdate );
780 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
782 solutionUpdated_ =
true;
789 newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
793 MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
798 RCP<MV> rightPrecUpdate =
799 MVT::Clone (*update, MVT::GetNumberVecs (*update));
801 #ifdef BELOS_TEUCHOS_TIME_MONITOR 804 OPT::Apply( *RP_, *update, *rightPrecUpdate );
807 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
814 template <
class ScalarType,
class MV,
class OP>
817 if (label != label_) {
820 if (timerOp_ != Teuchos::null) {
821 std::string opLabel = label_ +
": Operation Op*x";
822 #ifdef BELOS_TEUCHOS_TIME_MONITOR 826 if (timerPrec_ != Teuchos::null) {
827 std::string precLabel = label_ +
": Operation Prec*x";
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 835 template <
class ScalarType,
class MV,
class OP>
842 if (timerOp_ == Teuchos::null) {
843 std::string opLabel = label_ +
": Operation Op*x";
844 #ifdef BELOS_TEUCHOS_TIME_MONITOR 848 if (timerPrec_ == Teuchos::null) {
849 std::string precLabel = label_ +
": Operation Prec*x";
850 #ifdef BELOS_TEUCHOS_TIME_MONITOR 856 if (newX != Teuchos::null)
858 if (newB != Teuchos::null)
863 curX_ = Teuchos::null;
864 curB_ = Teuchos::null;
868 if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
876 solutionUpdated_ =
false;
880 if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
881 R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
883 computeCurrResVec( &*R0_, &*X_, &*B_ );
885 if (LP_!=Teuchos::null) {
886 if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) {
887 PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
890 #ifdef BELOS_TEUCHOS_TIME_MONITOR 893 OPT::Apply( *LP_, *R0_, *PR0_ );
902 if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
904 computeCurrResVec( &*helper, &*X_, &*B_ );
905 R0_user_ = Teuchos::null;
909 if (LP_!=Teuchos::null) {
912 if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (PR0_user_==R0_user_)
913 || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) {
917 #ifdef BELOS_TEUCHOS_TIME_MONITOR 920 OPT::Apply( *LP_, *getInitResVec(), *helper );
922 PR0_user_ = Teuchos::null;
929 if (R0_user_!=Teuchos::null)
931 PR0_user_ = R0_user_;
935 PR0_user_ = Teuchos::null;
948 template <
class ScalarType,
class MV,
class OP>
957 template <
class ScalarType,
class MV,
class OP>
966 template <
class ScalarType,
class MV,
class OP>
973 return Teuchos::null;
977 template <
class ScalarType,
class MV,
class OP>
984 return Teuchos::null;
988 template <
class ScalarType,
class MV,
class OP>
994 const bool leftPrec = LP_ != null;
995 const bool rightPrec = RP_ != null;
1001 RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
1006 if (!leftPrec && !rightPrec){
1007 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1010 OPT::Apply( *A_, x, y );
1015 else if( leftPrec && rightPrec )
1018 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1021 OPT::Apply( *RP_, x, y );
1024 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1027 OPT::Apply( *A_, y, *ytemp );
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1033 OPT::Apply( *LP_, *ytemp, y );
1042 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1045 OPT::Apply( *A_, x, *ytemp );
1048 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1051 OPT::Apply( *LP_, *ytemp, y );
1060 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1063 OPT::Apply( *RP_, x, *ytemp );
1066 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1069 OPT::Apply( *A_, *ytemp, y );
1074 template <
class ScalarType,
class MV,
class OP>
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1080 OPT::Apply( *A_,x, y);
1088 template <
class ScalarType,
class MV,
class OP>
1090 if (LP_!=Teuchos::null) {
1091 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1094 return ( OPT::Apply( *LP_,x, y) );
1102 template <
class ScalarType,
class MV,
class OP>
1104 if (RP_!=Teuchos::null) {
1105 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1108 return ( OPT::Apply( *RP_,x, y) );
1116 template <
class ScalarType,
class MV,
class OP>
1122 if (LP_!=Teuchos::null)
1126 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1129 OPT::Apply( *A_, *X, *R_temp );
1131 MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1133 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1136 OPT::Apply( *LP_, *R_temp, *R );
1142 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1145 OPT::Apply( *A_, *X, *R );
1147 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1163 if (LP_!=Teuchos::null)
1165 Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1167 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1170 OPT::Apply( *A_, *localX, *R_temp );
1172 MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1174 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1177 OPT::Apply( *LP_, *R_temp, *R );
1183 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1186 OPT::Apply( *A_, *localX, *R );
1188 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1195 template <
class ScalarType,
class MV,
class OP>
1202 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1205 OPT::Apply( *A_, *X, *R );
1207 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1223 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1226 OPT::Apply( *A_, *localX, *R );
1228 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
Exception thrown to signal error with the Belos::LinearProblem object.
void setHermitian(bool isSym=true)
Tell the linear problem that the (preconditioned) operator is Hermitian.
const std::vector< int > getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
Teuchos::RCP< const MV > getCurrRHSVec()
Get a pointer to the current right-hand side of the linear system.
Teuchos::RCP< MV > getCurrLHSVec()
Get a pointer to the current left-hand side (solution) of the linear system.
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
bool nonnull(const std::shared_ptr< T > &p)
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.
bool is_null(const std::shared_ptr< T > &p)
bool isRightPrec() const
Whether the linear system is being preconditioned on the right.
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.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
The timers for this object.
bool isHermitian() const
Whether the (preconditioned) operator is Hermitian.
void applyLeftPrec(const MV &x, MV &y) const
Apply ONLY the left preconditioner to x, returning y.
LinearProblem(void)
Default constructor.
void setOperator(const Teuchos::RCP< const OP > &A)
Set the operator A of the linear problem .
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.
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.
void setRHS(const Teuchos::RCP< const MV > &B)
Set right-hand-side B of linear problem .
Traits class which defines basic operations on multivectors.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setLeftPrec(const Teuchos::RCP< const OP > &LP)
Set left preconditioner (LP) of linear problem .
Teuchos::RCP< MV > getLHS() const
A pointer to the left-hand side X.
A linear system to solve, and its associated information.
void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
void applyOp(const MV &x, MV &y) const
Apply ONLY the operator to x, returning y.
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
void setLabel(const std::string &label)
Set the label prefix used by the timers in this object.
int getLSNumber() const
The number of linear systems that have been set.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
bool isProblemSet() const
Whether the problem has been set.
LinearProblemError(const std::string &what_arg)
bool isSolutionUpdated() const
Has the current approximate solution been updated?
void setCurrLS()
Tell the linear problem that the solver is finished with the current linear system.
bool isLeftPrec() const
Whether the linear system is being preconditioned on the left.
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. ...
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Teuchos::RCP< const OP > getOperator() const
A pointer to the (unpreconditioned) operator A.
void setRightPrec(const Teuchos::RCP< const OP > &RP)
Set right preconditioner (RP) of linear problem .
void setInitResVec(const Teuchos::RCP< const MV > &R0)
Set the user-defined residual of linear problem .
void setLHS(const Teuchos::RCP< MV > &X)
Set left-hand-side X of linear problem .
void apply(const MV &x, MV &y) const
Apply the composite operator of this linear problem to x, returning y.
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
void setInitPrecResVec(const Teuchos::RCP< const MV > &PR0)
Set the user-defined preconditioned residual of linear problem .
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. ...