Belos Package Browser (Single Doxygen Collection)  Development
BelosLinearProblem.hpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // Belos: Block Linear Solvers Package
6 // Copyright 2004 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
43 #ifndef BELOS_LINEAR_PROBLEM_HPP
44 #define BELOS_LINEAR_PROBLEM_HPP
45 
49 #include "BelosMultiVecTraits.hpp"
50 #include "BelosOperatorTraits.hpp"
51 #include "Teuchos_ParameterList.hpp"
52 #include "Teuchos_TimeMonitor.hpp"
53 
54 namespace Belos {
55 
57 
58 
61  class LinearProblemError : public BelosError {
62  public:
63  LinearProblemError (const std::string& what_arg) :
64  BelosError(what_arg) {}
65  };
66 
68 
81  template <class ScalarType, class MV, class OP>
82  class LinearProblem {
83  public:
84 
86 
87 
94  LinearProblem (void);
95 
103  LinearProblem (const Teuchos::RCP<const OP> &A,
104  const Teuchos::RCP<MV> &X,
105  const Teuchos::RCP<const MV> &B);
106 
110  LinearProblem (const LinearProblem<ScalarType,MV,OP>& Problem);
111 
113  virtual ~LinearProblem (void);
114 
116 
118 
119 
123  void setOperator (const Teuchos::RCP<const OP> &A) {
124  A_ = A;
125  isSet_=false;
126  }
127 
134  void setLHS (const Teuchos::RCP<MV> &X) {
135  X_ = X;
136  isSet_=false;
137  }
138 
143  void setRHS (const Teuchos::RCP<const MV> &B) {
144  B_ = B;
145  isSet_=false;
146  }
147 
151  void setInitResVec(const Teuchos::RCP<const MV> &R0) {
152  R0_user_ = R0;
153  isSet_=false;
154  }
155 
159  void setInitPrecResVec(const Teuchos::RCP<const MV> &PR0) {
160  PR0_user_ = PR0;
161  isSet_=false;
162  }
163 
167  void setLeftPrec(const Teuchos::RCP<const OP> &LP) { LP_ = LP; }
168 
172  void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }
173 
182  void setCurrLS ();
183 
193  void setLSIndex (const std::vector<int>& index);
194 
205  void setHermitian( bool isSym = true ) { isHermitian_ = isSym; }
206 
213  void setLabel (const std::string& label);
214 
251  Teuchos::RCP<MV>
252  updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
253  bool updateLP = false,
254  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
255 
273  Teuchos::RCP<MV> updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
274  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
275  { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
276 
278 
280 
281 
307  bool
308  setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null,
309  const Teuchos::RCP<const MV> &newB = Teuchos::null);
310 
312 
314 
315 
317  Teuchos::RCP<const OP> getOperator() const { return(A_); }
318 
320  Teuchos::RCP<MV> getLHS() const { return(X_); }
321 
323  Teuchos::RCP<const MV> getRHS() const { return(B_); }
324 
326  Teuchos::RCP<const MV> getInitResVec() const;
327 
332  Teuchos::RCP<const MV> getInitPrecResVec() const;
333 
348  Teuchos::RCP<MV> getCurrLHSVec();
349 
364  Teuchos::RCP<const MV> getCurrRHSVec();
365 
367  Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
368 
370  Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
371 
393  const std::vector<int> getLSIndex() const { return(rhsIndex_); }
394 
401  int getLSNumber() const { return(lsNum_); }
402 
409  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
410  return Teuchos::tuple(timerOp_,timerPrec_);
411  }
412 
413 
415 
417 
418 
426  bool isSolutionUpdated() const { return(solutionUpdated_); }
427 
429  bool isProblemSet() const { return(isSet_); }
430 
436  bool isHermitian() const { return(isHermitian_); }
437 
439  bool isLeftPrec() const { return(LP_!=Teuchos::null); }
440 
442  bool isRightPrec() const { return(RP_!=Teuchos::null); }
443 
445 
447 
448 
450 
457  void apply( const MV& x, MV& y ) const;
458 
466  void applyOp( const MV& x, MV& y ) const;
467 
474  void applyLeftPrec( const MV& x, MV& y ) const;
475 
482  void applyRightPrec( const MV& x, MV& y ) const;
483 
485 
489  void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
490 
492 
496  void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
497 
499 
500  private:
501 
503  Teuchos::RCP<const OP> A_;
504 
506  Teuchos::RCP<MV> X_;
507 
509  Teuchos::RCP<MV> curX_;
510 
512  Teuchos::RCP<const MV> B_;
513 
515  Teuchos::RCP<const MV> curB_;
516 
518  Teuchos::RCP<MV> R0_;
519 
521  Teuchos::RCP<MV> PR0_;
522 
524  Teuchos::RCP<const MV> R0_user_;
525 
527  Teuchos::RCP<const MV> PR0_user_;
528 
530  Teuchos::RCP<const OP> LP_;
531 
533  Teuchos::RCP<const OP> RP_;
534 
536  mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
537 
540 
543 
545  std::vector<int> rhsIndex_;
546 
548  int lsNum_;
549 
551 
552 
555 
558 
560  bool isSet_;
561 
565 
568 
570 
572  std::string label_;
573 
576  };
577 
578  //--------------------------------------------
579  // Constructor Implementations
580  //--------------------------------------------
581 
582  template <class ScalarType, class MV, class OP>
584  blocksize_(0),
585  num2Solve_(0),
586  rhsIndex_(0),
587  lsNum_(0),
588  Left_Scale_(false),
589  Right_Scale_(false),
590  isSet_(false),
591  isHermitian_(false),
592  solutionUpdated_(false),
593  label_("Belos")
594  {
595  }
596 
597  template <class ScalarType, class MV, class OP>
598  LinearProblem<ScalarType,MV,OP>::LinearProblem(const Teuchos::RCP<const OP> &A,
599  const Teuchos::RCP<MV> &X,
600  const Teuchos::RCP<const MV> &B
601  ) :
602  A_(A),
603  X_(X),
604  B_(B),
605  blocksize_(0),
606  num2Solve_(0),
607  rhsIndex_(0),
608  lsNum_(0),
609  Left_Scale_(false),
610  Right_Scale_(false),
611  isSet_(false),
612  isHermitian_(false),
613  solutionUpdated_(false),
614  label_("Belos")
615  {
616  }
617 
618  template <class ScalarType, class MV, class OP>
620  A_(Problem.A_),
621  X_(Problem.X_),
622  curX_(Problem.curX_),
623  B_(Problem.B_),
624  curB_(Problem.curB_),
625  R0_(Problem.R0_),
626  PR0_(Problem.PR0_),
627  R0_user_(Problem.R0_user_),
628  PR0_user_(Problem.PR0_user_),
629  LP_(Problem.LP_),
630  RP_(Problem.RP_),
631  timerOp_(Problem.timerOp_),
632  timerPrec_(Problem.timerPrec_),
633  blocksize_(Problem.blocksize_),
634  num2Solve_(Problem.num2Solve_),
635  rhsIndex_(Problem.rhsIndex_),
636  lsNum_(Problem.lsNum_),
637  Left_Scale_(Problem.Left_Scale_),
638  Right_Scale_(Problem.Right_Scale_),
639  isSet_(Problem.isSet_),
640  isHermitian_(Problem.isHermitian_),
641  solutionUpdated_(Problem.solutionUpdated_),
642  label_(Problem.label_)
643  {
644  }
645 
646  template <class ScalarType, class MV, class OP>
648  {}
649 
650  template <class ScalarType, class MV, class OP>
651  void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
652  {
653  // Set new linear systems using the indices in index.
654  rhsIndex_ = index;
655 
656  // Compute the new block linear system.
657  // ( first clean up old linear system )
658  curB_ = Teuchos::null;
659  curX_ = Teuchos::null;
660 
661  // Create indices for the new linear system.
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;
671  validIdx++;
672  }
673  else {
674  iIndex[ivalidIdx] = i;
675  ivalidIdx++;
676  }
677  }
678  vldIndex.resize(validIdx);
679  newIndex.resize(validIdx);
680  iIndex.resize(ivalidIdx);
681  num2Solve_ = validIdx;
682 
683  // Create the new linear system using index
684  if (num2Solve_ != blocksize_) {
685  newIndex.resize(num2Solve_);
686  vldIndex.resize(num2Solve_);
687  //
688  // First create multivectors of blocksize.
689  // Fill the RHS with random vectors LHS with zero vectors.
690  curX_ = MVT::Clone( *X_, blocksize_ );
691  MVT::MvInit(*curX_);
692  Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
693  MVT::MvRandom(*tmpCurB);
694  //
695  // Now put in the part of B into curB
696  Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
697  MVT::SetBlock( *tptr, newIndex, *tmpCurB );
698  curB_ = tmpCurB;
699  //
700  // Now put in the part of X into curX
701  tptr = MVT::CloneView( *X_, vldIndex );
702  MVT::SetBlock( *tptr, newIndex, *curX_ );
703  //
704  solutionUpdated_ = false;
705  }
706  else {
707  curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
708  curB_ = MVT::CloneView( *B_, rhsIndex_ );
709  }
710  //
711  // Increment the number of linear systems that have been loaded into this object.
712  //
713  lsNum_++;
714  }
715 
716 
717  template <class ScalarType, class MV, class OP>
719  {
720  //
721  // We only need to copy the solutions back if the linear systems of
722  // interest are less than the block size.
723  //
724  if (num2Solve_ < blocksize_) {
725  //
726  // Get a view of the current solutions and correction std::vector.
727  //
728  int validIdx = 0;
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;
735  validIdx++;
736  }
737  }
738  Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
739  MVT::SetBlock( *tptr, vldIndex, *X_ );
740  }
741  //
742  // Clear the current vectors of this linear system so that any future calls
743  // to get the vectors for this system return null pointers.
744  //
745  curX_ = Teuchos::null;
746  curB_ = Teuchos::null;
747  rhsIndex_.resize(0);
748  }
749 
750 
751  template <class ScalarType, class MV, class OP>
752  Teuchos::RCP<MV>
754  updateSolution (const Teuchos::RCP<MV>& update,
755  bool updateLP,
756  ScalarType scale)
757  {
758  using Teuchos::RCP;
759  using Teuchos::null;
760 
761  RCP<MV> newSoln;
762  if (update.is_null())
763  { // The caller didn't supply an update vector, so we assume
764  // that the current solution curX_ is unchanged, and return a
765  // pointer to it.
766  newSoln = curX_;
767  }
768  else // the update vector is NOT null.
769  {
770  if (updateLP) // The caller wants us to update the linear problem.
771  {
772  if (RP_.is_null())
773  { // There is no right preconditioner.
774  // curX_ := curX_ + scale * update.
775  MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
776  }
777  else
778  { // There is indeed a right preconditioner, so apply it
779  // before computing the new solution.
780  RCP<MV> rightPrecUpdate =
781  MVT::Clone (*update, MVT::GetNumberVecs (*update));
782  {
783 #ifdef BELOS_TEUCHOS_TIME_MONITOR
784  Teuchos::TimeMonitor PrecTimer (*timerPrec_);
785 #endif
786  OPT::Apply( *RP_, *update, *rightPrecUpdate );
787  }
788  // curX_ := curX_ + scale * rightPrecUpdate.
789  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
790  }
791  solutionUpdated_ = true;
792  newSoln = curX_;
793  }
794  else
795  { // Rather than updating our stored current solution curX_,
796  // we make a copy and compute the new solution in the
797  // copy, without modifying curX_.
798  newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
799  if (RP_.is_null())
800  { // There is no right preconditioner.
801  // newSoln := curX_ + scale * update.
802  MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
803  }
804  else
805  { // There is indeed a right preconditioner, so apply it
806  // before computing the new solution.
807  RCP<MV> rightPrecUpdate =
808  MVT::Clone (*update, MVT::GetNumberVecs (*update));
809  {
810 #ifdef BELOS_TEUCHOS_TIME_MONITOR
811  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
812 #endif
813  OPT::Apply( *RP_, *update, *rightPrecUpdate );
814  }
815  // newSoln := curX_ + scale * rightPrecUpdate.
816  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
817  }
818  }
819  }
820  return newSoln;
821  }
822 
823  template <class ScalarType, class MV, class OP>
824  void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
825  {
826  if (label != label_) {
827  label_ = label;
828  // Create new timers if they have already been created.
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 );
833 #endif
834  }
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 );
839 #endif
840  }
841  }
842  }
843 
844  template <class ScalarType, class MV, class OP>
845  bool
847  setProblem (const Teuchos::RCP<MV> &newX,
848  const Teuchos::RCP<const MV> &newB)
849  {
850  // Create timers if the haven't been created yet.
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 );
855 #endif
856  }
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 );
861 #endif
862  }
863 
864  // Set the linear system using the arguments newX and newB
865  if (newX != Teuchos::null)
866  X_ = newX;
867  if (newB != Teuchos::null)
868  B_ = newB;
869 
870  // Invalidate the current linear system indices and multivectors.
871  rhsIndex_.resize(0);
872  curX_ = Teuchos::null;
873  curB_ = Teuchos::null;
874 
875  // If we didn't set a matrix A, a left-hand side X, or a
876  // right-hand side B, then we didn't set the problem.
877  if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
878  isSet_ = false;
879  return isSet_;
880  }
881 
882  // Reset whether the solution has been updated. (We're just
883  // setting the problem now, so of course we haven't updated the
884  // solution yet.)
885  solutionUpdated_ = false;
886 
887  // Compute the initial residuals.
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_ ) );
891  }
892  computeCurrResVec( &*R0_, &*X_, &*B_ );
893 
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_ ) );
897  }
898  {
899 #ifdef BELOS_TEUCHOS_TIME_MONITOR
900  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
901 #endif
902  OPT::Apply( *LP_, *R0_, *PR0_ );
903  }
904  }
905  else {
906  PR0_ = R0_;
907  }
908  }
909  else { // User specified the residuals
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_ );
913  R0_user_ = helper;
914  }
915 
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_ ) );
919  {
920 #ifdef BELOS_TEUCHOS_TIME_MONITOR
921  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
922 #endif
923  OPT::Apply( *LP_, *R0_user_, *helper );
924  }
925  PR0_user_ = helper;
926  }
927  }
928  else {
929  PR0_user_ = R0_user_;
930  }
931  }
932 
933  // The problem has been set and is ready for use.
934  isSet_ = true;
935 
936  // Return isSet.
937  return isSet_;
938  }
939 
940  template <class ScalarType, class MV, class OP>
941  Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getInitResVec() const
942  {
943  if(Teuchos::nonnull(R0_user_)) {
944  return R0_user_;
945  }
946  return(R0_);
947  }
948 
949  template <class ScalarType, class MV, class OP>
951  {
952  if(Teuchos::nonnull(PR0_user_)) {
953  return PR0_user_;
954  }
955  return(PR0_);
956  }
957 
958  template <class ScalarType, class MV, class OP>
960  {
961  if (isSet_) {
962  return curX_;
963  }
964  else {
965  return Teuchos::null;
966  }
967  }
968 
969  template <class ScalarType, class MV, class OP>
971  {
972  if (isSet_) {
973  return curB_;
974  }
975  else {
976  return Teuchos::null;
977  }
978  }
979 
980  template <class ScalarType, class MV, class OP>
981  void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
982  {
983  using Teuchos::null;
984  using Teuchos::RCP;
985 
986  const bool leftPrec = LP_ != null;
987  const bool rightPrec = RP_ != null;
988 
989  // We only need a temporary vector for intermediate results if
990  // there is a left or right preconditioner. We really should just
991  // keep this temporary vector around instead of allocating it each
992  // time.
993  RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
994 
995  //
996  // No preconditioning.
997  //
998  if (!leftPrec && !rightPrec){
999 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1000  Teuchos::TimeMonitor OpTimer(*timerOp_);
1001 #endif
1002  OPT::Apply( *A_, x, y );
1003  }
1004  //
1005  // Preconditioning is being done on both sides
1006  //
1007  else if( leftPrec && rightPrec )
1008  {
1009  {
1010 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1011  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1012 #endif
1013  OPT::Apply( *RP_, x, y );
1014  }
1015  {
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1017  Teuchos::TimeMonitor OpTimer(*timerOp_);
1018 #endif
1019  OPT::Apply( *A_, y, *ytemp );
1020  }
1021  {
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1023  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1024 #endif
1025  OPT::Apply( *LP_, *ytemp, y );
1026  }
1027  }
1028  //
1029  // Preconditioning is only being done on the left side
1030  //
1031  else if( leftPrec )
1032  {
1033  {
1034 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1035  Teuchos::TimeMonitor OpTimer(*timerOp_);
1036 #endif
1037  OPT::Apply( *A_, x, *ytemp );
1038  }
1039  {
1040 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1041  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1042 #endif
1043  OPT::Apply( *LP_, *ytemp, y );
1044  }
1045  }
1046  //
1047  // Preconditioning is only being done on the right side
1048  //
1049  else
1050  {
1051  {
1052 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1053  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1054 #endif
1055  OPT::Apply( *RP_, x, *ytemp );
1056  }
1057  {
1058 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1059  Teuchos::TimeMonitor OpTimer(*timerOp_);
1060 #endif
1061  OPT::Apply( *A_, *ytemp, y );
1062  }
1063  }
1064  }
1065 
1066  template <class ScalarType, class MV, class OP>
1067  void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
1068  if (A_.get()) {
1069 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1070  Teuchos::TimeMonitor OpTimer(*timerOp_);
1071 #endif
1072  OPT::Apply( *A_,x, y);
1073  }
1074  else {
1075  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1076  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1077  }
1078  }
1079 
1080  template <class ScalarType, class MV, class OP>
1081  void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
1082  if (LP_!=Teuchos::null) {
1083 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1084  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1085 #endif
1086  return ( OPT::Apply( *LP_,x, y) );
1087  }
1088  else {
1089  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1090  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1091  }
1092  }
1093 
1094  template <class ScalarType, class MV, class OP>
1095  void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
1096  if (RP_!=Teuchos::null) {
1097 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1098  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1099 #endif
1100  return ( OPT::Apply( *RP_,x, y) );
1101  }
1102  else {
1103  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1104  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1105  }
1106  }
1107 
1108  template <class ScalarType, class MV, class OP>
1109  void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
1110 
1111  if (R) {
1112  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1113  {
1114  if (LP_!=Teuchos::null)
1115  {
1116  Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1117  {
1118 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1119  Teuchos::TimeMonitor OpTimer(*timerOp_);
1120 #endif
1121  OPT::Apply( *A_, *X, *R_temp );
1122  }
1123  MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1124  {
1125 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1126  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1127 #endif
1128  OPT::Apply( *LP_, *R_temp, *R );
1129  }
1130  }
1131  else
1132  {
1133  {
1134 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1135  Teuchos::TimeMonitor OpTimer(*timerOp_);
1136 #endif
1137  OPT::Apply( *A_, *X, *R );
1138  }
1139  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1140  }
1141  }
1142  else {
1143  // The solution and right-hand side may not be specified, check and use which ones exist.
1144  Teuchos::RCP<const MV> localB, localX;
1145  if (B)
1146  localB = Teuchos::rcp( B, false );
1147  else
1148  localB = curB_;
1149 
1150  if (X)
1151  localX = Teuchos::rcp( X, false );
1152  else
1153  localX = curX_;
1154 
1155  if (LP_!=Teuchos::null)
1156  {
1157  Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1158  {
1159 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1160  Teuchos::TimeMonitor OpTimer(*timerOp_);
1161 #endif
1162  OPT::Apply( *A_, *localX, *R_temp );
1163  }
1164  MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1165  {
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1167  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1168 #endif
1169  OPT::Apply( *LP_, *R_temp, *R );
1170  }
1171  }
1172  else
1173  {
1174  {
1175 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1176  Teuchos::TimeMonitor OpTimer(*timerOp_);
1177 #endif
1178  OPT::Apply( *A_, *localX, *R );
1179  }
1180  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1181  }
1182  }
1183  }
1184  }
1185 
1186 
1187  template <class ScalarType, class MV, class OP>
1188  void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
1189 
1190  if (R) {
1191  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1192  {
1193  {
1194 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1195  Teuchos::TimeMonitor OpTimer(*timerOp_);
1196 #endif
1197  OPT::Apply( *A_, *X, *R );
1198  }
1199  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1200  }
1201  else {
1202  // The solution and right-hand side may not be specified, check and use which ones exist.
1203  Teuchos::RCP<const MV> localB, localX;
1204  if (B)
1205  localB = Teuchos::rcp( B, false );
1206  else
1207  localB = curB_;
1208 
1209  if (X)
1210  localX = Teuchos::rcp( X, false );
1211  else
1212  localX = curX_;
1213 
1214  {
1215 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1216  Teuchos::TimeMonitor OpTimer(*timerOp_);
1217 #endif
1218  OPT::Apply( *A_, *localX, *R );
1219  }
1220  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1221  }
1222  }
1223  }
1224 
1225 } // end Belos namespace
1226 
1227 #endif /* BELOS_LINEAR_PROBLEM_HPP */
1228 
1229 
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.
Definition: BelosTypes.hpp:60
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. ...