Belos Version of the Day
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 
00049 #include "BelosMultiVecTraits.hpp"
00050 #include "BelosOperatorTraits.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00054 namespace Belos {
00055 
00057 
00058   
00061   class LinearProblemError : public BelosError {
00062   public: 
00063     LinearProblemError (const std::string& what_arg) : 
00064       BelosError(what_arg) {}
00065   };
00066   
00068 
00081   template <class ScalarType, class MV, class OP>
00082   class LinearProblem {
00083   public:
00084     
00086 
00087 
00094     LinearProblem (void);
00095     
00103     LinearProblem (const Teuchos::RCP<const OP> &A, 
00104        const Teuchos::RCP<MV> &X, 
00105        const Teuchos::RCP<const MV> &B);
00106     
00110     LinearProblem (const LinearProblem<ScalarType,MV,OP>& Problem);
00111     
00113     virtual ~LinearProblem (void);
00114 
00116     
00118 
00119     
00123     void setOperator (const Teuchos::RCP<const OP> &A) { 
00124       A_ = A; 
00125       isSet_=false; 
00126     }
00127     
00133     void setLHS (const Teuchos::RCP<MV> &X) { 
00134       X_ = X; 
00135       isSet_=false; 
00136     }
00137     
00142     void setRHS (const Teuchos::RCP<const MV> &B) { 
00143       B_ = B; 
00144       isSet_=false; 
00145     }
00146     
00150     void setLeftPrec(const Teuchos::RCP<const OP> &LP) {  LP_ = LP; }
00151 
00155     void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }
00156 
00165     void setCurrLS ();
00166 
00176     void setLSIndex (const std::vector<int>& index); 
00177     
00184     void setHermitian() { isHermitian_ = true; }
00185    
00192     void setLabel (const std::string& label) { label_ = label; }
00193 
00230     Teuchos::RCP<MV> 
00231     updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
00232         bool updateLP = false,
00233         ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());    
00234 
00251     Teuchos::RCP<MV> updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
00252                                     ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
00253     { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }
00254 
00256     
00258 
00259     
00282     bool 
00283     setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null, 
00284     const Teuchos::RCP<const MV> &newB = Teuchos::null);
00285 
00287     
00289 
00290     
00292     Teuchos::RCP<const OP> getOperator() const { return(A_); }
00293     
00295     Teuchos::RCP<MV> getLHS() const { return(X_); }
00296     
00298     Teuchos::RCP<const MV> getRHS() const { return(B_); }
00299     
00301     Teuchos::RCP<const MV> getInitResVec() const { return(R0_); }
00302     
00307     Teuchos::RCP<const MV> getInitPrecResVec() const { return(PR0_); }
00308     
00310 
00319     Teuchos::RCP<MV> getCurrLHSVec();
00320     
00322 
00331     Teuchos::RCP<const MV> getCurrRHSVec();
00332     
00334     Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
00335     
00337     Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
00338     
00357     const std::vector<int> getLSIndex() const { return(rhsIndex_); }
00358 
00364     int getLSNumber() const { return(lsNum_); }
00365 
00372     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00373       return Teuchos::tuple(timerOp_,timerPrec_);
00374     }
00375 
00376 
00378     
00380 
00381     
00389     bool isSolutionUpdated() const { return(solutionUpdated_); }
00390 
00392     bool isProblemSet() const { return(isSet_); }
00393     
00396     bool isHermitian() const { return(isHermitian_); }
00397     
00399     bool isLeftPrec() const { return(LP_!=Teuchos::null); }
00400 
00402     bool isRightPrec() const { return(RP_!=Teuchos::null); }
00403  
00405     
00407 
00408     
00410 
00417     void apply( const MV& x, MV& y ) const;
00418     
00420 
00427     void applyOp( const MV& x, MV& y ) const;
00428     
00430 
00434     void applyLeftPrec( const MV& x, MV& y ) const;
00435     
00437 
00441     void applyRightPrec( const MV& x, MV& y ) const;
00442     
00444 
00448     void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;
00449 
00451 
00455     void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
00456     
00458     
00459   private:
00460     
00462     Teuchos::RCP<const OP> A_;
00463     
00465     Teuchos::RCP<MV> X_;
00466     
00468     Teuchos::RCP<MV> curX_;
00469     
00471     Teuchos::RCP<const MV> B_;
00472     
00474     Teuchos::RCP<const MV> curB_;
00475     
00477     Teuchos::RCP<MV> R0_;
00478    
00480     Teuchos::RCP<MV> PR0_;
00481  
00483     Teuchos::RCP<const OP> LP_;  
00484     
00486     Teuchos::RCP<const OP> RP_;
00487     
00489     mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;
00490 
00492     int blocksize_;
00493 
00495     int num2Solve_;
00496     
00498     std::vector<int> rhsIndex_;    
00499 
00501     int lsNum_;
00502 
00504 
00505 
00507     bool Left_Scale_;
00508 
00510     bool Right_Scale_;
00511 
00513     bool isSet_;
00514 
00517     bool isHermitian_;
00518 
00520     bool solutionUpdated_;    
00521 
00523    
00525     std::string label_;
00526  
00527     typedef MultiVecTraits<ScalarType,MV>  MVT;
00528     typedef OperatorTraits<ScalarType,MV,OP>  OPT;
00529   };
00530   
00531   //--------------------------------------------
00532   //  Constructor Implementations
00533   //--------------------------------------------
00534   
00535   template <class ScalarType, class MV, class OP>
00536   LinearProblem<ScalarType,MV,OP>::LinearProblem(void) : 
00537     blocksize_(0),
00538     num2Solve_(0),
00539     rhsIndex_(0),
00540     lsNum_(0),
00541     Left_Scale_(false),
00542     Right_Scale_(false),
00543     isSet_(false),
00544     isHermitian_(false),
00545     solutionUpdated_(false),
00546     label_("Belos")
00547   {
00548   }
00549   
00550   template <class ScalarType, class MV, class OP>
00551   LinearProblem<ScalarType,MV,OP>::LinearProblem(const Teuchos::RCP<const OP> &A, 
00552              const Teuchos::RCP<MV> &X, 
00553              const Teuchos::RCP<const MV> &B
00554              ) :
00555     A_(A),
00556     X_(X),
00557     B_(B),
00558     blocksize_(0),
00559     num2Solve_(0),
00560     rhsIndex_(0),
00561     lsNum_(0),
00562     Left_Scale_(false),
00563     Right_Scale_(false),
00564     isSet_(false),
00565     isHermitian_(false),
00566     solutionUpdated_(false),
00567     label_("Belos")
00568   {
00569   }
00570   
00571   template <class ScalarType, class MV, class OP>
00572   LinearProblem<ScalarType,MV,OP>::LinearProblem(const LinearProblem<ScalarType,MV,OP>& Problem) :
00573     A_(Problem.A_),
00574     X_(Problem.X_),
00575     curX_(Problem.curX_),
00576     B_(Problem.B_),
00577     curB_(Problem.curB_),
00578     R0_(Problem.R0_),
00579     PR0_(Problem.PR0_),
00580     LP_(Problem.LP_),
00581     RP_(Problem.RP_),
00582     timerOp_(Problem.timerOp_),
00583     timerPrec_(Problem.timerPrec_),
00584     blocksize_(Problem.blocksize_),
00585     num2Solve_(Problem.num2Solve_),
00586     rhsIndex_(Problem.rhsIndex_),
00587     lsNum_(Problem.lsNum_),
00588     Left_Scale_(Problem.Left_Scale_),
00589     Right_Scale_(Problem.Right_Scale_),
00590     isSet_(Problem.isSet_),
00591     isHermitian_(Problem.isHermitian_),
00592     solutionUpdated_(Problem.solutionUpdated_),
00593     label_(Problem.label_)
00594   {
00595   }
00596   
00597   template <class ScalarType, class MV, class OP>
00598   LinearProblem<ScalarType,MV,OP>::~LinearProblem(void)
00599   {}
00600   
00601   template <class ScalarType, class MV, class OP>
00602   void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
00603   {
00604     // Set new linear systems using the indices in index.
00605     rhsIndex_ = index;
00606     
00607     // Compute the new block linear system.
00608     // ( first clean up old linear system )
00609     curB_ = Teuchos::null;
00610     curX_ = Teuchos::null;
00611    
00612     // Create indices for the new linear system.
00613     int validIdx = 0, ivalidIdx = 0;
00614     blocksize_ = rhsIndex_.size();
00615     std::vector<int> vldIndex( blocksize_ );
00616     std::vector<int> newIndex( blocksize_ );
00617     std::vector<int> iIndex( blocksize_ );
00618     for (int i=0; i<blocksize_; ++i) {
00619       if (rhsIndex_[i] > -1) {
00620         vldIndex[validIdx] = rhsIndex_[i];
00621         newIndex[validIdx] = i;
00622         validIdx++;
00623       }
00624       else {
00625         iIndex[ivalidIdx] = i;
00626         ivalidIdx++;
00627       }
00628     }
00629     vldIndex.resize(validIdx);
00630     newIndex.resize(validIdx);   
00631     iIndex.resize(ivalidIdx);
00632     num2Solve_ = validIdx;
00633 
00634     // Create the new linear system using index
00635     if (num2Solve_ != blocksize_) {
00636       newIndex.resize(num2Solve_);
00637       vldIndex.resize(num2Solve_);
00638       //
00639       // First create multivectors of blocksize.
00640       // Fill the RHS with random vectors LHS with zero vectors.
00641       curX_ = MVT::Clone( *X_, blocksize_ );
00642       MVT::MvInit(*curX_);
00643       Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
00644       MVT::MvRandom(*tmpCurB);
00645       //
00646       // Now put in the part of B into curB 
00647       Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
00648       MVT::SetBlock( *tptr, newIndex, *tmpCurB );
00649       curB_ = tmpCurB;
00650       //
00651       // Now put in the part of X into curX
00652       tptr = MVT::CloneView( *X_, vldIndex );
00653       MVT::SetBlock( *tptr, newIndex, *curX_ );
00654       //
00655       solutionUpdated_ = false;
00656     }
00657     else {
00658       curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
00659       curB_ = MVT::CloneView( *B_, rhsIndex_ );
00660     }
00661     //
00662     // Increment the number of linear systems that have been loaded into this object.
00663     //
00664     lsNum_++;
00665   }
00666 
00667 
00668   template <class ScalarType, class MV, class OP>
00669   void LinearProblem<ScalarType,MV,OP>::setCurrLS() 
00670   { 
00671     //
00672     // We only need to copy the solutions back if the linear systems of
00673     // interest are less than the block size.
00674     //
00675     if (num2Solve_ < blocksize_) {
00676       //
00677       // Get a view of the current solutions and correction std::vector.
00678       //
00679       int validIdx = 0;
00680       std::vector<int> newIndex( num2Solve_ );
00681       std::vector<int> vldIndex( num2Solve_ );
00682       for (int i=0; i<blocksize_; ++i) {
00683         if ( rhsIndex_[i] > -1 ) { 
00684           vldIndex[validIdx] = rhsIndex_[i];
00685     newIndex[validIdx] = i;
00686           validIdx++;
00687         } 
00688       }
00689       Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
00690       MVT::SetBlock( *tptr, vldIndex, *X_ );
00691     }
00692     //
00693     // Clear the current vectors of this linear system so that any future calls
00694     // to get the vectors for this system return null pointers.
00695     //
00696     curX_ = Teuchos::null;
00697     curB_ = Teuchos::null;
00698     rhsIndex_.resize(0);
00699   }
00700   
00701 
00702   template <class ScalarType, class MV, class OP>
00703   Teuchos::RCP<MV> 
00704   LinearProblem<ScalarType,MV,OP>::
00705   updateSolution (const Teuchos::RCP<MV>& update, 
00706       bool updateLP,
00707       ScalarType scale)
00708   { 
00709     using Teuchos::RCP;
00710     using Teuchos::null;
00711 
00712     RCP<MV> newSoln;
00713     if (update.is_null())
00714       { // The caller didn't supply an update vector, so we assume
00715   // that the current solution curX_ is unchanged, and return a
00716   // pointer to it.
00717   newSoln = curX_;
00718       }
00719     else // the update vector is NOT null.
00720       { 
00721   if (updateLP) // The caller wants us to update the linear problem.
00722     { 
00723       if (RP_.is_null()) 
00724         { // There is no right preconditioner.
00725     // curX_ := curX_ + scale * update.
00726     MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ ); 
00727         }
00728       else 
00729         { // There is indeed a right preconditioner, so apply it
00730     // before computing the new solution.
00731     RCP<MV> rightPrecUpdate = 
00732       MVT::Clone (*update, MVT::GetNumberVecs (*update));
00733     {
00734 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00735       Teuchos::TimeMonitor PrecTimer (*timerPrec_);
00736 #endif
00737       OPT::Apply( *RP_, *update, *rightPrecUpdate ); 
00738     }
00739     // curX_ := curX_ + scale * rightPrecUpdate.
00740     MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ ); 
00741         } 
00742       solutionUpdated_ = true; 
00743       newSoln = curX_;
00744     }
00745   else 
00746     { // Rather than updating our stored current solution curX_,
00747       // we make a copy and compute the new solution in the
00748       // copy, without modifying curX_.
00749       newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
00750       if (RP_.is_null())
00751         { // There is no right preconditioner.
00752     // newSoln := curX_ + scale * update.
00753     MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln ); 
00754         }
00755       else 
00756         { // There is indeed a right preconditioner, so apply it
00757     // before computing the new solution.
00758     RCP<MV> rightPrecUpdate = 
00759       MVT::Clone (*update, MVT::GetNumberVecs (*update));
00760     {
00761 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00762       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00763 #endif
00764       OPT::Apply( *RP_, *update, *rightPrecUpdate ); 
00765     }
00766     // newSoln := curX_ + scale * rightPrecUpdate.
00767     MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln ); 
00768         } 
00769     }
00770       }
00771     return newSoln;
00772   }
00773   
00774 
00775   template <class ScalarType, class MV, class OP>
00776   bool 
00777   LinearProblem<ScalarType,MV,OP>::
00778   setProblem (const Teuchos::RCP<MV> &newX, 
00779         const Teuchos::RCP<const MV> &newB)
00780   {
00781     // Create timers if the haven't been created yet.
00782     if (timerOp_ == Teuchos::null) {
00783       std::string opLabel = label_ + ": Operation Op*x";
00784 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00785       timerOp_ = Teuchos::TimeMonitor::getNewTimer( opLabel );
00786 #endif
00787     }
00788     if (timerPrec_ == Teuchos::null) {
00789       std::string precLabel = label_ + ": Operation Prec*x";
00790 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00791       timerPrec_ = Teuchos::TimeMonitor::getNewTimer( precLabel );
00792 #endif
00793     }
00794 
00795     // Set the linear system using the arguments newX and newB
00796     if (newX != Teuchos::null)
00797       X_ = newX;
00798     if (newB != Teuchos::null)
00799       B_ = newB;
00800 
00801     // Invalidate the current linear system indices and multivectors.
00802     rhsIndex_.resize(0);
00803     curX_ = Teuchos::null;
00804     curB_ = Teuchos::null;
00805 
00806     // If we didn't set a matrix A, a left-hand side X, or a
00807     // right-hand side B, then we didn't set the problem.
00808     if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
00809       isSet_ = false;
00810       return isSet_;
00811     }
00812 
00813     // Reset whether the solution has been updated.  (We're just
00814     // setting the problem now, so of course we haven't updated the
00815     // solution yet.)
00816     solutionUpdated_ = false;
00817     
00818     // Compute the initial residuals.
00819     if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
00820       R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
00821     }
00822     computeCurrResVec( &*R0_, &*X_, &*B_ );
00823 
00824     if (LP_!=Teuchos::null) {
00825       if (PR0_==Teuchos::null || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *B_ )) {
00826         PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
00827       }
00828       {
00829 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00830         Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00831 #endif
00832         OPT::Apply( *LP_, *R0_, *PR0_ );
00833       }
00834     } 
00835     else {
00836       PR0_ = R0_;
00837     }    
00838 
00839     // The problem has been set and is ready for use.
00840     isSet_ = true;
00841     
00842     // Return isSet.
00843     return isSet_;
00844   }
00845   
00846   template <class ScalarType, class MV, class OP>
00847   Teuchos::RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrLHSVec()
00848   {
00849     if (isSet_) {
00850       return curX_;
00851     }
00852     else {
00853       return Teuchos::null;
00854     }
00855   }
00856   
00857   template <class ScalarType, class MV, class OP>
00858   Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getCurrRHSVec()
00859   {
00860     if (isSet_) {
00861       return curB_;
00862     }
00863     else {
00864       return Teuchos::null;
00865     }
00866   }
00867   
00868   template <class ScalarType, class MV, class OP>
00869   void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
00870   {
00871     using Teuchos::null;
00872     using Teuchos::RCP;
00873 
00874     const bool leftPrec = LP_ != null;
00875     const bool rightPrec = RP_ != null;
00876 
00877     // We only need a temporary vector for intermediate results if
00878     // there is a left or right preconditioner.  We really should just
00879     // keep this temporary vector around instead of allocating it each
00880     // time.
00881     RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
00882 
00883     //
00884     // No preconditioning.
00885     // 
00886     if (!leftPrec && !rightPrec){ 
00887 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00888       Teuchos::TimeMonitor OpTimer(*timerOp_);
00889 #endif
00890       OPT::Apply( *A_, x, y );
00891     }
00892     //
00893     // Preconditioning is being done on both sides
00894     //
00895     else if( leftPrec && rightPrec ) 
00896       {
00897         {
00898 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00899           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00900 #endif
00901     OPT::Apply( *RP_, x, y );   
00902         }
00903         {
00904 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00905           Teuchos::TimeMonitor OpTimer(*timerOp_);
00906 #endif
00907     OPT::Apply( *A_, y, *ytemp );
00908         }
00909         {
00910 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00911           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00912 #endif
00913     OPT::Apply( *LP_, *ytemp, y );
00914         }
00915       }
00916     //
00917     // Preconditioning is only being done on the left side
00918     //
00919     else if( leftPrec ) 
00920       {
00921         {
00922 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00923           Teuchos::TimeMonitor OpTimer(*timerOp_);
00924 #endif
00925     OPT::Apply( *A_, x, *ytemp );
00926         }
00927         {
00928 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00929           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00930 #endif
00931     OPT::Apply( *LP_, *ytemp, y );
00932         }
00933       }
00934     //
00935     // Preconditioning is only being done on the right side
00936     //
00937     else 
00938       {
00939         {
00940 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00941           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00942 #endif
00943     OPT::Apply( *RP_, x, *ytemp );
00944         }
00945         {
00946 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00947           Teuchos::TimeMonitor OpTimer(*timerOp_);
00948 #endif
00949           OPT::Apply( *A_, *ytemp, y );
00950         }
00951       }  
00952   }
00953   
00954   template <class ScalarType, class MV, class OP>
00955   void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
00956     if (A_.get()) {
00957 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00958       Teuchos::TimeMonitor OpTimer(*timerOp_);
00959 #endif
00960       OPT::Apply( *A_,x, y);   
00961     }
00962     else {
00963       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00964         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00965     }
00966   }
00967   
00968   template <class ScalarType, class MV, class OP>
00969   void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
00970     if (LP_!=Teuchos::null) {
00971 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00972       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00973 #endif
00974       return ( OPT::Apply( *LP_,x, y) );
00975     }
00976     else {
00977       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00978         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00979     }
00980   }
00981   
00982   template <class ScalarType, class MV, class OP>
00983   void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
00984     if (RP_!=Teuchos::null) {
00985 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00986       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00987 #endif
00988       return ( OPT::Apply( *RP_,x, y) );
00989     }
00990     else {
00991       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
00992         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
00993     }
00994   }
00995   
00996   template <class ScalarType, class MV, class OP>
00997   void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
00998 
00999     if (R) {
01000       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
01001   {
01002     if (LP_!=Teuchos::null)
01003       {
01004         Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
01005               {
01006 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01007                 Teuchos::TimeMonitor OpTimer(*timerOp_);
01008 #endif
01009           OPT::Apply( *A_, *X, *R_temp );
01010               }
01011         MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
01012               {
01013 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01014                 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01015 #endif
01016           OPT::Apply( *LP_, *R_temp, *R );
01017               }
01018       }
01019     else 
01020       {
01021               {
01022 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01023                 Teuchos::TimeMonitor OpTimer(*timerOp_);
01024 #endif
01025           OPT::Apply( *A_, *X, *R );
01026               }
01027         MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
01028       }
01029   }
01030       else { 
01031   // The solution and right-hand side may not be specified, check and use which ones exist.
01032   Teuchos::RCP<const MV> localB, localX;
01033   if (B)
01034     localB = Teuchos::rcp( B, false );
01035   else
01036     localB = curB_;
01037   
01038   if (X)
01039     localX = Teuchos::rcp( X, false );
01040   else
01041     localX = curX_;
01042   
01043   if (LP_!=Teuchos::null)
01044     {
01045       Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
01046             {
01047 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01048               Teuchos::TimeMonitor OpTimer(*timerOp_);
01049 #endif
01050         OPT::Apply( *A_, *localX, *R_temp );
01051             }
01052       MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
01053             {
01054 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01055               Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01056 #endif
01057         OPT::Apply( *LP_, *R_temp, *R );
01058             }
01059     }
01060   else 
01061     {
01062             {
01063 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01064               Teuchos::TimeMonitor OpTimer(*timerOp_);
01065 #endif
01066           OPT::Apply( *A_, *localX, *R );
01067             }
01068       MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
01069     }
01070       }    
01071     } 
01072   }
01073   
01074   
01075   template <class ScalarType, class MV, class OP>
01076   void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
01077 
01078     if (R) {
01079       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
01080   {
01081           {
01082 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01083             Teuchos::TimeMonitor OpTimer(*timerOp_);
01084 #endif
01085       OPT::Apply( *A_, *X, *R );
01086           }
01087     MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
01088   }
01089       else { 
01090   // The solution and right-hand side may not be specified, check and use which ones exist.
01091   Teuchos::RCP<const MV> localB, localX;
01092   if (B)
01093     localB = Teuchos::rcp( B, false );
01094   else
01095     localB = curB_;
01096   
01097   if (X)
01098     localX = Teuchos::rcp( X, false );
01099   else
01100     localX = curX_;
01101     
01102         {
01103 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01104           Teuchos::TimeMonitor OpTimer(*timerOp_);
01105 #endif
01106     OPT::Apply( *A_, *localX, *R );
01107         }
01108   MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
01109       }    
01110     }
01111   }
01112   
01113 } // end Belos namespace
01114 
01115 #endif /* BELOS_LINEAR_PROBLEM_HPP */
01116 
01117 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines