Belos  Version of the Day
 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  void setCurrLS ();
186 
196  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 
255  updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
256  bool updateLP = false,
257  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
258 
276  Teuchos::RCP<MV> updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
277  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
278  { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
279 
281 
283 
284 
310  bool
311  setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null,
312  const Teuchos::RCP<const MV> &newB = Teuchos::null);
313 
315 
317 
318 
320  Teuchos::RCP<const OP> getOperator() const { return(A_); }
321 
323  Teuchos::RCP<MV> getLHS() const { return(X_); }
324 
326  Teuchos::RCP<const MV> getRHS() const { return(B_); }
327 
330 
336 
352 
368 
370  Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
371 
373  Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
374 
396  const std::vector<int> getLSIndex() const { return(rhsIndex_); }
397 
404  int getLSNumber() const { return(lsNum_); }
405 
413  return Teuchos::tuple(timerOp_,timerPrec_);
414  }
415 
416 
418 
420 
421 
429  bool isSolutionUpdated() const { return(solutionUpdated_); }
430 
432  bool isProblemSet() const { return(isSet_); }
433 
439  bool isHermitian() const { return(isHermitian_); }
440 
442  bool isLeftPrec() const { return(LP_!=Teuchos::null); }
443 
445  bool isRightPrec() const { return(RP_!=Teuchos::null); }
446 
448 
450 
451 
453 
460  void apply( const MV& x, MV& y ) const;
461 
469  void applyOp( const MV& x, MV& y ) const;
470 
477  void applyLeftPrec( const MV& x, MV& y ) const;
478 
485  void applyRightPrec( const MV& x, MV& y ) const;
486 
488 
492  void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
493 
495 
499  void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
500 
502 
503  private:
504 
507 
509  Teuchos::RCP<MV> X_;
510 
512  Teuchos::RCP<MV> curX_;
513 
516 
519 
521  Teuchos::RCP<MV> R0_;
522 
524  Teuchos::RCP<MV> PR0_;
525 
527  Teuchos::RCP<const MV> R0_user_;
528 
530  Teuchos::RCP<const MV> PR0_user_;
531 
534 
537 
539  mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
540 
542  int blocksize_;
543 
545  int num2Solve_;
546 
548  std::vector<int> rhsIndex_;
549 
551  int lsNum_;
552 
554 
555 
557  bool isSet_;
558 
561  bool isHermitian_;
562 
564  bool solutionUpdated_;
565 
567 
569  std::string label_;
570 
571  typedef MultiVecTraits<ScalarType,MV> MVT;
573  };
574 
575  //--------------------------------------------
576  // Constructor Implementations
577  //--------------------------------------------
578 
579  template <class ScalarType, class MV, class OP>
581  blocksize_(0),
582  num2Solve_(0),
583  rhsIndex_(0),
584  lsNum_(0),
585  isSet_(false),
586  isHermitian_(false),
587  solutionUpdated_(false),
588  label_("Belos")
589  {
590  }
591 
592  template <class ScalarType, class MV, class OP>
594  const Teuchos::RCP<MV> &X,
595  const Teuchos::RCP<const MV> &B
596  ) :
597  A_(A),
598  X_(X),
599  B_(B),
600  blocksize_(0),
601  num2Solve_(0),
602  rhsIndex_(0),
603  lsNum_(0),
604  isSet_(false),
605  isHermitian_(false),
606  solutionUpdated_(false),
607  label_("Belos")
608  {
609  }
610 
611  template <class ScalarType, class MV, class OP>
613  A_(Problem.A_),
614  X_(Problem.X_),
615  curX_(Problem.curX_),
616  B_(Problem.B_),
617  curB_(Problem.curB_),
618  R0_(Problem.R0_),
619  PR0_(Problem.PR0_),
620  R0_user_(Problem.R0_user_),
621  PR0_user_(Problem.PR0_user_),
622  LP_(Problem.LP_),
623  RP_(Problem.RP_),
624  timerOp_(Problem.timerOp_),
625  timerPrec_(Problem.timerPrec_),
626  blocksize_(Problem.blocksize_),
627  num2Solve_(Problem.num2Solve_),
628  rhsIndex_(Problem.rhsIndex_),
629  lsNum_(Problem.lsNum_),
630  isSet_(Problem.isSet_),
631  isHermitian_(Problem.isHermitian_),
632  solutionUpdated_(Problem.solutionUpdated_),
633  label_(Problem.label_)
634  {
635  }
636 
637  template <class ScalarType, class MV, class OP>
639  {}
640 
641  template <class ScalarType, class MV, class OP>
642  void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
643  {
644  // Set new linear systems using the indices in index.
645  rhsIndex_ = index;
646 
647  // Compute the new block linear system.
648  // ( first clean up old linear system )
649  curB_ = Teuchos::null;
650  curX_ = Teuchos::null;
651 
652  // Create indices for the new linear system.
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;
662  validIdx++;
663  }
664  else {
665  iIndex[ivalidIdx] = i;
666  ivalidIdx++;
667  }
668  }
669  vldIndex.resize(validIdx);
670  newIndex.resize(validIdx);
671  iIndex.resize(ivalidIdx);
672  num2Solve_ = validIdx;
673 
674  // Create the new linear system using index
675  if (num2Solve_ != blocksize_) {
676  newIndex.resize(num2Solve_);
677  vldIndex.resize(num2Solve_);
678  //
679  // First create multivectors of blocksize.
680  // Fill the RHS with random vectors LHS with zero vectors.
681  curX_ = MVT::Clone( *X_, blocksize_ );
682  MVT::MvInit(*curX_);
683  Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
684  MVT::MvRandom(*tmpCurB);
685  //
686  // Now put in the part of B into curB
687  Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
688  MVT::SetBlock( *tptr, newIndex, *tmpCurB );
689  curB_ = tmpCurB;
690  //
691  // Now put in the part of X into curX
692  tptr = MVT::CloneView( *X_, vldIndex );
693  MVT::SetBlock( *tptr, newIndex, *curX_ );
694  //
695  solutionUpdated_ = false;
696  }
697  else {
698  curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
699  curB_ = MVT::CloneView( *B_, rhsIndex_ );
700  }
701  //
702  // Increment the number of linear systems that have been loaded into this object.
703  //
704  lsNum_++;
705  }
706 
707 
708  template <class ScalarType, class MV, class OP>
710  {
711  //
712  // We only need to copy the solutions back if the linear systems of
713  // interest are less than the block size.
714  //
715  if (num2Solve_ < blocksize_) {
716  //
717  // Get a view of the current solutions and correction std::vector.
718  //
719  int validIdx = 0;
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;
726  validIdx++;
727  }
728  }
729  Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
730  MVT::SetBlock( *tptr, vldIndex, *X_ );
731  }
732  //
733  // Clear the current vectors of this linear system so that any future calls
734  // to get the vectors for this system return null pointers.
735  //
736  curX_ = Teuchos::null;
737  curB_ = Teuchos::null;
738  rhsIndex_.resize(0);
739  }
740 
741 
742  template <class ScalarType, class MV, class OP>
746  bool updateLP,
747  ScalarType scale)
748  {
749  using Teuchos::RCP;
750  using Teuchos::null;
751 
752  RCP<MV> newSoln;
753  if (update.is_null())
754  { // The caller didn't supply an update vector, so we assume
755  // that the current solution curX_ is unchanged, and return a
756  // pointer to it.
757  newSoln = curX_;
758  }
759  else // the update vector is NOT null.
760  {
761  if (updateLP) // The caller wants us to update the linear problem.
762  {
763  if (RP_.is_null())
764  { // There is no right preconditioner.
765  // curX_ := curX_ + scale * update.
766  MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
767  }
768  else
769  { // There is indeed a right preconditioner, so apply it
770  // before computing the new solution.
771  RCP<MV> rightPrecUpdate =
772  MVT::Clone (*update, MVT::GetNumberVecs (*update));
773  {
774 #ifdef BELOS_TEUCHOS_TIME_MONITOR
775  Teuchos::TimeMonitor PrecTimer (*timerPrec_);
776 #endif
777  OPT::Apply( *RP_, *update, *rightPrecUpdate );
778  }
779  // curX_ := curX_ + scale * rightPrecUpdate.
780  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
781  }
782  solutionUpdated_ = true;
783  newSoln = curX_;
784  }
785  else
786  { // Rather than updating our stored current solution curX_,
787  // we make a copy and compute the new solution in the
788  // copy, without modifying curX_.
789  newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
790  if (RP_.is_null())
791  { // There is no right preconditioner.
792  // newSoln := curX_ + scale * update.
793  MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
794  }
795  else
796  { // There is indeed a right preconditioner, so apply it
797  // before computing the new solution.
798  RCP<MV> rightPrecUpdate =
799  MVT::Clone (*update, MVT::GetNumberVecs (*update));
800  {
801 #ifdef BELOS_TEUCHOS_TIME_MONITOR
802  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
803 #endif
804  OPT::Apply( *RP_, *update, *rightPrecUpdate );
805  }
806  // newSoln := curX_ + scale * rightPrecUpdate.
807  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
808  }
809  }
810  }
811  return newSoln;
812  }
813 
814  template <class ScalarType, class MV, class OP>
815  void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
816  {
817  if (label != label_) {
818  label_ = label;
819  // Create new timers if they have already been created.
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 );
824 #endif
825  }
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 );
830 #endif
831  }
832  }
833  }
834 
835  template <class ScalarType, class MV, class OP>
836  bool
839  const Teuchos::RCP<const MV> &newB)
840  {
841  // Create timers if the haven't been created yet.
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 );
846 #endif
847  }
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 );
852 #endif
853  }
854 
855  // Set the linear system using the arguments newX and newB
856  if (newX != Teuchos::null)
857  X_ = newX;
858  if (newB != Teuchos::null)
859  B_ = newB;
860 
861  // Invalidate the current linear system indices and multivectors.
862  rhsIndex_.resize(0);
863  curX_ = Teuchos::null;
864  curB_ = Teuchos::null;
865 
866  // If we didn't set a matrix A, a left-hand side X, or a
867  // right-hand side B, then we didn't set the problem.
868  if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
869  isSet_ = false;
870  return isSet_;
871  }
872 
873  // Reset whether the solution has been updated. (We're just
874  // setting the problem now, so of course we haven't updated the
875  // solution yet.)
876  solutionUpdated_ = false;
877 
878  // Compute the initial residuals.
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_ ) );
882  }
883  computeCurrResVec( &*R0_, &*X_, &*B_ );
884 
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_ ) );
888  }
889  {
890 #ifdef BELOS_TEUCHOS_TIME_MONITOR
891  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
892 #endif
893  OPT::Apply( *LP_, *R0_, *PR0_ );
894  }
895  }
896  else {
897  PR0_ = R0_;
898  }
899  }
900  else { // User specified the residuals
901  // If the user did not specify the right sized residual, create one and set R0_user_ to null.
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;
906  R0_ = helper;
907  }
908 
909  if (LP_!=Teuchos::null) {
910  // If the user provided preconditioned residual is the wrong size or pointing at
911  // the wrong object, create one and set the PR0_user_ to 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_ ) );
915  {
916  // Get the initial residual from getInitResVec because R0_user_ may be null from above.
917 #ifdef BELOS_TEUCHOS_TIME_MONITOR
918  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
919 #endif
920  OPT::Apply( *LP_, *getInitResVec(), *helper );
921  }
922  PR0_user_ = Teuchos::null;
923  PR0_ = helper;
924  }
925  }
926  else {
927  // The preconditioned initial residual vector is the residual vector.
928  // NOTE: R0_user_ could be null if incompatible.
929  if (R0_user_!=Teuchos::null)
930  {
931  PR0_user_ = R0_user_;
932  }
933  else
934  {
935  PR0_user_ = Teuchos::null;
936  PR0_ = R0_;
937  }
938  }
939  }
940 
941  // The problem has been set and is ready for use.
942  isSet_ = true;
943 
944  // Return isSet.
945  return isSet_;
946  }
947 
948  template <class ScalarType, class MV, class OP>
950  {
951  if(Teuchos::nonnull(R0_user_)) {
952  return R0_user_;
953  }
954  return(R0_);
955  }
956 
957  template <class ScalarType, class MV, class OP>
959  {
960  if(Teuchos::nonnull(PR0_user_)) {
961  return PR0_user_;
962  }
963  return(PR0_);
964  }
965 
966  template <class ScalarType, class MV, class OP>
968  {
969  if (isSet_) {
970  return curX_;
971  }
972  else {
973  return Teuchos::null;
974  }
975  }
976 
977  template <class ScalarType, class MV, class OP>
979  {
980  if (isSet_) {
981  return curB_;
982  }
983  else {
984  return Teuchos::null;
985  }
986  }
987 
988  template <class ScalarType, class MV, class OP>
989  void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
990  {
991  using Teuchos::null;
992  using Teuchos::RCP;
993 
994  const bool leftPrec = LP_ != null;
995  const bool rightPrec = RP_ != null;
996 
997  // We only need a temporary vector for intermediate results if
998  // there is a left or right preconditioner. We really should just
999  // keep this temporary vector around instead of allocating it each
1000  // time.
1001  RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
1002 
1003  //
1004  // No preconditioning.
1005  //
1006  if (!leftPrec && !rightPrec){
1007 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1008  Teuchos::TimeMonitor OpTimer(*timerOp_);
1009 #endif
1010  OPT::Apply( *A_, x, y );
1011  }
1012  //
1013  // Preconditioning is being done on both sides
1014  //
1015  else if( leftPrec && rightPrec )
1016  {
1017  {
1018 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1019  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1020 #endif
1021  OPT::Apply( *RP_, x, y );
1022  }
1023  {
1024 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1025  Teuchos::TimeMonitor OpTimer(*timerOp_);
1026 #endif
1027  OPT::Apply( *A_, y, *ytemp );
1028  }
1029  {
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1031  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1032 #endif
1033  OPT::Apply( *LP_, *ytemp, y );
1034  }
1035  }
1036  //
1037  // Preconditioning is only being done on the left side
1038  //
1039  else if( leftPrec )
1040  {
1041  {
1042 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1043  Teuchos::TimeMonitor OpTimer(*timerOp_);
1044 #endif
1045  OPT::Apply( *A_, x, *ytemp );
1046  }
1047  {
1048 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1049  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1050 #endif
1051  OPT::Apply( *LP_, *ytemp, y );
1052  }
1053  }
1054  //
1055  // Preconditioning is only being done on the right side
1056  //
1057  else
1058  {
1059  {
1060 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1061  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1062 #endif
1063  OPT::Apply( *RP_, x, *ytemp );
1064  }
1065  {
1066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1067  Teuchos::TimeMonitor OpTimer(*timerOp_);
1068 #endif
1069  OPT::Apply( *A_, *ytemp, y );
1070  }
1071  }
1072  }
1073 
1074  template <class ScalarType, class MV, class OP>
1075  void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
1076  if (A_.get()) {
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1078  Teuchos::TimeMonitor OpTimer(*timerOp_);
1079 #endif
1080  OPT::Apply( *A_,x, y);
1081  }
1082  else {
1083  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1085  }
1086  }
1087 
1088  template <class ScalarType, class MV, class OP>
1089  void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
1090  if (LP_!=Teuchos::null) {
1091 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1092  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1093 #endif
1094  return ( OPT::Apply( *LP_,x, y) );
1095  }
1096  else {
1097  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1099  }
1100  }
1101 
1102  template <class ScalarType, class MV, class OP>
1103  void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
1104  if (RP_!=Teuchos::null) {
1105 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1106  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1107 #endif
1108  return ( OPT::Apply( *RP_,x, y) );
1109  }
1110  else {
1111  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1113  }
1114  }
1115 
1116  template <class ScalarType, class MV, class OP>
1117  void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
1118 
1119  if (R) {
1120  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1121  {
1122  if (LP_!=Teuchos::null)
1123  {
1124  Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1125  {
1126 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1127  Teuchos::TimeMonitor OpTimer(*timerOp_);
1128 #endif
1129  OPT::Apply( *A_, *X, *R_temp );
1130  }
1131  MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1132  {
1133 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1134  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1135 #endif
1136  OPT::Apply( *LP_, *R_temp, *R );
1137  }
1138  }
1139  else
1140  {
1141  {
1142 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1143  Teuchos::TimeMonitor OpTimer(*timerOp_);
1144 #endif
1145  OPT::Apply( *A_, *X, *R );
1146  }
1147  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1148  }
1149  }
1150  else {
1151  // The solution and right-hand side may not be specified, check and use which ones exist.
1152  Teuchos::RCP<const MV> localB, localX;
1153  if (B)
1154  localB = Teuchos::rcp( B, false );
1155  else
1156  localB = curB_;
1157 
1158  if (X)
1159  localX = Teuchos::rcp( X, false );
1160  else
1161  localX = curX_;
1162 
1163  if (LP_!=Teuchos::null)
1164  {
1165  Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1166  {
1167 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1168  Teuchos::TimeMonitor OpTimer(*timerOp_);
1169 #endif
1170  OPT::Apply( *A_, *localX, *R_temp );
1171  }
1172  MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1173  {
1174 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1175  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1176 #endif
1177  OPT::Apply( *LP_, *R_temp, *R );
1178  }
1179  }
1180  else
1181  {
1182  {
1183 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1184  Teuchos::TimeMonitor OpTimer(*timerOp_);
1185 #endif
1186  OPT::Apply( *A_, *localX, *R );
1187  }
1188  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1189  }
1190  }
1191  }
1192  }
1193 
1194 
1195  template <class ScalarType, class MV, class OP>
1196  void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
1197 
1198  if (R) {
1199  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1200  {
1201  {
1202 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1203  Teuchos::TimeMonitor OpTimer(*timerOp_);
1204 #endif
1205  OPT::Apply( *A_, *X, *R );
1206  }
1207  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1208  }
1209  else {
1210  // The solution and right-hand side may not be specified, check and use which ones exist.
1211  Teuchos::RCP<const MV> localB, localX;
1212  if (B)
1213  localB = Teuchos::rcp( B, false );
1214  else
1215  localB = curB_;
1216 
1217  if (X)
1218  localX = Teuchos::rcp( X, false );
1219  else
1220  localX = curX_;
1221 
1222  {
1223 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1224  Teuchos::TimeMonitor OpTimer(*timerOp_);
1225 #endif
1226  OPT::Apply( *A_, *localX, *R );
1227  }
1228  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1229  }
1230  }
1231  }
1232 
1233 } // end Belos namespace
1234 
1235 #endif /* BELOS_LINEAR_PROBLEM_HPP */
1236 
1237 
Teuchos::RCP< MV > getLHS() const
A pointer to the left-hand side X.
Exception thrown to signal error with the Belos::LinearProblem object.
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.
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)
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)
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::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
LinearProblem(void)
Default constructor.
void setOperator(const Teuchos::RCP< const OP > &A)
Set the operator A of the linear problem .
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 > getRHS() const
A pointer to the right-hand side B.
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.
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.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSolutionUpdated() const
Has the current approximate solution been updated?
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.
bool isLeftPrec() const
Whether the linear system is being preconditioned on the left.
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
void applyRightPrec(const MV &x, MV &y) const
Apply ONLY the right preconditioner to x, returning y.
void setLabel(const std::string &label)
Set the label prefix used by the timers in this object.
void applyLeftPrec(const MV &x, MV &y) const
Apply ONLY the left preconditioner to x, returning y.
Teuchos::RCP< const OP > getLeftPrec() const
Get a pointer to the left preconditioner.
const std::vector< int > getLSIndex() const
(Zero-based) indices of the linear system(s) currently being solved.
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)
Teuchos::RCP< const OP > getRightPrec() const
Get a pointer to the right preconditioner.
void setCurrLS()
Tell the linear problem that the solver is finished with the current 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 .
void setInitPrecResVec(const Teuchos::RCP< const MV > &PR0)
Set the user-defined preconditioned residual of linear problem .
bool is_null() const

Generated on Fri Jun 5 2020 10:20:41 for Belos by doxygen 1.8.5