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     
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( bool isSym = true ) { isHermitian_ = isSym; }
00190    
00197     void setLabel (const std::string& 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   template <class ScalarType, class MV, class OP>
00800   void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
00801   {
00802     if (label != label_) {
00803       label_ = label;
00804       // Create new timers if they have already been created.
00805       if (timerOp_ != Teuchos::null) {
00806         std::string opLabel = label_ + ": Operation Op*x";
00807 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00808         timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
00809 #endif
00810       }
00811       if (timerPrec_ != Teuchos::null) {
00812         std::string precLabel = label_ + ": Operation Prec*x";
00813 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00814         timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
00815 #endif
00816       }
00817     }
00818   }
00819 
00820   template <class ScalarType, class MV, class OP>
00821   bool 
00822   LinearProblem<ScalarType,MV,OP>::
00823   setProblem (const Teuchos::RCP<MV> &newX, 
00824         const Teuchos::RCP<const MV> &newB)
00825   {
00826     // Create timers if the haven't been created yet.
00827     if (timerOp_ == Teuchos::null) {
00828       std::string opLabel = label_ + ": Operation Op*x";
00829 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00830       timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
00831 #endif
00832     }
00833     if (timerPrec_ == Teuchos::null) {
00834       std::string precLabel = label_ + ": Operation Prec*x";
00835 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00836       timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
00837 #endif
00838     }
00839 
00840     // Set the linear system using the arguments newX and newB
00841     if (newX != Teuchos::null)
00842       X_ = newX;
00843     if (newB != Teuchos::null)
00844       B_ = newB;
00845 
00846     // Invalidate the current linear system indices and multivectors.
00847     rhsIndex_.resize(0);
00848     curX_ = Teuchos::null;
00849     curB_ = Teuchos::null;
00850 
00851     // If we didn't set a matrix A, a left-hand side X, or a
00852     // right-hand side B, then we didn't set the problem.
00853     if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
00854       isSet_ = false;
00855       return isSet_;
00856     }
00857 
00858     // Reset whether the solution has been updated.  (We're just
00859     // setting the problem now, so of course we haven't updated the
00860     // solution yet.)
00861     solutionUpdated_ = false;
00862     
00863     // Compute the initial residuals.
00864     if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
00865       R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
00866     }
00867     computeCurrResVec( &*R0_, &*X_, &*B_ );
00868 
00869     if (LP_!=Teuchos::null) {
00870       if (PR0_==Teuchos::null || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *B_ )) {
00871         PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
00872       }
00873       {
00874 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00875         Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00876 #endif
00877         OPT::Apply( *LP_, *R0_, *PR0_ );
00878       }
00879     } 
00880     else {
00881       PR0_ = R0_;
00882     }    
00883 
00884     // The problem has been set and is ready for use.
00885     isSet_ = true;
00886     
00887     // Return isSet.
00888     return isSet_;
00889   }
00890   
00891   template <class ScalarType, class MV, class OP>
00892   Teuchos::RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrLHSVec()
00893   {
00894     if (isSet_) {
00895       return curX_;
00896     }
00897     else {
00898       return Teuchos::null;
00899     }
00900   }
00901   
00902   template <class ScalarType, class MV, class OP>
00903   Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getCurrRHSVec()
00904   {
00905     if (isSet_) {
00906       return curB_;
00907     }
00908     else {
00909       return Teuchos::null;
00910     }
00911   }
00912   
00913   template <class ScalarType, class MV, class OP>
00914   void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
00915   {
00916     using Teuchos::null;
00917     using Teuchos::RCP;
00918 
00919     const bool leftPrec = LP_ != null;
00920     const bool rightPrec = RP_ != null;
00921 
00922     // We only need a temporary vector for intermediate results if
00923     // there is a left or right preconditioner.  We really should just
00924     // keep this temporary vector around instead of allocating it each
00925     // time.
00926     RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;
00927 
00928     //
00929     // No preconditioning.
00930     // 
00931     if (!leftPrec && !rightPrec){ 
00932 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00933       Teuchos::TimeMonitor OpTimer(*timerOp_);
00934 #endif
00935       OPT::Apply( *A_, x, y );
00936     }
00937     //
00938     // Preconditioning is being done on both sides
00939     //
00940     else if( leftPrec && rightPrec ) 
00941       {
00942         {
00943 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00944           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00945 #endif
00946     OPT::Apply( *RP_, x, y );   
00947         }
00948         {
00949 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00950           Teuchos::TimeMonitor OpTimer(*timerOp_);
00951 #endif
00952     OPT::Apply( *A_, y, *ytemp );
00953         }
00954         {
00955 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00956           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00957 #endif
00958     OPT::Apply( *LP_, *ytemp, y );
00959         }
00960       }
00961     //
00962     // Preconditioning is only being done on the left side
00963     //
00964     else if( leftPrec ) 
00965       {
00966         {
00967 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00968           Teuchos::TimeMonitor OpTimer(*timerOp_);
00969 #endif
00970     OPT::Apply( *A_, x, *ytemp );
00971         }
00972         {
00973 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00974           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00975 #endif
00976     OPT::Apply( *LP_, *ytemp, y );
00977         }
00978       }
00979     //
00980     // Preconditioning is only being done on the right side
00981     //
00982     else 
00983       {
00984         {
00985 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00986           Teuchos::TimeMonitor PrecTimer(*timerPrec_);
00987 #endif
00988     OPT::Apply( *RP_, x, *ytemp );
00989         }
00990         {
00991 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00992           Teuchos::TimeMonitor OpTimer(*timerOp_);
00993 #endif
00994           OPT::Apply( *A_, *ytemp, y );
00995         }
00996       }  
00997   }
00998   
00999   template <class ScalarType, class MV, class OP>
01000   void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
01001     if (A_.get()) {
01002 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01003       Teuchos::TimeMonitor OpTimer(*timerOp_);
01004 #endif
01005       OPT::Apply( *A_,x, y);   
01006     }
01007     else {
01008       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
01009         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
01010     }
01011   }
01012   
01013   template <class ScalarType, class MV, class OP>
01014   void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
01015     if (LP_!=Teuchos::null) {
01016 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01017       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01018 #endif
01019       return ( OPT::Apply( *LP_,x, y) );
01020     }
01021     else {
01022       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
01023         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
01024     }
01025   }
01026   
01027   template <class ScalarType, class MV, class OP>
01028   void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
01029     if (RP_!=Teuchos::null) {
01030 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01031       Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01032 #endif
01033       return ( OPT::Apply( *RP_,x, y) );
01034     }
01035     else {
01036       MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
01037         Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
01038     }
01039   }
01040   
01041   template <class ScalarType, class MV, class OP>
01042   void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {
01043 
01044     if (R) {
01045       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
01046   {
01047     if (LP_!=Teuchos::null)
01048       {
01049         Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
01050               {
01051 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01052                 Teuchos::TimeMonitor OpTimer(*timerOp_);
01053 #endif
01054           OPT::Apply( *A_, *X, *R_temp );
01055               }
01056         MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
01057               {
01058 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01059                 Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01060 #endif
01061           OPT::Apply( *LP_, *R_temp, *R );
01062               }
01063       }
01064     else 
01065       {
01066               {
01067 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01068                 Teuchos::TimeMonitor OpTimer(*timerOp_);
01069 #endif
01070           OPT::Apply( *A_, *X, *R );
01071               }
01072         MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
01073       }
01074   }
01075       else { 
01076   // The solution and right-hand side may not be specified, check and use which ones exist.
01077   Teuchos::RCP<const MV> localB, localX;
01078   if (B)
01079     localB = Teuchos::rcp( B, false );
01080   else
01081     localB = curB_;
01082   
01083   if (X)
01084     localX = Teuchos::rcp( X, false );
01085   else
01086     localX = curX_;
01087   
01088   if (LP_!=Teuchos::null)
01089     {
01090       Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
01091             {
01092 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01093               Teuchos::TimeMonitor OpTimer(*timerOp_);
01094 #endif
01095         OPT::Apply( *A_, *localX, *R_temp );
01096             }
01097       MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
01098             {
01099 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01100               Teuchos::TimeMonitor PrecTimer(*timerPrec_);
01101 #endif
01102         OPT::Apply( *LP_, *R_temp, *R );
01103             }
01104     }
01105   else 
01106     {
01107             {
01108 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01109               Teuchos::TimeMonitor OpTimer(*timerOp_);
01110 #endif
01111           OPT::Apply( *A_, *localX, *R );
01112             }
01113       MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
01114     }
01115       }    
01116     } 
01117   }
01118   
01119   
01120   template <class ScalarType, class MV, class OP>
01121   void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {
01122 
01123     if (R) {
01124       if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
01125   {
01126           {
01127 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01128             Teuchos::TimeMonitor OpTimer(*timerOp_);
01129 #endif
01130       OPT::Apply( *A_, *X, *R );
01131           }
01132     MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
01133   }
01134       else { 
01135   // The solution and right-hand side may not be specified, check and use which ones exist.
01136   Teuchos::RCP<const MV> localB, localX;
01137   if (B)
01138     localB = Teuchos::rcp( B, false );
01139   else
01140     localB = curB_;
01141   
01142   if (X)
01143     localX = Teuchos::rcp( X, false );
01144   else
01145     localX = curX_;
01146     
01147         {
01148 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01149           Teuchos::TimeMonitor OpTimer(*timerOp_);
01150 #endif
01151     OPT::Apply( *A_, *localX, *R );
01152         }
01153   MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
01154       }    
01155     }
01156   }
01157   
01158 } // end Belos namespace
01159 
01160 #endif /* BELOS_LINEAR_PROBLEM_HPP */
01161 
01162 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines