Belos Package Browser (Single Doxygen Collection) Development
BelosLinearProblem.hpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //                 Belos: Block Linear Solvers Package
00006 //                  Copyright 2004 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #ifndef BELOS_LINEAR_PROBLEM_HPP
00044 #define BELOS_LINEAR_PROBLEM_HPP
00045 
00050 #include "BelosMultiVecTraits.hpp"
00051 #include "BelosOperatorTraits.hpp"
00052 #include "Teuchos_ParameterList.hpp"
00053 #include "Teuchos_TimeMonitor.hpp"
00054 
00055 using Teuchos::RCP;
00056 using Teuchos::rcp;
00057 using Teuchos::null;
00058 using Teuchos::rcp_const_cast;
00059 using Teuchos::ParameterList;
00060 
00066 namespace Belos {
00067 
00069 
00070   
00073   class LinearProblemError : public BelosError
00074   {public: LinearProblemError(const std::string& what_arg) : BelosError(what_arg) {}};
00075   
00077   
00078 
00079   template <class ScalarType, class MV, class OP>
00080   class LinearProblem {
00081    
00082   public:
00083     
00085 
00086 
00087 
00091     LinearProblem(void);
00092     
00094 
00099     LinearProblem(const RCP<const OP> &A, 
00100       const RCP<MV> &X, 
00101       const RCP<const MV> &B
00102       );
00103     
00105 
00107     LinearProblem(const LinearProblem<ScalarType,MV,OP>& Problem);
00108     
00110 
00112     virtual ~LinearProblem(void);
00114     
00116 
00117     
00119 
00121     void setOperator(const RCP<const OP> &A) { A_ = A; isSet_=false; }
00122     
00124 
00126     void setLHS(const RCP<MV> &X) { X_ = X; isSet_=false; }
00127     
00129 
00131     void setRHS(const RCP<const MV> &B) { B_ = B; isSet_=false; }
00132     
00134 
00136     void setLeftPrec(const RCP<const OP> &LP) {  LP_ = LP; }
00137     
00139 
00141     void setRightPrec(const RCP<const OP> &RP) { RP_ = RP; }
00142     
00144 
00149     void setCurrLS();
00150     
00152 
00158     void setLSIndex(std::vector<int>& index); 
00159     
00161 
00165     void setHermitian(){ isHermitian_ = true; }
00166    
00168 
00171     void setLabel(const std::string& label) { label_ = label; }
00172 
00174 
00179     RCP<MV> updateSolution( const RCP<MV>& update = null,
00180             bool updateLP = false,
00181                                     ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() );    
00182 
00184     RCP<MV> updateSolution( const RCP<MV>& update = null,
00185                                     ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
00186     { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
00187 
00189     
00191 
00192     
00194 
00199     bool setProblem( const RCP<MV> &newX = null, const RCP<const MV> &newB = null );
00200 
00202     
00204 
00205     
00207     RCP<const OP> getOperator() const { return(A_); }
00208     
00210     RCP<MV> getLHS() const { return(X_); }
00211     
00213     RCP<const MV> getRHS() const { return(B_); }
00214     
00216 
00218     RCP<const MV> getInitResVec() const { return(R0_); }
00219     
00221 
00223     RCP<const MV> getInitPrecResVec() const { return(PR0_); }
00224     
00226 
00233     RCP<MV> getCurrLHSVec();
00234     
00236 
00243     RCP<MV> getCurrRHSVec();
00244     
00246     RCP<const OP> getLeftPrec() const { return(LP_); };
00247     
00249     RCP<const OP> getRightPrec() const { return(RP_); };
00250     
00252 
00263     const std::vector<int> getLSIndex() const { return(rhsIndex_); }
00264 
00266     /* This can be used by status test classes to determine if the solver manager has advanced 
00267        and is solving another linear system.
00268     */
00269     int getLSNumber() const { return(lsNum_); }
00270 
00277     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00278       return tuple(timerOp_,timerPrec_);
00279     }
00280 
00281 
00283     
00285 
00286     
00288 
00292     bool isSolutionUpdated() const { return(solutionUpdated_); }
00293 
00295     bool isProblemSet() const { return(isSet_); }
00296     
00298     bool isHermitian() const { return(isHermitian_); }
00299     
00301     bool isLeftPrec() const { return(LP_!=null); }
00302 
00304     bool isRightPrec() const { return(RP_!=null); }
00305  
00307     
00309 
00310     
00312 
00319     void apply( const MV& x, MV& y ) const;
00320     
00322 
00329     void applyOp( const MV& x, MV& y ) const;
00330     
00332 
00336     void applyLeftPrec( const MV& x, MV& y ) const;
00337     
00339 
00343     void applyRightPrec( const MV& x, MV& y ) const;
00344     
00346 
00350     void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
00351 
00353 
00357     void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
00358     
00360     
00361   private:
00362     
00364     RCP<const OP> A_;
00365     
00367     RCP<MV> X_;
00368     
00370     RCP<MV> curX_;
00371     
00373     RCP<const MV> B_;
00374     
00376     RCP<MV> curB_;
00377     
00379     RCP<MV> R0_;
00380    
00382     RCP<MV> PR0_;
00383  
00385     RCP<const OP> LP_;  
00386     
00388     RCP<const OP> RP_;
00389     
00391     mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
00392 
00394     int blocksize_;
00395 
00397     int num2Solve_;
00398     
00400     std::vector<int> rhsIndex_;    
00401 
00403     int lsNum_;
00404 
00406     bool Left_Scale_;
00407     bool Right_Scale_;
00408     bool isSet_;
00409     bool isHermitian_;
00410     bool solutionUpdated_;    
00411    
00413     std::string label_;
00414  
00415     typedef MultiVecTraits<ScalarType,MV>  MVT;
00416     typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00417   };
00418   
00419   //--------------------------------------------
00420   //  Constructor Implementations
00421   //--------------------------------------------
00422   
00423   template <class ScalarType, class MV, class OP>
00424   LinearProblem<ScalarType,MV,OP>::LinearProblem(void) : 
00425     blocksize_(0),
00426     num2Solve_(0),
00427     rhsIndex_(0),
00428     lsNum_(0),
00429     Left_Scale_(false),
00430     Right_Scale_(false),
00431     isSet_(false),
00432     isHermitian_(false),
00433     solutionUpdated_(false),
00434     label_("Belos")
00435   {
00436   }
00437   
00438   template <class ScalarType, class MV, class OP>
00439   LinearProblem<ScalarType,MV,OP>::LinearProblem(const RCP<const OP> &A, 
00440              const RCP<MV> &X, 
00441              const RCP<const MV> &B
00442              ) :
00443     A_(A),
00444     X_(X),
00445     B_(B),
00446     blocksize_(0),
00447     num2Solve_(0),
00448     rhsIndex_(0),
00449     lsNum_(0),
00450     Left_Scale_(false),
00451     Right_Scale_(false),
00452     isSet_(false),
00453     isHermitian_(false),
00454     solutionUpdated_(false),
00455     label_("Belos")
00456   {
00457   }
00458   
00459   template <class ScalarType, class MV, class OP>
00460   LinearProblem<ScalarType,MV,OP>::LinearProblem(const LinearProblem<ScalarType,MV,OP>& Problem) :
00461     A_(Problem.A_),
00462     X_(Problem.X_),
00463     curX_(Problem.curX_),
00464     B_(Problem.B_),
00465     curB_(Problem.curB_),
00466     R0_(Problem.R0_),
00467     PR0_(Problem.PR0_),
00468     LP_(Problem.LP_),
00469     RP_(Problem.RP_),
00470     timerOp_(Problem.timerOp_),
00471     timerPrec_(Problem.timerPrec_),
00472     blocksize_(Problem.blocksize_),
00473     num2Solve_(Problem.num2Solve_),
00474     rhsIndex_(Problem.rhsIndex_),
00475     lsNum_(Problem.lsNum_),
00476     Left_Scale_(Problem.Left_Scale_),
00477     Right_Scale_(Problem.Right_Scale_),
00478     isSet_(Problem.isSet_),
00479     isHermitian_(Problem.isHermitian_),
00480     solutionUpdated_(Problem.solutionUpdated_),
00481     label_(Problem.label_)
00482   {
00483   }
00484   
00485   template <class ScalarType, class MV, class OP>
00486   LinearProblem<ScalarType,MV,OP>::~LinearProblem(void)
00487   {}
00488   
00489   template <class ScalarType, class MV, class OP>
00490   void LinearProblem<ScalarType,MV,OP>::setLSIndex(std::vector<int>& index)
00491   {
00492     // Set new linear systems using the indices in index.
00493     rhsIndex_ = index;
00494     
00495     // Compute the new block linear system.
00496     // ( first clean up old linear system )
00497     curB_ = null;
00498     curX_ = null;
00499    
00500     // Create indices for the new linear system.
00501     int validIdx = 0, ivalidIdx = 0;
00502     blocksize_ = rhsIndex_.size();
00503     std::vector<int> vldIndex( blocksize_ );
00504     std::vector<int> newIndex( blocksize_ );
00505     std::vector<int> iIndex( blocksize_ );
00506     for (int i=0; i<blocksize_; ++i) {
00507       if (rhsIndex_[i] > -1) {
00508         vldIndex[validIdx] = rhsIndex_[i];
00509         newIndex[validIdx] = i;
00510         validIdx++;
00511       }
00512       else {
00513         iIndex[ivalidIdx] = i;
00514         ivalidIdx++;
00515       }
00516     }
00517     vldIndex.resize(validIdx);
00518     newIndex.resize(validIdx);   
00519     iIndex.resize(ivalidIdx);
00520     num2Solve_ = validIdx;
00521 
00522     // Create the new linear system using index
00523     if (num2Solve_ != blocksize_) {
00524       newIndex.resize(num2Solve_);
00525       vldIndex.resize(num2Solve_);
00526       //
00527       // First create multivectors of blocksize.
00528       // Fill the RHS with random vectors LHS with zero vectors.
00529       curX_ = MVT::Clone( *X_, blocksize_ );
00530       MVT::MvInit(*curX_);
00531       curB_ = MVT::Clone( *B_, blocksize_ );
00532       MVT::MvRandom(*curB_);
00533       //
00534       // Now put in the part of B into curB 
00535       RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
00536       MVT::SetBlock( *tptr, newIndex, *curB_ );
00537       //
00538       // Now put in the part of X into curX
00539       tptr = MVT::CloneView( *X_, vldIndex );
00540       MVT::SetBlock( *tptr, newIndex, *curX_ );
00541       //
00542       solutionUpdated_ = false;
00543     }
00544     else {
00545       curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
00546       curB_ = rcp_const_cast<MV>(MVT::CloneView( *B_, rhsIndex_ ));
00547     }
00548     //
00549     // Increment the number of linear systems that have been loaded into this object.
00550     //
00551     lsNum_++;
00552   }
00553 
00554 
00555   template <class ScalarType, class MV, class OP>
00556   void LinearProblem<ScalarType,MV,OP>::setCurrLS() 
00557   { 
00558     //
00559     // We only need to copy the solutions back if the linear systems of
00560     // interest are less than the block size.
00561     //
00562     if (num2Solve_ < blocksize_) {
00563       //
00564       // Get a view of the current solutions and correction std::vector.
00565       //
00566       int validIdx = 0;
00567       std::vector<int> newIndex( num2Solve_ );
00568       std::vector<int> vldIndex( num2Solve_ );
00569       for (int i=0; i<blocksize_; ++i) {
00570         if ( rhsIndex_[i] > -1 ) { 
00571           vldIndex[validIdx] = rhsIndex_[i];
00572     newIndex[validIdx] = i;
00573           validIdx++;
00574         } 
00575       }
00576       RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
00577       MVT::SetBlock( *tptr, vldIndex, *X_ );
00578     }
00579     //
00580     // Clear the current vectors of this linear system so that any future calls
00581     // to get the vectors for this system return null pointers.
00582     //
00583     curX_ = null;
00584     curB_ = null;
00585     rhsIndex_.resize(0);
00586   }
00587   
00588 
00589   template <class ScalarType, class MV, class OP>
00590   RCP<MV> LinearProblem<ScalarType,MV,OP>::updateSolution( const RCP<MV>& update, 
00591                  bool updateLP,
00592                  ScalarType scale )
00593   { 
00594     RCP<MV> newSoln;
00595     if (update != null) {
00596       if (updateLP == true) {
00597   if (RP_!=null) {
00598     //
00599     // Apply the right preconditioner before computing the current solution.
00600     RCP<MV> TrueUpdate = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
00601     {
00602       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00603       OPT::Apply( *RP_, *update, *TrueUpdate ); 
00604     }
00605     MVT::MvAddMv( 1.0, *curX_, scale, *TrueUpdate, *curX_ ); 
00606   } 
00607   else {
00608     MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ ); 
00609   }
00610   solutionUpdated_ = true; 
00611   newSoln = curX_;
00612       }
00613       else {
00614   newSoln = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
00615   if (RP_!=null) {
00616     //
00617     // Apply the right preconditioner before computing the current solution.
00618     RCP<MV> trueUpdate = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
00619     {
00620       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00621       OPT::Apply( *RP_, *update, *trueUpdate ); 
00622     }
00623     MVT::MvAddMv( 1.0, *curX_, scale, *trueUpdate, *newSoln ); 
00624   } 
00625   else {
00626     MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln ); 
00627   }
00628       }
00629     }
00630     else {
00631       newSoln = curX_;
00632     }
00633     return newSoln;
00634   }
00635   
00636 
00637   template <class ScalarType, class MV, class OP>
00638   bool LinearProblem<ScalarType,MV,OP>::setProblem( const RCP<MV> &newX, const RCP<const MV> &newB )
00639   {
00640     // Create timers if the haven't been created yet.
00641     if (timerOp_ == Teuchos::null) {
00642       std::string opLabel = label_ + ": Operation Op*x";
00643       timerOp_ = Teuchos::TimeMonitor::getNewTimer( opLabel );
00644     }
00645     if (timerPrec_ == Teuchos::null) {
00646       std::string precLabel = label_ + ": Operation Prec*x";
00647       timerPrec_ = Teuchos::TimeMonitor::getNewTimer( precLabel );
00648     }
00649 
00650     // Set the linear system using the arguments newX and newB
00651     if (newX != null)
00652       X_ = newX;
00653     if (newB != null)
00654       B_ = newB;
00655 
00656     // Invalidate the current linear system indices and multivectors.
00657     rhsIndex_.resize(0);
00658     curX_ = null;
00659     curB_ = null;
00660 
00661     // Check the validity of the linear problem object.
00662     // If no operator A exists, then throw an std::exception.
00663     if (A_ == null || X_ == null || B_ == null) {
00664       isSet_ = false;
00665       return isSet_;
00666     }
00667 
00668     // Initialize the state booleans
00669     solutionUpdated_ = false;
00670     
00671     // Compute the initial residuals.
00672     if (R0_==null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *X_ )) {
00673       R0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
00674     }
00675     computeCurrResVec( &*R0_, &*X_, &*B_ );
00676 
00677     if (LP_!=null) {
00678       if (PR0_==null || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *X_ )) {
00679         PR0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
00680       }
00681       {
00682         Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00683         OPT::Apply( *LP_, *R0_, *PR0_ );
00684       }
00685     } 
00686     else {
00687       PR0_ = R0_;
00688     }    
00689 
00690     // The problem has been set and is ready for use.
00691     isSet_ = true;
00692     
00693     // Return isSet.
00694     return isSet_;
00695   }
00696   
00697   template <class ScalarType, class MV, class OP>
00698   RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrLHSVec()
00699   {
00700     if (isSet_) {
00701       return curX_;
00702     }
00703     else {
00704       return Teuchos::null;
00705     }
00706   }
00707   
00708   template <class ScalarType, class MV, class OP>
00709   RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrRHSVec()
00710   {
00711     if (isSet_) {
00712       return curB_;
00713     }
00714     else {
00715       return Teuchos::null;
00716     }
00717   }
00718   
00719   template <class ScalarType, class MV, class OP>
00720   void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
00721   {
00722     RCP<MV> ytemp = MVT::Clone( y, MVT::GetNumberVecs( y ) );
00723     bool leftPrec = LP_!=null;
00724     bool rightPrec = RP_!=null;
00725     //
00726     // No preconditioning.
00727     // 
00728     if (!leftPrec && !rightPrec){ 
00729       Teuchos::TimeMonitor OpTimer(*timerOp_);
00730       OPT::Apply( *A_, x, y );
00731     }
00732     //
00733     // Preconditioning is being done on both sides
00734     //
00735     else if( leftPrec && rightPrec ) 
00736       {
00737         {
00738           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00739     OPT::Apply( *RP_, x, y );   
00740         }
00741         {
00742           Teuchos::TimeMonitor OpTimer(*timerOp_);
00743     OPT::Apply( *A_, y, *ytemp );
00744         }
00745         {
00746           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00747     OPT::Apply( *LP_, *ytemp, y );
00748         }
00749       }
00750     //
00751     // Preconditioning is only being done on the left side
00752     //
00753     else if( leftPrec ) 
00754       {
00755         {
00756           Teuchos::TimeMonitor OpTimer(*timerOp_);
00757     OPT::Apply( *A_, x, *ytemp );
00758         }
00759         {
00760           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00761     OPT::Apply( *LP_, *ytemp, y );
00762         }
00763       }
00764     //
00765     // Preconditioning is only being done on the right side
00766     //
00767     else 
00768       {
00769         {
00770           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00771     OPT::Apply( *RP_, x, *ytemp );
00772         }
00773         {
00774           Teuchos::TimeMonitor OpTimer(*timerOp_);
00775           OPT::Apply( *A_, *ytemp, y );
00776         }
00777       }  
00778   }
00779   
00780   template <class ScalarType, class MV, class OP>
00781   void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
00782     if (A_.get()) {
00783       Teuchos::TimeMonitor OpTimer(*timerOp_);
00784       OPT::Apply( *A_,x, y);   
00785     }
00786     else {
00787       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00788         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00789     }
00790   }
00791   
00792   template <class ScalarType, class MV, class OP>
00793   void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
00794     if (LP_!=null) {
00795       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00796       return ( OPT::Apply( *LP_,x, y) );
00797     }
00798     else {
00799       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00800         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00801     }
00802   }
00803   
00804   template <class ScalarType, class MV, class OP>
00805   void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
00806     if (RP_!=null) {
00807       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00808       return ( OPT::Apply( *RP_,x, y) );
00809     }
00810     else {
00811       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00812         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00813     }
00814   }
00815   
00816   template <class ScalarType, class MV, class OP>
00817   void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
00818 
00819     if (R) {
00820       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
00821   {
00822     if (LP_!=null)
00823       {
00824         RCP<MV> R_temp = MVT::Clone( *X, MVT::GetNumberVecs( *X ) );
00825               {
00826                 Teuchos::TimeMonitor OpTimer(*timerOp_);
00827           OPT::Apply( *A_, *X, *R_temp );
00828               }
00829         MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
00830               {
00831                 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00832           OPT::Apply( *LP_, *R_temp, *R );
00833               }
00834       }
00835     else 
00836       {
00837               {
00838                 Teuchos::TimeMonitor OpTimer(*timerOp_);
00839           OPT::Apply( *A_, *X, *R );
00840               }
00841         MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
00842       }
00843   }
00844       else { 
00845   // The solution and right-hand side may not be specified, check and use which ones exist.
00846   RCP<const MV> localB, localX;
00847   if (B)
00848     localB = rcp( B, false );
00849   else
00850     localB = curB_;
00851   
00852   if (X)
00853     localX = rcp( X, false );
00854   else
00855     localX = curX_;
00856   
00857   if (LP_!=null)
00858     {
00859       RCP<MV> R_temp = MVT::Clone( *localX, MVT::GetNumberVecs( *localX ) );
00860             {
00861               Teuchos::TimeMonitor OpTimer(*timerOp_);
00862         OPT::Apply( *A_, *localX, *R_temp );
00863             }
00864       MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
00865             {
00866               Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00867         OPT::Apply( *LP_, *R_temp, *R );
00868             }
00869     }
00870   else 
00871     {
00872             {
00873               Teuchos::TimeMonitor OpTimer(*timerOp_);
00874           OPT::Apply( *A_, *localX, *R );
00875             }
00876       MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
00877     }
00878       }    
00879     } 
00880   }
00881   
00882   
00883   template <class ScalarType, class MV, class OP>
00884   void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
00885 
00886     if (R) {
00887       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
00888   {
00889           {
00890             Teuchos::TimeMonitor OpTimer(*timerOp_);
00891       OPT::Apply( *A_, *X, *R );
00892           }
00893     MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
00894   }
00895       else { 
00896   // The solution and right-hand side may not be specified, check and use which ones exist.
00897   RCP<const MV> localB, localX;
00898   if (B)
00899     localB = rcp( B, false );
00900   else
00901     localB = curB_;
00902   
00903   if (X)
00904     localX = rcp( X, false );
00905   else
00906     localX = curX_;
00907     
00908         {
00909           Teuchos::TimeMonitor OpTimer(*timerOp_);
00910     OPT::Apply( *A_, *localX, *R );
00911         }
00912   MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
00913       }    
00914     }
00915   }
00916   
00917 } // end Belos namespace
00918 
00919 #endif /* BELOS_LINEAR_PROBLEM_HPP */
00920 
00921 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines