Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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  const Teuchos::RCP<MV> &X,
104  const Teuchos::RCP<const MV> &B);
105 
109  LinearProblem (const LinearProblem<ScalarType,MV,OP>& Problem);
110 
112  virtual ~LinearProblem (void) {}
113 
115 
117 
118 
123  A_ = A;
124  isSet_=false;
125  }
126 
133  void setLHS (const Teuchos::RCP<MV> &X) {
134  X_ = X;
135  isSet_=false;
136  }
137 
142  void setRHS (const Teuchos::RCP<const MV> &B) {
143  B_ = B;
144  isSet_=false;
145  }
146 
153  R0_user_ = R0;
154  isSet_=false;
155  }
156 
163  PR0_user_ = PR0;
164  isSet_=false;
165  }
166 
170  void setLeftPrec(const Teuchos::RCP<const OP> &LP) { LP_ = LP; }
171 
175  void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }
176 
185  virtual void setCurrLS ();
186 
196  virtual void setLSIndex (const std::vector<int>& index);
197 
208  void setHermitian( bool isSym = true ) { isHermitian_ = isSym; }
209 
216  void setLabel (const std::string& label);
217 
254  virtual
257  bool updateLP = false,
258  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
259 
277  virtual
280  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
281  { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
282 
284 
286 
287 
313  virtual
314  bool
315  setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null,
316  const Teuchos::RCP<const MV> &newB = Teuchos::null);
317 
319 
321 
322 
324  Teuchos::RCP<const OP> getOperator() const { return(A_); }
325 
327  Teuchos::RCP<MV> getLHS() const { return(X_); }
328 
330  Teuchos::RCP<const MV> getRHS() const { return(B_); }
331 
334 
340 
356 
372 
374  Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
375 
378 
400  const std::vector<int>& getLSIndex() const { return(rhsIndex_); }
401 
408  int getLSNumber() const { return(lsNum_); }
409 
417  return Teuchos::tuple(timerOp_,timerPrec_);
418  }
419 
420 
422 
424 
425 
433  bool isSolutionUpdated() const { return(solutionUpdated_); }
434 
436  bool isProblemSet() const { return(isSet_); }
437 
443  bool isHermitian() const { return(isHermitian_); }
444 
446  bool isLeftPrec() const { return(LP_!=Teuchos::null); }
447 
449  bool isRightPrec() const { return(RP_!=Teuchos::null); }
450 
452 
454 
455 
457 
464  virtual void apply( const MV& x, MV& y ) const;
465 
473  virtual void applyOp( const MV& x, MV& y ) const;
474 
481  virtual void applyLeftPrec( const MV& x, MV& y ) const;
482 
489  virtual void applyRightPrec( const MV& x, MV& y ) const;
490 
492 
496  virtual void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
497 
499 
503  virtual void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
504 
506 
507  protected:
508 
511 
514 
517 
520 
523 
526 
529 
532 
535 
538 
541 
544 
547 
550 
552  std::vector<int> rhsIndex_;
553 
555  int lsNum_;
556 
558 
559 
561  bool isSet_;
562 
566 
569 
571 
573  std::string label_;
574 
575  private:
576 
579  };
580 
581  //--------------------------------------------
582  // Constructor Implementations
583  //--------------------------------------------
584 
585  template <class ScalarType, class MV, class OP>
587  blocksize_(0),
588  num2Solve_(0),
589  rhsIndex_(0),
590  lsNum_(0),
591  isSet_(false),
592  isHermitian_(false),
593  solutionUpdated_(false),
594  label_("Belos")
595  {
596  }
597 
598  template <class ScalarType, class MV, class OP>
600  const Teuchos::RCP<MV> &X,
602  ) :
603  A_(A),
604  X_(X),
605  B_(B),
606  blocksize_(0),
607  num2Solve_(0),
608  rhsIndex_(0),
609  lsNum_(0),
610  isSet_(false),
611  isHermitian_(false),
612  solutionUpdated_(false),
613  label_("Belos")
614  {
615  }
616 
617  template <class ScalarType, class MV, class OP>
619  A_(Problem.A_),
620  X_(Problem.X_),
621  curX_(Problem.curX_),
622  B_(Problem.B_),
623  curB_(Problem.curB_),
624  R0_(Problem.R0_),
625  PR0_(Problem.PR0_),
626  R0_user_(Problem.R0_user_),
627  PR0_user_(Problem.PR0_user_),
628  LP_(Problem.LP_),
629  RP_(Problem.RP_),
630  timerOp_(Problem.timerOp_),
631  timerPrec_(Problem.timerPrec_),
632  blocksize_(Problem.blocksize_),
633  num2Solve_(Problem.num2Solve_),
634  rhsIndex_(Problem.rhsIndex_),
635  lsNum_(Problem.lsNum_),
636  isSet_(Problem.isSet_),
637  isHermitian_(Problem.isHermitian_),
638  solutionUpdated_(Problem.solutionUpdated_),
639  label_(Problem.label_)
640  {
641  }
642 
643  template <class ScalarType, class MV, class OP>
644  void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
645  {
646  // Set new linear systems using the indices in index.
647  rhsIndex_ = index;
648 
649  // Compute the new block linear system.
650  // ( first clean up old linear system )
651  curB_ = Teuchos::null;
652  curX_ = Teuchos::null;
653 
654  // Create indices for the new linear system.
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;
664  validIdx++;
665  }
666  else {
667  iIndex[ivalidIdx] = i;
668  ivalidIdx++;
669  }
670  }
671  vldIndex.resize(validIdx);
672  newIndex.resize(validIdx);
673  iIndex.resize(ivalidIdx);
674  num2Solve_ = validIdx;
675 
676  // Create the new linear system using index
677  if (num2Solve_ != blocksize_) {
678  newIndex.resize(num2Solve_);
679  vldIndex.resize(num2Solve_);
680  //
681  // First create multivectors of blocksize.
682  // Fill the RHS with random vectors LHS with zero vectors.
683  curX_ = MVT::Clone( *X_, blocksize_ );
684  MVT::MvInit(*curX_);
685  Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
686  MVT::MvRandom(*tmpCurB);
687  //
688  // Now put in the part of B into curB
689  Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
690  MVT::SetBlock( *tptr, newIndex, *tmpCurB );
691  curB_ = tmpCurB;
692  //
693  // Now put in the part of X into curX
694  tptr = MVT::CloneView( *X_, vldIndex );
695  MVT::SetBlock( *tptr, newIndex, *curX_ );
696  //
697  solutionUpdated_ = false;
698  }
699  else {
700  curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
701  curB_ = MVT::CloneView( *B_, rhsIndex_ );
702  }
703  //
704  // Increment the number of linear systems that have been loaded into this object.
705  //
706  lsNum_++;
707  }
708 
709 
710  template <class ScalarType, class MV, class OP>
712  {
713  //
714  // We only need to copy the solutions back if the linear systems of
715  // interest are less than the block size.
716  //
717  if (num2Solve_ < blocksize_) {
718  //
719  // Get a view of the current solutions and correction std::vector.
720  //
721  int validIdx = 0;
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;
728  validIdx++;
729  }
730  }
731  Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
732  MVT::SetBlock( *tptr, vldIndex, *X_ );
733  }
734  //
735  // Clear the current vectors of this linear system so that any future calls
736  // to get the vectors for this system return null pointers.
737  //
738  curX_ = Teuchos::null;
739  curB_ = Teuchos::null;
740  rhsIndex_.resize(0);
741  }
742 
743 
744  template <class ScalarType, class MV, class OP>
748  bool updateLP,
749  ScalarType scale)
750  {
751  using Teuchos::RCP;
752  using Teuchos::null;
753 
754  RCP<MV> newSoln;
755  if (update.is_null())
756  { // The caller didn't supply an update vector, so we assume
757  // that the current solution curX_ is unchanged, and return a
758  // pointer to it.
759  newSoln = curX_;
760  }
761  else // the update vector is NOT null.
762  {
763  if (updateLP) // The caller wants us to update the linear problem.
764  {
765  if (RP_.is_null())
766  { // There is no right preconditioner.
767  // curX_ := curX_ + scale * update.
768  MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
769  }
770  else
771  { // There is indeed a right preconditioner, so apply it
772  // before computing the new solution.
773  RCP<MV> rightPrecUpdate =
774  MVT::Clone (*update, MVT::GetNumberVecs (*update));
775  applyRightPrec( *update, *rightPrecUpdate );
776  // curX_ := curX_ + scale * rightPrecUpdate.
777  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
778  }
779  solutionUpdated_ = true;
780  newSoln = curX_;
781  }
782  else
783  { // Rather than updating our stored current solution curX_,
784  // we make a copy and compute the new solution in the
785  // copy, without modifying curX_.
786  newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
787  if (RP_.is_null())
788  { // There is no right preconditioner.
789  // newSoln := curX_ + scale * update.
790  MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
791  }
792  else
793  { // There is indeed a right preconditioner, so apply it
794  // before computing the new solution.
795  RCP<MV> rightPrecUpdate =
796  MVT::Clone (*update, MVT::GetNumberVecs (*update));
797  applyRightPrec( *update, *rightPrecUpdate );
798  // newSoln := curX_ + scale * rightPrecUpdate.
799  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
800  }
801  }
802  }
803  return newSoln;
804  }
805 
806  template <class ScalarType, class MV, class OP>
807  void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
808  {
809  if (label != label_) {
810  label_ = label;
811  // Create new timers if they have already been created.
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 );
816 #endif
817  }
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 );
822 #endif
823  }
824  }
825  }
826 
827  template <class ScalarType, class MV, class OP>
828  bool
831  const Teuchos::RCP<const MV> &newB)
832  {
833  // Create timers if the haven't been created yet.
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 );
838 #endif
839  }
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 );
844 #endif
845  }
846 
847  // Set the linear system using the arguments newX and newB
848  if (newX != Teuchos::null)
849  X_ = newX;
850  if (newB != Teuchos::null)
851  B_ = newB;
852 
853  // Invalidate the current linear system indices and multivectors.
854  rhsIndex_.resize(0);
855  curX_ = Teuchos::null;
856  curB_ = Teuchos::null;
857 
858  // If we didn't set a matrix A, a left-hand side X, or a
859  // right-hand side B, then we didn't set the problem.
860  if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
861  isSet_ = false;
862  return isSet_;
863  }
864 
865  // Reset whether the solution has been updated. (We're just
866  // setting the problem now, so of course we haven't updated the
867  // solution yet.)
868  solutionUpdated_ = false;
869 
870  // Compute the initial residuals.
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_ ) );
874  }
875  computeCurrResVec( &*R0_, &*X_, &*B_ );
876 
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_ ) );
880  }
881  applyLeftPrec( *R0_, *PR0_ );
882  }
883  else {
884  PR0_ = R0_;
885  }
886  }
887  else { // User specified the residuals
888  // If the user did not specify the right sized residual, create one and set R0_user_ to null.
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;
893  R0_ = helper;
894  }
895 
896  if (LP_!=Teuchos::null) {
897  // If the user provided preconditioned residual is the wrong size or pointing at
898  // the wrong object, create one and set the PR0_user_ to 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_ ) );
902 
903  // Get the initial residual from getInitResVec because R0_user_ may be null from above.
904  applyLeftPrec( *getInitResVec(), *helper );
905  PR0_user_ = Teuchos::null;
906  PR0_ = helper;
907  }
908  }
909  else {
910  // The preconditioned initial residual vector is the residual vector.
911  // NOTE: R0_user_ could be null if incompatible.
912  if (R0_user_!=Teuchos::null)
913  {
914  PR0_user_ = R0_user_;
915  }
916  else
917  {
918  PR0_user_ = Teuchos::null;
919  PR0_ = R0_;
920  }
921  }
922  }
923 
924  // The problem has been set and is ready for use.
925  isSet_ = true;
926 
927  // Return isSet.
928  return isSet_;
929  }
930 
931  template <class ScalarType, class MV, class OP>
933  {
934  if(Teuchos::nonnull(R0_user_)) {
935  return R0_user_;
936  }
937  return(R0_);
938  }
939 
940  template <class ScalarType, class MV, class OP>
942  {
943  if(Teuchos::nonnull(PR0_user_)) {
944  return PR0_user_;
945  }
946  return(PR0_);
947  }
948 
949  template <class ScalarType, class MV, class OP>
951  {
952  if (isSet_) {
953  return curX_;
954  }
955  else {
956  return Teuchos::null;
957  }
958  }
959 
960  template <class ScalarType, class MV, class OP>
962  {
963  if (isSet_) {
964  return curB_;
965  }
966  else {
967  return Teuchos::null;
968  }
969  }
970 
971  template <class ScalarType, class MV, class OP>
972  void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
973  {
974  using Teuchos::null;
975  using Teuchos::RCP;
976 
977  const bool leftPrec = LP_ != null;
978  const bool rightPrec = RP_ != null;
979 
980  // We only need a temporary vector for intermediate results if
981  // there is a left or right preconditioner. We really should just
982  // keep this temporary vector around instead of allocating it each
983  // time.
984  RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
985 
986  //
987  // No preconditioning.
988  //
989  if (!leftPrec && !rightPrec) {
990  applyOp( x, y );
991  }
992  //
993  // Preconditioning is being done on both sides
994  //
995  else if( leftPrec && rightPrec )
996  {
997  applyRightPrec( x, y );
998  applyOp( y, *ytemp );
999  applyLeftPrec( *ytemp, y );
1000  }
1001  //
1002  // Preconditioning is only being done on the left side
1003  //
1004  else if( leftPrec )
1005  {
1006  applyOp( x, *ytemp );
1007  applyLeftPrec( *ytemp, y );
1008  }
1009  //
1010  // Preconditioning is only being done on the right side
1011  //
1012  else
1013  {
1014  applyRightPrec( x, *ytemp );
1015  applyOp( *ytemp, y );
1016  }
1017  }
1018 
1019  template <class ScalarType, class MV, class OP>
1020  void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
1021  if (A_.get()) {
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1023  Teuchos::TimeMonitor OpTimer(*timerOp_);
1024 #endif
1025  OPT::Apply( *A_,x, y);
1026  }
1027  else {
1028  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1030  }
1031  }
1032 
1033  template <class ScalarType, class MV, class OP>
1034  void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
1035  if (LP_!=Teuchos::null) {
1036 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1037  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1038 #endif
1039  OPT::Apply( *LP_,x, y);
1040  }
1041  else {
1042  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1044  }
1045  }
1046 
1047  template <class ScalarType, class MV, class OP>
1048  void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
1049  if (RP_!=Teuchos::null) {
1050 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1051  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1052 #endif
1053  OPT::Apply( *RP_,x, y);
1054  }
1055  else {
1056  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1058  }
1059  }
1060 
1061  template <class ScalarType, class MV, class OP>
1062  void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
1063 
1064  if (R) {
1065  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1066  {
1067  if (LP_!=Teuchos::null)
1068  {
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 );
1073  }
1074  else
1075  {
1076  applyOp( *X, *R );
1077  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1078  }
1079  }
1080  else {
1081  // The solution and right-hand side may not be specified, check and use which ones exist.
1082  Teuchos::RCP<const MV> localB, localX;
1083  if (B)
1084  localB = Teuchos::rcp( B, false );
1085  else
1086  localB = curB_;
1087 
1088  if (X)
1089  localX = Teuchos::rcp( X, false );
1090  else
1091  localX = curX_;
1092 
1093  if (LP_!=Teuchos::null)
1094  {
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 );
1099  }
1100  else
1101  {
1102  applyOp( *localX, *R );
1103  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1104  }
1105  }
1106  }
1107  }
1108 
1109 
1110  template <class ScalarType, class MV, class OP>
1111  void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
1112 
1113  if (R) {
1114  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1115  {
1116  applyOp( *X, *R );
1117  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1118  }
1119  else {
1120  // The solution and right-hand side may not be specified, check and use which ones exist.
1121  Teuchos::RCP<const MV> localB, localX;
1122  if (B)
1123  localB = Teuchos::rcp( B, false );
1124  else
1125  localB = curB_;
1126 
1127  if (X)
1128  localX = Teuchos::rcp( X, false );
1129  else
1130  localX = curX_;
1131 
1132  applyOp( *localX, *R );
1133  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1134  }
1135  }
1136  }
1137 
1138 } // end Belos namespace
1139 
1140 #endif /* BELOS_LINEAR_PROBLEM_HPP */
1141 
1142 
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?
Teuchos::RCP< MV > getLHS() const
A pointer to the left-hand side X.
Exception thrown to signal error with the Belos::LinearProblem object.
virtual void apply(const MV &x, MV &y) const
Apply the composite operator of this linear problem to x, returning y.
void setHermitian(bool isSym=true)
Tell the linear problem that the (preconditioned) operator is Hermitian.
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.
Teuchos::RCP< const MV > B_
Right-hand side of linear system.
bool isProblemSet() const
Whether the problem has been set.
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.
bool nonnull(const std::shared_ptr< T > &p)
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.
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< const MV > R0_user_
User-defined initial residual of the linear system.
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::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
Teuchos::RCP< MV > curX_
Current solution vector of the linear system.
Teuchos::RCP< const OP > LP_
Left preconditioning operator of linear system.
LinearProblem(void)
Default constructor.
void setOperator(const Teuchos::RCP< const OP > &A)
Set the operator A of the linear problem .
Teuchos::RCP< const MV > getRHS() const
A pointer to the right-hand side B.
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.
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. ...
void setRHS(const Teuchos::RCP< const MV > &B)
Set right-hand-side B of linear problem .
bool isRightPrec() const
Whether the linear system is being preconditioned on the right.
Traits class which defines basic operations on multivectors.
virtual void applyOp(const MV &x, MV &y) const
Apply ONLY the operator to x, returning y.
Teuchos::RCP< const OP > getOperator() const
A pointer to the (unpreconditioned) operator A.
int blocksize_
Current block size of linear system.
const std::vector< int > & getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSolutionUpdated() const
Has the current approximate solution been updated?
bool is_null(const RCP< T > &p)
void setLeftPrec(const Teuchos::RCP< const OP > &LP)
Set left preconditioner (LP) of linear problem .
A linear system to solve, and its associated information.
int getLSNumber() const
The number of linear systems that have been set.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
The timers for this object.
MultiVecTraits< ScalarType, MV > MVT
bool isLeftPrec() const
Whether the linear system is being preconditioned on the left.
OperatorTraits< ScalarType, MV, OP > OPT
std::string label_
Linear problem label that prefixes the timer labels.
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
virtual void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
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.
Teuchos::RCP< const MV > PR0_user_
User-defined preconditioned initial residual of the linear system.
Teuchos::RCP< MV > R0_
Initial residual of the linear system.
Teuchos::RCP< const OP > A_
Operator of linear system.
virtual void applyLeftPrec(const MV &x, MV &y) const
Apply ONLY the left preconditioner to x, returning y.
Teuchos::RCP< MV > X_
Solution vector of linear system.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
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. ...
LinearProblemError(const std::string &what_arg)
bool solutionUpdated_
Has the current approximate solution been updated?
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
virtual void setCurrLS()
Tell the linear problem that the solver is finished with the current linear system.
Teuchos::RCP< MV > PR0_
Preconditioned initial residual of the linear system.
Class which defines basic traits for the operator type.
bool isHermitian() const
Whether the (preconditioned) operator is Hermitian.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
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.
Teuchos::RCP< const OP > RP_
Right preconditioning operator of linear system.
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 ( &lt;= blocksize_ )
void setInitPrecResVec(const Teuchos::RCP< const MV > &PR0)
Set the user-defined preconditioned residual of linear problem .
Teuchos::RCP< Teuchos::Time > timerPrec_