Belos  Version of the Day
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 
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 setLeftPrec(const Teuchos::RCP<const OP> &LP) { LP_ = LP; }
152 
156  void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }
157 
166  void setCurrLS ();
167 
177  void setLSIndex (const std::vector<int>& index);
178 
189  void setHermitian( bool isSym = true ) { isHermitian_ = isSym; }
190 
197  void setLabel (const std::string& label);
198 
235  Teuchos::RCP<MV>
236  updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
237  bool updateLP = false,
238  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());
239 
257  Teuchos::RCP<MV> updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
258  ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
259  { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
260 
262 
264 
265 
291  bool
292  setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null,
293  const Teuchos::RCP<const MV> &newB = Teuchos::null);
294 
296 
298 
299 
301  Teuchos::RCP<const OP> getOperator() const { return(A_); }
302 
304  Teuchos::RCP<MV> getLHS() const { return(X_); }
305 
307  Teuchos::RCP<const MV> getRHS() const { return(B_); }
308 
310  Teuchos::RCP<const MV> getInitResVec() const { return(R0_); }
311 
316  Teuchos::RCP<const MV> getInitPrecResVec() const { return(PR0_); }
317 
332  Teuchos::RCP<MV> getCurrLHSVec();
333 
348  Teuchos::RCP<const MV> getCurrRHSVec();
349 
351  Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
352 
354  Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
355 
377  const std::vector<int> getLSIndex() const { return(rhsIndex_); }
378 
385  int getLSNumber() const { return(lsNum_); }
386 
393  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
394  return Teuchos::tuple(timerOp_,timerPrec_);
395  }
396 
397 
399 
401 
402 
410  bool isSolutionUpdated() const { return(solutionUpdated_); }
411 
413  bool isProblemSet() const { return(isSet_); }
414 
420  bool isHermitian() const { return(isHermitian_); }
421 
423  bool isLeftPrec() const { return(LP_!=Teuchos::null); }
424 
426  bool isRightPrec() const { return(RP_!=Teuchos::null); }
427 
429 
431 
432 
434 
441  void apply( const MV& x, MV& y ) const;
442 
450  void applyOp( const MV& x, MV& y ) const;
451 
458  void applyLeftPrec( const MV& x, MV& y ) const;
459 
466  void applyRightPrec( const MV& x, MV& y ) const;
467 
469 
473  void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
474 
476 
480  void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
481 
483 
484  private:
485 
487  Teuchos::RCP<const OP> A_;
488 
490  Teuchos::RCP<MV> X_;
491 
493  Teuchos::RCP<MV> curX_;
494 
496  Teuchos::RCP<const MV> B_;
497 
499  Teuchos::RCP<const MV> curB_;
500 
502  Teuchos::RCP<MV> R0_;
503 
505  Teuchos::RCP<MV> PR0_;
506 
508  Teuchos::RCP<const OP> LP_;
509 
511  Teuchos::RCP<const OP> RP_;
512 
514  mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
515 
517  int blocksize_;
518 
520  int num2Solve_;
521 
523  std::vector<int> rhsIndex_;
524 
526  int lsNum_;
527 
529 
530 
532  bool Left_Scale_;
533 
535  bool Right_Scale_;
536 
538  bool isSet_;
539 
542  bool isHermitian_;
543 
545  bool solutionUpdated_;
546 
548 
550  std::string label_;
551 
554  };
555 
556  //--------------------------------------------
557  // Constructor Implementations
558  //--------------------------------------------
559 
560  template <class ScalarType, class MV, class OP>
562  blocksize_(0),
563  num2Solve_(0),
564  rhsIndex_(0),
565  lsNum_(0),
566  Left_Scale_(false),
567  Right_Scale_(false),
568  isSet_(false),
569  isHermitian_(false),
570  solutionUpdated_(false),
571  label_("Belos")
572  {
573  }
574 
575  template <class ScalarType, class MV, class OP>
576  LinearProblem<ScalarType,MV,OP>::LinearProblem(const Teuchos::RCP<const OP> &A,
577  const Teuchos::RCP<MV> &X,
578  const Teuchos::RCP<const MV> &B
579  ) :
580  A_(A),
581  X_(X),
582  B_(B),
583  blocksize_(0),
584  num2Solve_(0),
585  rhsIndex_(0),
586  lsNum_(0),
587  Left_Scale_(false),
588  Right_Scale_(false),
589  isSet_(false),
590  isHermitian_(false),
591  solutionUpdated_(false),
592  label_("Belos")
593  {
594  }
595 
596  template <class ScalarType, class MV, class OP>
598  A_(Problem.A_),
599  X_(Problem.X_),
600  curX_(Problem.curX_),
601  B_(Problem.B_),
602  curB_(Problem.curB_),
603  R0_(Problem.R0_),
604  PR0_(Problem.PR0_),
605  LP_(Problem.LP_),
606  RP_(Problem.RP_),
607  timerOp_(Problem.timerOp_),
608  timerPrec_(Problem.timerPrec_),
609  blocksize_(Problem.blocksize_),
610  num2Solve_(Problem.num2Solve_),
611  rhsIndex_(Problem.rhsIndex_),
612  lsNum_(Problem.lsNum_),
613  Left_Scale_(Problem.Left_Scale_),
614  Right_Scale_(Problem.Right_Scale_),
615  isSet_(Problem.isSet_),
616  isHermitian_(Problem.isHermitian_),
617  solutionUpdated_(Problem.solutionUpdated_),
618  label_(Problem.label_)
619  {
620  }
621 
622  template <class ScalarType, class MV, class OP>
624  {}
625 
626  template <class ScalarType, class MV, class OP>
627  void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
628  {
629  // Set new linear systems using the indices in index.
630  rhsIndex_ = index;
631 
632  // Compute the new block linear system.
633  // ( first clean up old linear system )
634  curB_ = Teuchos::null;
635  curX_ = Teuchos::null;
636 
637  // Create indices for the new linear system.
638  int validIdx = 0, ivalidIdx = 0;
639  blocksize_ = rhsIndex_.size();
640  std::vector<int> vldIndex( blocksize_ );
641  std::vector<int> newIndex( blocksize_ );
642  std::vector<int> iIndex( blocksize_ );
643  for (int i=0; i<blocksize_; ++i) {
644  if (rhsIndex_[i] > -1) {
645  vldIndex[validIdx] = rhsIndex_[i];
646  newIndex[validIdx] = i;
647  validIdx++;
648  }
649  else {
650  iIndex[ivalidIdx] = i;
651  ivalidIdx++;
652  }
653  }
654  vldIndex.resize(validIdx);
655  newIndex.resize(validIdx);
656  iIndex.resize(ivalidIdx);
657  num2Solve_ = validIdx;
658 
659  // Create the new linear system using index
660  if (num2Solve_ != blocksize_) {
661  newIndex.resize(num2Solve_);
662  vldIndex.resize(num2Solve_);
663  //
664  // First create multivectors of blocksize.
665  // Fill the RHS with random vectors LHS with zero vectors.
666  curX_ = MVT::Clone( *X_, blocksize_ );
667  MVT::MvInit(*curX_);
668  Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
669  MVT::MvRandom(*tmpCurB);
670  //
671  // Now put in the part of B into curB
672  Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
673  MVT::SetBlock( *tptr, newIndex, *tmpCurB );
674  curB_ = tmpCurB;
675  //
676  // Now put in the part of X into curX
677  tptr = MVT::CloneView( *X_, vldIndex );
678  MVT::SetBlock( *tptr, newIndex, *curX_ );
679  //
680  solutionUpdated_ = false;
681  }
682  else {
683  curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
684  curB_ = MVT::CloneView( *B_, rhsIndex_ );
685  }
686  //
687  // Increment the number of linear systems that have been loaded into this object.
688  //
689  lsNum_++;
690  }
691 
692 
693  template <class ScalarType, class MV, class OP>
695  {
696  //
697  // We only need to copy the solutions back if the linear systems of
698  // interest are less than the block size.
699  //
700  if (num2Solve_ < blocksize_) {
701  //
702  // Get a view of the current solutions and correction std::vector.
703  //
704  int validIdx = 0;
705  std::vector<int> newIndex( num2Solve_ );
706  std::vector<int> vldIndex( num2Solve_ );
707  for (int i=0; i<blocksize_; ++i) {
708  if ( rhsIndex_[i] > -1 ) {
709  vldIndex[validIdx] = rhsIndex_[i];
710  newIndex[validIdx] = i;
711  validIdx++;
712  }
713  }
714  Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
715  MVT::SetBlock( *tptr, vldIndex, *X_ );
716  }
717  //
718  // Clear the current vectors of this linear system so that any future calls
719  // to get the vectors for this system return null pointers.
720  //
721  curX_ = Teuchos::null;
722  curB_ = Teuchos::null;
723  rhsIndex_.resize(0);
724  }
725 
726 
727  template <class ScalarType, class MV, class OP>
728  Teuchos::RCP<MV>
730  updateSolution (const Teuchos::RCP<MV>& update,
731  bool updateLP,
732  ScalarType scale)
733  {
734  using Teuchos::RCP;
735  using Teuchos::null;
736 
737  RCP<MV> newSoln;
738  if (update.is_null())
739  { // The caller didn't supply an update vector, so we assume
740  // that the current solution curX_ is unchanged, and return a
741  // pointer to it.
742  newSoln = curX_;
743  }
744  else // the update vector is NOT null.
745  {
746  if (updateLP) // The caller wants us to update the linear problem.
747  {
748  if (RP_.is_null())
749  { // There is no right preconditioner.
750  // curX_ := curX_ + scale * update.
751  MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ );
752  }
753  else
754  { // There is indeed a right preconditioner, so apply it
755  // before computing the new solution.
756  RCP<MV> rightPrecUpdate =
757  MVT::Clone (*update, MVT::GetNumberVecs (*update));
758  {
759 #ifdef BELOS_TEUCHOS_TIME_MONITOR
760  Teuchos::TimeMonitor PrecTimer (*timerPrec_);
761 #endif
762  OPT::Apply( *RP_, *update, *rightPrecUpdate );
763  }
764  // curX_ := curX_ + scale * rightPrecUpdate.
765  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ );
766  }
767  solutionUpdated_ = true;
768  newSoln = curX_;
769  }
770  else
771  { // Rather than updating our stored current solution curX_,
772  // we make a copy and compute the new solution in the
773  // copy, without modifying curX_.
774  newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
775  if (RP_.is_null())
776  { // There is no right preconditioner.
777  // newSoln := curX_ + scale * update.
778  MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln );
779  }
780  else
781  { // There is indeed a right preconditioner, so apply it
782  // before computing the new solution.
783  RCP<MV> rightPrecUpdate =
784  MVT::Clone (*update, MVT::GetNumberVecs (*update));
785  {
786 #ifdef BELOS_TEUCHOS_TIME_MONITOR
787  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
788 #endif
789  OPT::Apply( *RP_, *update, *rightPrecUpdate );
790  }
791  // newSoln := curX_ + scale * rightPrecUpdate.
792  MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln );
793  }
794  }
795  }
796  return newSoln;
797  }
798 
799  template <class ScalarType, class MV, class OP>
800  void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
801  {
802  if (label != label_) {
803  label_ = label;
804  // Create new timers if they have already been created.
805  if (timerOp_ != Teuchos::null) {
806  std::string opLabel = label_ + ": Operation Op*x";
807 #ifdef BELOS_TEUCHOS_TIME_MONITOR
808  timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
809 #endif
810  }
811  if (timerPrec_ != Teuchos::null) {
812  std::string precLabel = label_ + ": Operation Prec*x";
813 #ifdef BELOS_TEUCHOS_TIME_MONITOR
814  timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
815 #endif
816  }
817  }
818  }
819 
820  template <class ScalarType, class MV, class OP>
821  bool
823  setProblem (const Teuchos::RCP<MV> &newX,
824  const Teuchos::RCP<const MV> &newB)
825  {
826  // Create timers if the haven't been created yet.
827  if (timerOp_ == Teuchos::null) {
828  std::string opLabel = label_ + ": Operation Op*x";
829 #ifdef BELOS_TEUCHOS_TIME_MONITOR
830  timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
831 #endif
832  }
833  if (timerPrec_ == Teuchos::null) {
834  std::string precLabel = label_ + ": Operation Prec*x";
835 #ifdef BELOS_TEUCHOS_TIME_MONITOR
836  timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
837 #endif
838  }
839 
840  // Set the linear system using the arguments newX and newB
841  if (newX != Teuchos::null)
842  X_ = newX;
843  if (newB != Teuchos::null)
844  B_ = newB;
845 
846  // Invalidate the current linear system indices and multivectors.
847  rhsIndex_.resize(0);
848  curX_ = Teuchos::null;
849  curB_ = Teuchos::null;
850 
851  // If we didn't set a matrix A, a left-hand side X, or a
852  // right-hand side B, then we didn't set the problem.
853  if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
854  isSet_ = false;
855  return isSet_;
856  }
857 
858  // Reset whether the solution has been updated. (We're just
859  // setting the problem now, so of course we haven't updated the
860  // solution yet.)
861  solutionUpdated_ = false;
862 
863  // Compute the initial residuals.
864  if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
865  R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
866  }
867  computeCurrResVec( &*R0_, &*X_, &*B_ );
868 
869  if (LP_!=Teuchos::null) {
870  if (PR0_==Teuchos::null || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *B_ )) {
871  PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
872  }
873  {
874 #ifdef BELOS_TEUCHOS_TIME_MONITOR
875  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
876 #endif
877  OPT::Apply( *LP_, *R0_, *PR0_ );
878  }
879  }
880  else {
881  PR0_ = R0_;
882  }
883 
884  // The problem has been set and is ready for use.
885  isSet_ = true;
886 
887  // Return isSet.
888  return isSet_;
889  }
890 
891  template <class ScalarType, class MV, class OP>
893  {
894  if (isSet_) {
895  return curX_;
896  }
897  else {
898  return Teuchos::null;
899  }
900  }
901 
902  template <class ScalarType, class MV, class OP>
904  {
905  if (isSet_) {
906  return curB_;
907  }
908  else {
909  return Teuchos::null;
910  }
911  }
912 
913  template <class ScalarType, class MV, class OP>
914  void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
915  {
916  using Teuchos::null;
917  using Teuchos::RCP;
918 
919  const bool leftPrec = LP_ != null;
920  const bool rightPrec = RP_ != null;
921 
922  // We only need a temporary vector for intermediate results if
923  // there is a left or right preconditioner. We really should just
924  // keep this temporary vector around instead of allocating it each
925  // time.
926  RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
927 
928  //
929  // No preconditioning.
930  //
931  if (!leftPrec && !rightPrec){
932 #ifdef BELOS_TEUCHOS_TIME_MONITOR
933  Teuchos::TimeMonitor OpTimer(*timerOp_);
934 #endif
935  OPT::Apply( *A_, x, y );
936  }
937  //
938  // Preconditioning is being done on both sides
939  //
940  else if( leftPrec && rightPrec )
941  {
942  {
943 #ifdef BELOS_TEUCHOS_TIME_MONITOR
944  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
945 #endif
946  OPT::Apply( *RP_, x, y );
947  }
948  {
949 #ifdef BELOS_TEUCHOS_TIME_MONITOR
950  Teuchos::TimeMonitor OpTimer(*timerOp_);
951 #endif
952  OPT::Apply( *A_, y, *ytemp );
953  }
954  {
955 #ifdef BELOS_TEUCHOS_TIME_MONITOR
956  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
957 #endif
958  OPT::Apply( *LP_, *ytemp, y );
959  }
960  }
961  //
962  // Preconditioning is only being done on the left side
963  //
964  else if( leftPrec )
965  {
966  {
967 #ifdef BELOS_TEUCHOS_TIME_MONITOR
968  Teuchos::TimeMonitor OpTimer(*timerOp_);
969 #endif
970  OPT::Apply( *A_, x, *ytemp );
971  }
972  {
973 #ifdef BELOS_TEUCHOS_TIME_MONITOR
974  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
975 #endif
976  OPT::Apply( *LP_, *ytemp, y );
977  }
978  }
979  //
980  // Preconditioning is only being done on the right side
981  //
982  else
983  {
984  {
985 #ifdef BELOS_TEUCHOS_TIME_MONITOR
986  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
987 #endif
988  OPT::Apply( *RP_, x, *ytemp );
989  }
990  {
991 #ifdef BELOS_TEUCHOS_TIME_MONITOR
992  Teuchos::TimeMonitor OpTimer(*timerOp_);
993 #endif
994  OPT::Apply( *A_, *ytemp, y );
995  }
996  }
997  }
998 
999  template <class ScalarType, class MV, class OP>
1000  void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
1001  if (A_.get()) {
1002 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1003  Teuchos::TimeMonitor OpTimer(*timerOp_);
1004 #endif
1005  OPT::Apply( *A_,x, y);
1006  }
1007  else {
1008  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1009  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1010  }
1011  }
1012 
1013  template <class ScalarType, class MV, class OP>
1014  void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
1015  if (LP_!=Teuchos::null) {
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1017  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1018 #endif
1019  return ( OPT::Apply( *LP_,x, y) );
1020  }
1021  else {
1022  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1023  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1024  }
1025  }
1026 
1027  template <class ScalarType, class MV, class OP>
1028  void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
1029  if (RP_!=Teuchos::null) {
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1031  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1032 #endif
1033  return ( OPT::Apply( *RP_,x, y) );
1034  }
1035  else {
1036  MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x,
1037  Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
1038  }
1039  }
1040 
1041  template <class ScalarType, class MV, class OP>
1042  void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
1043 
1044  if (R) {
1045  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1046  {
1047  if (LP_!=Teuchos::null)
1048  {
1049  Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
1050  {
1051 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1052  Teuchos::TimeMonitor OpTimer(*timerOp_);
1053 #endif
1054  OPT::Apply( *A_, *X, *R_temp );
1055  }
1056  MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
1057  {
1058 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1059  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1060 #endif
1061  OPT::Apply( *LP_, *R_temp, *R );
1062  }
1063  }
1064  else
1065  {
1066  {
1067 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1068  Teuchos::TimeMonitor OpTimer(*timerOp_);
1069 #endif
1070  OPT::Apply( *A_, *X, *R );
1071  }
1072  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1073  }
1074  }
1075  else {
1076  // The solution and right-hand side may not be specified, check and use which ones exist.
1077  Teuchos::RCP<const MV> localB, localX;
1078  if (B)
1079  localB = Teuchos::rcp( B, false );
1080  else
1081  localB = curB_;
1082 
1083  if (X)
1084  localX = Teuchos::rcp( X, false );
1085  else
1086  localX = curX_;
1087 
1088  if (LP_!=Teuchos::null)
1089  {
1090  Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
1091  {
1092 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1093  Teuchos::TimeMonitor OpTimer(*timerOp_);
1094 #endif
1095  OPT::Apply( *A_, *localX, *R_temp );
1096  }
1097  MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
1098  {
1099 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1100  Teuchos::TimeMonitor PrecTimer(*timerPrec_);
1101 #endif
1102  OPT::Apply( *LP_, *R_temp, *R );
1103  }
1104  }
1105  else
1106  {
1107  {
1108 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1109  Teuchos::TimeMonitor OpTimer(*timerOp_);
1110 #endif
1111  OPT::Apply( *A_, *localX, *R );
1112  }
1113  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1114  }
1115  }
1116  }
1117  }
1118 
1119 
1120  template <class ScalarType, class MV, class OP>
1121  void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
1122 
1123  if (R) {
1124  if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
1125  {
1126  {
1127 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1128  Teuchos::TimeMonitor OpTimer(*timerOp_);
1129 #endif
1130  OPT::Apply( *A_, *X, *R );
1131  }
1132  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
1133  }
1134  else {
1135  // The solution and right-hand side may not be specified, check and use which ones exist.
1136  Teuchos::RCP<const MV> localB, localX;
1137  if (B)
1138  localB = Teuchos::rcp( B, false );
1139  else
1140  localB = curB_;
1141 
1142  if (X)
1143  localX = Teuchos::rcp( X, false );
1144  else
1145  localX = curX_;
1146 
1147  {
1148 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1149  Teuchos::TimeMonitor OpTimer(*timerOp_);
1150 #endif
1151  OPT::Apply( *A_, *localX, *R );
1152  }
1153  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
1154  }
1155  }
1156  }
1157 
1158 } // end Belos namespace
1159 
1160 #endif /* BELOS_LINEAR_PROBLEM_HPP */
1161 
1162 
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.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Teuchos::RCP< MV > getCurrLHSVec()
Get a pointer to the current left-hand side (solution) of the linear system.
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 Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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.
LinearProblem(void)
Default constructor.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
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.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
bool isSolutionUpdated() const
Has the current approximate solution been updated?
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
void setLeftPrec(const Teuchos::RCP< const OP > &LP)
Set left preconditioner (LP) of linear problem .
Teuchos::RCP< const MV > getInitResVec() const
A pointer to the initial unpreconditioned residual vector.
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.
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.
Teuchos::RCP< const MV > getInitPrecResVec() const
A pointer to the preconditioned initial residual vector.
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 setLHS(const Teuchos::RCP< MV > &X)
Set left-hand-side X of linear problem .

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11