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) {
193 void setLSIndex (
const std::vector<int>& index);
213 void setLabel (
const std::string& label);
253 bool updateLP =
false,
254 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
274 ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() )
const 308 setProblem (
const Teuchos::RCP<MV> &newX = Teuchos::null,
309 const Teuchos::RCP<const MV> &newB = Teuchos::null);
323 Teuchos::RCP<const MV>
getRHS()
const {
return(
B_); }
409 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
457 void apply(
const MV& x, MV& y )
const;
466 void applyOp(
const MV& x, MV& y )
const;
503 Teuchos::RCP<const OP>
A_;
512 Teuchos::RCP<const MV>
B_;
530 Teuchos::RCP<const OP>
LP_;
533 Teuchos::RCP<const OP>
RP_;
582 template <
class ScalarType,
class MV,
class OP>
592 solutionUpdated_(false),
597 template <
class ScalarType,
class MV,
class OP>
599 const Teuchos::RCP<MV> &X,
600 const Teuchos::RCP<const MV> &B
613 solutionUpdated_(false),
618 template <
class ScalarType,
class MV,
class OP>
632 timerPrec_(
Problem.timerPrec_),
633 blocksize_(
Problem.blocksize_),
634 num2Solve_(
Problem.num2Solve_),
637 Left_Scale_(
Problem.Left_Scale_),
638 Right_Scale_(
Problem.Right_Scale_),
640 isHermitian_(
Problem.isHermitian_),
641 solutionUpdated_(
Problem.solutionUpdated_),
646 template <
class ScalarType,
class MV,
class OP>
650 template <
class ScalarType,
class MV,
class OP>
658 curB_ = Teuchos::null;
659 curX_ = Teuchos::null;
662 int validIdx = 0, ivalidIdx = 0;
663 blocksize_ = rhsIndex_.size();
664 std::vector<int> vldIndex( blocksize_ );
665 std::vector<int> newIndex( blocksize_ );
666 std::vector<int> iIndex( blocksize_ );
667 for (
int i=0; i<blocksize_; ++i) {
668 if (rhsIndex_[i] > -1) {
669 vldIndex[validIdx] = rhsIndex_[i];
670 newIndex[validIdx] = i;
674 iIndex[ivalidIdx] = i;
678 vldIndex.resize(validIdx);
679 newIndex.resize(validIdx);
680 iIndex.resize(ivalidIdx);
681 num2Solve_ = validIdx;
684 if (num2Solve_ != blocksize_) {
685 newIndex.resize(num2Solve_);
686 vldIndex.resize(num2Solve_);
690 curX_ = MVT::Clone( *X_, blocksize_ );
692 Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
693 MVT::MvRandom(*tmpCurB);
696 Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
697 MVT::SetBlock( *tptr, newIndex, *tmpCurB );
701 tptr = MVT::CloneView( *X_, vldIndex );
702 MVT::SetBlock( *tptr, newIndex, *curX_ );
704 solutionUpdated_ =
false;
707 curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
708 curB_ = MVT::CloneView( *B_, rhsIndex_ );
717 template <
class ScalarType,
class MV,
class OP>
724 if (num2Solve_ < blocksize_) {
729 std::vector<int> newIndex( num2Solve_ );
730 std::vector<int> vldIndex( num2Solve_ );
731 for (
int i=0; i<blocksize_; ++i) {
732 if ( rhsIndex_[i] > -1 ) {
733 vldIndex[validIdx] = rhsIndex_[i];
734 newIndex[validIdx] = i;
738 Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
739 MVT::SetBlock( *tptr, vldIndex, *X_ );
745 curX_ = Teuchos::null;
746 curB_ = Teuchos::null;
751 template <
class ScalarType,
class MV,
class OP>
762 if (update.is_null())
775 MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
780 RCP<MV> rightPrecUpdate =
781 MVT::Clone (*update, MVT::GetNumberVecs (*update));
783 #ifdef BELOS_TEUCHOS_TIME_MONITOR 784 Teuchos::TimeMonitor PrecTimer (*timerPrec_);
786 OPT::Apply( *RP_, *update, *rightPrecUpdate );
789 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
791 solutionUpdated_ =
true;
798 newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
802 MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
807 RCP<MV> rightPrecUpdate =
808 MVT::Clone (*update, MVT::GetNumberVecs (*update));
810 #ifdef BELOS_TEUCHOS_TIME_MONITOR 811 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
813 OPT::Apply( *RP_, *update, *rightPrecUpdate );
816 MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
823 template <
class ScalarType,
class MV,
class OP>
826 if (label != label_) {
829 if (timerOp_ != Teuchos::null) {
830 std::string opLabel = label_ +
": Operation Op*x";
831 #ifdef BELOS_TEUCHOS_TIME_MONITOR 832 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
835 if (timerPrec_ != Teuchos::null) {
836 std::string precLabel = label_ +
": Operation Prec*x";
837 #ifdef BELOS_TEUCHOS_TIME_MONITOR 838 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
844 template <
class ScalarType,
class MV,
class OP>
848 const Teuchos::RCP<const MV> &newB)
851 if (timerOp_ == Teuchos::null) {
852 std::string opLabel = label_ +
": Operation Op*x";
853 #ifdef BELOS_TEUCHOS_TIME_MONITOR 854 timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
857 if (timerPrec_ == Teuchos::null) {
858 std::string precLabel = label_ +
": Operation Prec*x";
859 #ifdef BELOS_TEUCHOS_TIME_MONITOR 860 timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
865 if (newX != Teuchos::null)
867 if (newB != Teuchos::null)
872 curX_ = Teuchos::null;
873 curB_ = Teuchos::null;
877 if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
885 solutionUpdated_ =
false;
888 if(Teuchos::is_null(R0_user_)) {
889 if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
890 R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
892 computeCurrResVec( &*R0_, &*X_, &*B_ );
894 if (LP_!=Teuchos::null) {
895 if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) {
896 PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
899 #ifdef BELOS_TEUCHOS_TIME_MONITOR 900 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
902 OPT::Apply( *LP_, *R0_, *PR0_ );
910 if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
911 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
912 computeCurrResVec( &*helper, &*X_, &*B_ );
916 if (LP_!=Teuchos::null) {
917 if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) {
918 Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
920 #ifdef BELOS_TEUCHOS_TIME_MONITOR 921 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
923 OPT::Apply( *LP_, *R0_user_, *helper );
929 PR0_user_ = R0_user_;
940 template <
class ScalarType,
class MV,
class OP>
943 if(Teuchos::nonnull(R0_user_)) {
949 template <
class ScalarType,
class MV,
class OP>
952 if(Teuchos::nonnull(PR0_user_)) {
958 template <
class ScalarType,
class MV,
class OP>
965 return Teuchos::null;
969 template <
class ScalarType,
class MV,
class OP>
976 return Teuchos::null;
980 template <
class ScalarType,
class MV,
class OP>
986 const bool leftPrec = LP_ != null;
987 const bool rightPrec = RP_ != null;
993 RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
998 if (!leftPrec && !rightPrec){
999 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1000 Teuchos::TimeMonitor OpTimer(*timerOp_);
1002 OPT::Apply( *A_, x, y );
1007 else if( leftPrec && rightPrec )
1010 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1011 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1013 OPT::Apply( *RP_, x, y );
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1017 Teuchos::TimeMonitor OpTimer(*timerOp_);
1019 OPT::Apply( *A_, y, *ytemp );
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1023 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1025 OPT::Apply( *LP_, *ytemp, y );
1034 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1035 Teuchos::TimeMonitor OpTimer(*timerOp_);
1037 OPT::Apply( *A_, x, *ytemp );
1040 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1041 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1043 OPT::Apply( *LP_, *ytemp, y );
1052 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1053 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1055 OPT::Apply( *RP_, x, *ytemp );
1058 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1059 Teuchos::TimeMonitor OpTimer(*timerOp_);
1061 OPT::Apply( *A_, *ytemp, y );
1066 template <
class ScalarType,
class MV,
class OP>
1069 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1070 Teuchos::TimeMonitor OpTimer(*timerOp_);
1072 OPT::Apply( *A_,x, y);
1075 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1076 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1080 template <
class ScalarType,
class MV,
class OP>
1082 if (LP_!=Teuchos::null) {
1083 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1084 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1086 return ( OPT::Apply( *LP_,x, y) );
1089 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1090 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1094 template <
class ScalarType,
class MV,
class OP>
1096 if (RP_!=Teuchos::null) {
1097 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1098 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1100 return ( OPT::Apply( *RP_,x, y) );
1103 MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1104 Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1108 template <
class ScalarType,
class MV,
class OP>
1114 if (LP_!=Teuchos::null)
1116 Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1118 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1119 Teuchos::TimeMonitor OpTimer(*timerOp_);
1121 OPT::Apply( *A_, *X, *R_temp );
1123 MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1125 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1126 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1128 OPT::Apply( *LP_, *R_temp, *R );
1134 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1135 Teuchos::TimeMonitor OpTimer(*timerOp_);
1137 OPT::Apply( *A_, *X, *R );
1139 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1144 Teuchos::RCP<const MV> localB, localX;
1146 localB = Teuchos::rcp( B,
false );
1151 localX = Teuchos::rcp( X,
false );
1155 if (LP_!=Teuchos::null)
1157 Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1159 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1160 Teuchos::TimeMonitor OpTimer(*timerOp_);
1162 OPT::Apply( *A_, *localX, *R_temp );
1164 MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1167 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1169 OPT::Apply( *LP_, *R_temp, *R );
1175 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1176 Teuchos::TimeMonitor OpTimer(*timerOp_);
1178 OPT::Apply( *A_, *localX, *R );
1180 MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1187 template <
class ScalarType,
class MV,
class OP>
1194 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1195 Teuchos::TimeMonitor OpTimer(*timerOp_);
1197 OPT::Apply( *A_, *X, *R );
1199 MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1203 Teuchos::RCP<const MV> localB, localX;
1205 localB = Teuchos::rcp( B,
false );
1210 localX = Teuchos::rcp( X,
false );
1215 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1216 Teuchos::TimeMonitor OpTimer(*timerOp_);
1218 OPT::Apply( *A_, *localX, *R );
1220 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.
bool Right_Scale_
Is there a right scaling?
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 .
bool Left_Scale_
Is there a left scaling?
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. ...