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>
103 const Teuchos::RCP<MV> &X,
104 const Teuchos::RCP<const MV> &B);
133 void setLHS (
const Teuchos::RCP<MV> &X) {
142 void setRHS (
const Teuchos::RCP<const MV> &B) {
196 virtual void setLSIndex (
const std::vector<int>& index);
216 void setLabel (
const std::string& label);
257 bool updateLP =
false,
258 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
280 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() )
const 315 setProblem (
const Teuchos::RCP<MV> &newX = Teuchos::null,
316 const Teuchos::RCP<const MV> &newB = Teuchos::null);
330 Teuchos::RCP<const MV>
getRHS()
const {
return(
B_); }
416 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
464 virtual void apply(
const MV& x, MV& y )
const;
473 virtual void applyOp(
const MV& x, MV& y )
const;
496 virtual void computeCurrResVec( MV* R ,
const MV* X = 0,
const MV* B = 0 )
const;
510 Teuchos::RCP<const OP>
A_;
519 Teuchos::RCP<const MV>
B_;
537 Teuchos::RCP<const OP>
LP_;
540 Teuchos::RCP<const OP>
RP_;
585 template <
class ScalarType,
class MV,
class OP>
593 solutionUpdated_(false),
598 template <
class ScalarType,
class MV,
class OP>
600 const Teuchos::RCP<MV> &X,
601 const Teuchos::RCP<const MV> &B
612 solutionUpdated_(false),
617 template <
class ScalarType,
class MV,
class OP>
631 timerPrec_(
Problem.timerPrec_),
632 blocksize_(
Problem.blocksize_),
633 num2Solve_(
Problem.num2Solve_),
637 isHermitian_(
Problem.isHermitian_),
638 solutionUpdated_(
Problem.solutionUpdated_),
643 template <
class ScalarType,
class MV,
class OP>
651 curB_ = Teuchos::null;
652 curX_ = Teuchos::null;
655 int validIdx = 0, ivalidIdx = 0;
656 blocksize_ = rhsIndex_.size();
657 std::vector<int> vldIndex( blocksize_ );
658 std::vector<int> newIndex( blocksize_ );
659 std::vector<int> iIndex( blocksize_ );
660 for (
int i=0; i<blocksize_; ++i) {
661 if (rhsIndex_[i] > -1) {
662 vldIndex[validIdx] = rhsIndex_[i];
663 newIndex[validIdx] = i;
667 iIndex[ivalidIdx] = i;
671 vldIndex.resize(validIdx);
672 newIndex.resize(validIdx);
673 iIndex.resize(ivalidIdx);
674 num2Solve_ = validIdx;
677 if (num2Solve_ != blocksize_) {
678 newIndex.resize(num2Solve_);
679 vldIndex.resize(num2Solve_);
683 curX_ = MVT::Clone( *X_, blocksize_ );
685 Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
686 MVT::MvRandom(*tmpCurB);
689 Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
690 MVT::SetBlock( *tptr, newIndex, *tmpCurB );
694 tptr = MVT::CloneView( *X_, vldIndex );
695 MVT::SetBlock( *tptr, newIndex, *curX_ );
697 solutionUpdated_ =
false;
700 curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
701 curB_ = MVT::CloneView( *B_, rhsIndex_ );
710 template <
class ScalarType,
class MV,
class OP>
717 if (num2Solve_ < blocksize_) {
722 std::vector<int> newIndex( num2Solve_ );
723 std::vector<int> vldIndex( num2Solve_ );
724 for (
int i=0; i<blocksize_; ++i) {
725 if ( rhsIndex_[i] > -1 ) {
726 vldIndex[validIdx] = rhsIndex_[i];
727 newIndex[validIdx] = i;
731 Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
732 MVT::SetBlock( *tptr, vldIndex, *X_ );
738 curX_ = Teuchos::null;
739 curB_ = Teuchos::null;
744 template <
class ScalarType,
class MV,
class OP>
755 if (update.is_null())
768 MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
773 RCP<MV> rightPrecUpdate =
774 MVT::Clone (*update, MVT::GetNumberVecs (*update));
775 applyRightPrec( *update, *rightPrecUpdate );
777 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
779 solutionUpdated_ =
true;
786 newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
790 MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
795 RCP<MV> rightPrecUpdate =
796 MVT::Clone (*update, MVT::GetNumberVecs (*update));
797 applyRightPrec( *update, *rightPrecUpdate );
799 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
806 template <
class ScalarType,
class MV,
class OP>
809 if (label != label_) {
812 if (timerOp_ != Teuchos::null) {
813 std::string opLabel = label_ +
": Operation Op*x";
814 #ifdef BELOS_TEUCHOS_TIME_MONITOR 815 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
818 if (timerPrec_ != Teuchos::null) {
819 std::string precLabel = label_ +
": Operation Prec*x";
820 #ifdef BELOS_TEUCHOS_TIME_MONITOR 821 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
827 template <
class ScalarType,
class MV,
class OP>
831 const Teuchos::RCP<const MV> &newB)
834 if (timerOp_ == Teuchos::null) {
835 std::string opLabel = label_ +
": Operation Op*x";
836 #ifdef BELOS_TEUCHOS_TIME_MONITOR 837 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
840 if (timerPrec_ == Teuchos::null) {
841 std::string precLabel = label_ +
": Operation Prec*x";
842 #ifdef BELOS_TEUCHOS_TIME_MONITOR 843 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
848 if (newX != Teuchos::null)
850 if (newB != Teuchos::null)
855 curX_ = Teuchos::null;
856 curB_ = Teuchos::null;
860 if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
868 solutionUpdated_ =
false;
871 if(Teuchos::is_null(R0_user_)) {
872 if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
873 R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
875 computeCurrResVec( &*R0_, &*X_, &*B_ );
877 if (LP_!=Teuchos::null) {
878 if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) {
879 PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
881 applyLeftPrec( *R0_, *PR0_ );
889 if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
890 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
891 computeCurrResVec( &*helper, &*X_, &*B_ );
892 R0_user_ = Teuchos::null;
896 if (LP_!=Teuchos::null) {
899 if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (PR0_user_==R0_user_)
900 || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) {
901 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
904 applyLeftPrec( *getInitResVec(), *helper );
905 PR0_user_ = Teuchos::null;
912 if (R0_user_!=Teuchos::null)
914 PR0_user_ = R0_user_;
918 PR0_user_ = Teuchos::null;
931 template <
class ScalarType,
class MV,
class OP>
934 if(Teuchos::nonnull(R0_user_)) {
940 template <
class ScalarType,
class MV,
class OP>
943 if(Teuchos::nonnull(PR0_user_)) {
949 template <
class ScalarType,
class MV,
class OP>
956 return Teuchos::null;
960 template <
class ScalarType,
class MV,
class OP>
967 return Teuchos::null;
971 template <
class ScalarType,
class MV,
class OP>
977 const bool leftPrec = LP_ != null;
978 const bool rightPrec = RP_ != null;
984 RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
989 if (!leftPrec && !rightPrec) {
995 else if( leftPrec && rightPrec )
997 applyRightPrec( x, y );
998 applyOp( y, *ytemp );
999 applyLeftPrec( *ytemp, y );
1006 applyOp( x, *ytemp );
1007 applyLeftPrec( *ytemp, y );
1014 applyRightPrec( x, *ytemp );
1015 applyOp( *ytemp, y );
1019 template <
class ScalarType,
class MV,
class OP>
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1023 Teuchos::TimeMonitor OpTimer(*timerOp_);
1025 OPT::Apply( *A_,x, y);
1028 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1029 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1033 template <
class ScalarType,
class MV,
class OP>
1035 if (LP_!=Teuchos::null) {
1036 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1037 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1039 OPT::Apply( *LP_,x, y);
1042 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1043 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1047 template <
class ScalarType,
class MV,
class OP>
1049 if (RP_!=Teuchos::null) {
1050 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1051 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1053 OPT::Apply( *RP_,x, y);
1056 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1057 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1061 template <
class ScalarType,
class MV,
class OP>
1067 if (LP_!=Teuchos::null)
1069 Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1070 applyOp( *X, *R_temp );
1071 MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1072 applyLeftPrec( *R_temp, *R );
1077 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1082 Teuchos::RCP<const MV> localB, localX;
1084 localB = Teuchos::rcp( B,
false );
1089 localX = Teuchos::rcp( X,
false );
1093 if (LP_!=Teuchos::null)
1095 Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1096 applyOp( *localX, *R_temp );
1097 MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1098 applyLeftPrec( *R_temp, *R );
1102 applyOp( *localX, *R );
1103 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1110 template <
class ScalarType,
class MV,
class OP>
1117 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1121 Teuchos::RCP<const MV> localB, localX;
1123 localB = Teuchos::rcp( B,
false );
1128 localX = Teuchos::rcp( X,
false );
1132 applyOp( *localX, *R );
1133 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
bool isHermitian_
Whether the operator A is symmetric (in real arithmetic, or Hermitian in complex arithmetic).
bool isSet_
Has the linear problem to solve been set?
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.
Teuchos::RCP< const MV > B_
Right-hand side of linear system.
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.
virtual 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.
Teuchos::RCP< const MV > R0_user_
User-defined initial residual of the linear system.
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).
virtual 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.
Teuchos::RCP< MV > curX_
Current solution vector of the linear system.
bool isHermitian() const
Whether the (preconditioned) operator is Hermitian.
Teuchos::RCP< const OP > LP_
Left preconditioning operator of linear system.
virtual 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.
virtual void setLSIndex(const std::vector< int > &index)
Tell the linear problem which linear system(s) need to be solved next.
void setRHS(const Teuchos::RCP< const MV > &B)
Set right-hand-side B of linear problem .
Traits class which defines basic operations on multivectors.
int blocksize_
Current block size of linear system.
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.
virtual void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
virtual void applyOp(const MV &x, MV &y) const
Apply ONLY the operator to x, returning y.
std::string label_
Linear problem label that prefixes the timer labels.
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
int lsNum_
Number of linear systems that have been loaded in this linear problem object.
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 MV > PR0_user_
User-defined preconditioned initial residual of the linear system.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
Teuchos::RCP< MV > R0_
Initial residual of the linear system.
Teuchos::RCP< const OP > A_
Operator of linear system.
Teuchos::RCP< MV > X_
Solution vector of linear system.
bool isProblemSet() const
Whether the problem has been set.
const std::vector< int > & getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
LinearProblemError(const std::string &what_arg)
bool solutionUpdated_
Has the current approximate solution been updated?
bool isSolutionUpdated() const
Has the current approximate solution been updated?
virtual 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.
virtual 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. ...
Teuchos::RCP< MV > PR0_
Preconditioned initial residual of the linear system.
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 .
Teuchos::RCP< Teuchos::Time > timerOp_
Timers.
virtual void apply(const MV &x, MV &y) const
Apply the composite operator of this linear problem to x, returning y.
Teuchos::RCP< const OP > RP_
Right preconditioning operator of linear system.
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
std::vector< int > rhsIndex_
Indices of current linear systems being solver for.
Teuchos::RCP< const MV > curB_
Current right-hand side of the linear system.
virtual 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.
int num2Solve_
Number of linear systems that are currently being solver for ( <= blocksize_ )
void setInitPrecResVec(const Teuchos::RCP< const MV > &PR0)
Set the user-defined preconditioned residual of linear problem .
Teuchos::RCP< Teuchos::Time > timerPrec_
virtual 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. ...