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