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 void setLSIndex (
const std::vector<int>& index);
216 void setLabel (
const std::string& label);
256 bool updateLP =
false,
257 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
277 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() )
const 311 setProblem (
const Teuchos::RCP<MV> &newX = Teuchos::null,
312 const Teuchos::RCP<const MV> &newB = Teuchos::null);
326 Teuchos::RCP<const MV>
getRHS()
const {
return(
B_); }
412 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
460 void apply(
const MV& x, MV& y )
const;
469 void applyOp(
const MV& x, MV& y )
const;
506 Teuchos::RCP<const OP>
A_;
515 Teuchos::RCP<const MV>
B_;
533 Teuchos::RCP<const OP>
LP_;
536 Teuchos::RCP<const OP>
RP_;
579 template <
class ScalarType,
class MV,
class OP>
587 solutionUpdated_(false),
592 template <
class ScalarType,
class MV,
class OP>
594 const Teuchos::RCP<MV> &X,
595 const Teuchos::RCP<const MV> &B
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_ );
683 Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
684 MVT::MvRandom(*tmpCurB);
687 Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
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;
729 Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
730 MVT::SetBlock( *tptr, vldIndex, *X_ );
736 curX_ = Teuchos::null;
737 curB_ = Teuchos::null;
742 template <
class ScalarType,
class MV,
class OP>
753 if (update.is_null())
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 775 Teuchos::TimeMonitor PrecTimer (*timerPrec_);
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 802 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
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 823 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
826 if (timerPrec_ != Teuchos::null) {
827 std::string precLabel = label_ +
": Operation Prec*x";
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 829 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
835 template <
class ScalarType,
class MV,
class OP>
839 const Teuchos::RCP<const MV> &newB)
842 if (timerOp_ == Teuchos::null) {
843 std::string opLabel = label_ +
": Operation Op*x";
844 #ifdef BELOS_TEUCHOS_TIME_MONITOR 845 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
848 if (timerPrec_ == Teuchos::null) {
849 std::string precLabel = label_ +
": Operation Prec*x";
850 #ifdef BELOS_TEUCHOS_TIME_MONITOR 851 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
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;
879 if(Teuchos::is_null(R0_user_)) {
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 891 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
893 OPT::Apply( *LP_, *R0_, *PR0_ );
902 if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
903 Teuchos::RCP<MV> helper = MVT::Clone( *B_, 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_))) {
914 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
917 #ifdef BELOS_TEUCHOS_TIME_MONITOR 918 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
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>
951 if(Teuchos::nonnull(R0_user_)) {
957 template <
class ScalarType,
class MV,
class OP>
960 if(Teuchos::nonnull(PR0_user_)) {
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 1008 Teuchos::TimeMonitor OpTimer(*timerOp_);
1010 OPT::Apply( *A_, x, y );
1015 else if( leftPrec && rightPrec )
1018 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1019 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1021 OPT::Apply( *RP_, x, y );
1024 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1025 Teuchos::TimeMonitor OpTimer(*timerOp_);
1027 OPT::Apply( *A_, y, *ytemp );
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1031 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1033 OPT::Apply( *LP_, *ytemp, y );
1042 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1043 Teuchos::TimeMonitor OpTimer(*timerOp_);
1045 OPT::Apply( *A_, x, *ytemp );
1048 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1049 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1051 OPT::Apply( *LP_, *ytemp, y );
1060 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1061 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1063 OPT::Apply( *RP_, x, *ytemp );
1066 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1067 Teuchos::TimeMonitor OpTimer(*timerOp_);
1069 OPT::Apply( *A_, *ytemp, y );
1074 template <
class ScalarType,
class MV,
class OP>
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1078 Teuchos::TimeMonitor OpTimer(*timerOp_);
1080 OPT::Apply( *A_,x, y);
1083 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1084 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1088 template <
class ScalarType,
class MV,
class OP>
1090 if (LP_!=Teuchos::null) {
1091 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1092 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1094 return ( OPT::Apply( *LP_,x, y) );
1097 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1098 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1102 template <
class ScalarType,
class MV,
class OP>
1104 if (RP_!=Teuchos::null) {
1105 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1106 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1108 return ( OPT::Apply( *RP_,x, y) );
1111 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1112 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1116 template <
class ScalarType,
class MV,
class OP>
1122 if (LP_!=Teuchos::null)
1124 Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1126 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1127 Teuchos::TimeMonitor OpTimer(*timerOp_);
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 1134 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1136 OPT::Apply( *LP_, *R_temp, *R );
1142 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1143 Teuchos::TimeMonitor OpTimer(*timerOp_);
1145 OPT::Apply( *A_, *X, *R );
1147 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1152 Teuchos::RCP<const MV> localB, localX;
1154 localB = Teuchos::rcp( B,
false );
1159 localX = Teuchos::rcp( X,
false );
1163 if (LP_!=Teuchos::null)
1165 Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1167 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1168 Teuchos::TimeMonitor OpTimer(*timerOp_);
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 1175 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1177 OPT::Apply( *LP_, *R_temp, *R );
1183 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1184 Teuchos::TimeMonitor OpTimer(*timerOp_);
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 1203 Teuchos::TimeMonitor OpTimer(*timerOp_);
1205 OPT::Apply( *A_, *X, *R );
1207 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1211 Teuchos::RCP<const MV> localB, localX;
1213 localB = Teuchos::rcp( B,
false );
1218 localX = Teuchos::rcp( X,
false );
1223 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1224 Teuchos::TimeMonitor OpTimer(*timerOp_);
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.
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.
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.
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).
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.
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.
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.
void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
MultiVecTraits< ScalarType, MV > MVT
OperatorTraits< ScalarType, MV, OP > OPT
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.
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?
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. ...
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.
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.
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_
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. ...