Belos Version of the Day
BelosPCPGIter.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_PCPG_ITER_HPP
00043 #define BELOS_PCPG_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosMatOrthoManager.hpp"
00054 #include "BelosOutputManager.hpp"
00055 #include "BelosStatusTest.hpp"
00056 #include "BelosOperatorTraits.hpp"
00057 #include "BelosMultiVecTraits.hpp"
00058 
00059 #include "Teuchos_BLAS.hpp"
00060 #include "Teuchos_LAPACK.hpp"
00061 #include "Teuchos_SerialDenseMatrix.hpp"
00062 #include "Teuchos_SerialDenseVector.hpp"
00063 #include "Teuchos_ScalarTraits.hpp"
00064 #include "Teuchos_ParameterList.hpp"
00065 #include "Teuchos_TimeMonitor.hpp"
00066 
00078 namespace Belos {
00079   
00081 
00082   
00087   template <class ScalarType, class MV>
00088   struct PCPGIterState {
00094     int curDim;
00096     int prevUdim;
00097 
00099     Teuchos::RCP<MV> R;
00100 
00102     Teuchos::RCP<MV> Z;
00103 
00105     Teuchos::RCP<MV> P;
00106 
00108     Teuchos::RCP<MV> AP;
00109 
00111     Teuchos::RCP<MV> U;
00112 
00114     Teuchos::RCP<MV> C;
00115 
00120     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
00121 
00122     PCPGIterState() : curDim(0), 
00123                       prevUdim(0), 
00124                       R(Teuchos::null), Z(Teuchos::null), 
00125                       P(Teuchos::null), AP(Teuchos::null),
00126                       U(Teuchos::null), C(Teuchos::null),
00127           D(Teuchos::null)
00128     {}
00129   };
00130   
00132   
00134 
00135   
00147   class PCPGIterInitFailure : public BelosError {public:
00148     PCPGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00149     {}};
00150 
00155   class PCPGIterateFailure : public BelosError {public:
00156     PCPGIterateFailure(const std::string& what_arg) : BelosError(what_arg)
00157     {}};
00158 
00159 
00160   
00167   class PCPGIterOrthoFailure : public BelosError {public:
00168     PCPGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00169     {}};
00170   
00177   class PCPGIterLAPACKFailure : public BelosError {public:
00178     PCPGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00179     {}};
00180   
00182   
00183   
00184   template<class ScalarType, class MV, class OP>
00185   class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
00186     
00187   public:
00188     
00189     //
00190     // Convenience typedefs
00191     //
00192     typedef MultiVecTraits<ScalarType,MV> MVT;
00193     typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00194     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00195     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00196     typedef typename SCT::magnitudeType MagnitudeType;
00197     
00199 
00200     
00208     PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00209     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00210     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00211     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00212     Teuchos::ParameterList &params );
00213     
00215     virtual ~PCPGIter() {};
00217     
00218     
00220 
00221     
00240     void iterate();
00241     
00260     void initialize(PCPGIterState<ScalarType,MV> newstate);
00261     
00265     void initialize()
00266     {
00267       PCPGIterState<ScalarType,MV> empty;
00268       initialize(empty);
00269     }
00270     
00278     PCPGIterState<ScalarType,MV> getState() const {
00279       PCPGIterState<ScalarType,MV> state;
00280       state.Z = Z_;         // CG state
00281       state.P = P_;
00282       state.AP = AP_;
00283       state.R = R_;
00284       state.U = U_;         // seed state 
00285       state.C = C_; 
00286       state.D = D_;
00287       state.curDim = curDim_;
00288       state.prevUdim = prevUdim_;
00289       return state;
00290     }
00291     
00293     
00294     
00296 
00297     
00299     int getNumIters() const { return iter_; }
00300     
00302     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00303     
00306     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00307 
00309 
00312     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00313 
00315     int getCurSubspaceDim() const { 
00316       if (!initialized_) return 0;
00317       return curDim_;
00318     };
00319 
00321     int getPrevSubspaceDim() const { 
00322       if (!initialized_) return 0;
00323       return prevUdim_;
00324     };
00325     
00327     
00328     
00330 
00331     
00333     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00334     
00336     int getBlockSize() const { return 1; }
00337     
00339     int getNumRecycledBlocks() const { return savedBlocks_; }
00340 
00342     
00344     void setBlockSize(int blockSize) {
00345       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00346        "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
00347     }
00348 
00350     void setSize( int savedBlocks );
00351 
00353     bool isInitialized() { return initialized_; }
00354 
00356     void resetState(); 
00357     
00359     
00360   private:
00361     
00362     //
00363     // Internal methods
00364     //
00366     void setStateSize();
00367     
00368     //
00369     // Classes inputed through constructor that define the linear problem to be solved.
00370     //
00371     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00372     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00373     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00374     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00375 
00376     //
00377     // Algorithmic parameters
00378     // savedBlocks_ is the number of blocks allocated for the reused subspace
00379     int savedBlocks_; 
00380     //
00381     // 
00382     // Current solver state
00383     //
00384     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00385     // is capable of running; _initialize is controlled  by the initialize() member method
00386     // For the implications of the state of initialized_, please see documentation for initialize()
00387     bool initialized_;
00388     
00389     // stateStorageInitialized_ indicates that the state storage has be initialized to the current
00390     // savedBlocks_.  State storage initialization may be postponed if the linear problem was
00391     // generated without either the right-hand side or solution vectors.
00392     bool stateStorageInitialized_;
00393 
00394     // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots 
00395     bool keepDiagonal_;
00396 
00397     // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
00398     // out all entries before an iteration is started.
00399     bool initDiagonal_;
00400     
00401     // Current subspace dimension
00402     int curDim_;
00403 
00404     // Dimension of seed space used to solve current linear system
00405     int prevUdim_;
00406     
00407     // Number of iterations performed
00408     int iter_;
00409     // 
00410     // State Storage    ... of course this part is different for CG
00411     //
00412     // Residual
00413     Teuchos::RCP<MV> R_;
00414     //
00415     // Preconditioned residual
00416     Teuchos::RCP<MV> Z_;
00417     //
00418     // Direction std::vector
00419     Teuchos::RCP<MV> P_;
00420     //
00421     // Operator applied to direction std::vector
00422     Teuchos::RCP<MV> AP_;
00423     //
00424     // Recycled subspace vectors.
00425     Teuchos::RCP<MV> U_;
00426     // 
00427     // C = A * U,  linear system is Ax=b 
00428     Teuchos::RCP<MV> C_;
00429     //
00430     // Projected matrices
00431     // D_ : Diagonal matrix of pivots D = P'AP 
00432     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
00433   };
00434   
00436   // Constructor.
00437   template<class ScalarType, class MV, class OP>
00438   PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00439              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00440              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00441              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00442              Teuchos::ParameterList &params ):
00443     lp_(problem),
00444     om_(printer),
00445     stest_(tester),
00446     ortho_(ortho),
00447     savedBlocks_(0),
00448     initialized_(false),
00449     stateStorageInitialized_(false),
00450     keepDiagonal_(false), 
00451     initDiagonal_(false),
00452     curDim_(0),
00453     prevUdim_(0),
00454     iter_(0)
00455   {
00456     // Get the maximum number of blocks allowed for this Krylov subspace
00457 
00458     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
00459                        "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
00460     int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
00461 
00462     // Find out whether we are saving the Diagonal matrix.
00463     keepDiagonal_ = params.get("Keep Diagonal", false);
00464 
00465     // Find out whether we are initializing the Diagonal matrix.
00466     initDiagonal_ = params.get("Initialize Diagonal", false);
00467 
00468     // Set the number of blocks and allocate data
00469     setSize( rb );
00470   }
00471   
00473   // Set the block size and adjust as necessary
00474   template <class ScalarType, class MV, class OP>
00475   void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
00476   {
00477     // allocate space only; perform no computation
00478     // Any change in size invalidates the state of the solver as implemented here.
00479 
00480     TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
00481 
00482     if ( savedBlocks_ != savedBlocks) {
00483       stateStorageInitialized_ = false;
00484       savedBlocks_ = savedBlocks;
00485       initialized_ = false;
00486       curDim_ = 0;
00487       prevUdim_ = 0;
00488       setStateSize(); // Use the current savedBlocks_ to initialize the state storage.    
00489     }
00490   }
00491 
00493   // Enable the reuse of a single solver object for completely different linear systems
00494   template <class ScalarType, class MV, class OP>
00495   void PCPGIter<ScalarType,MV,OP>::resetState()
00496   {
00497       stateStorageInitialized_ = false;
00498       initialized_ = false;
00499       curDim_ = 0;
00500       prevUdim_ = 0;
00501       setStateSize();
00502   }
00503 
00505   // Setup the state storage.  Called by either initialize or, if savedBlocks_ changes, setSize.
00506   template <class ScalarType, class MV, class OP>
00507   void PCPGIter<ScalarType,MV,OP>::setStateSize ()
00508   {
00509     if (!stateStorageInitialized_) {
00510 
00511       // Check if there is any multivector to clone from.
00512       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00513       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00514       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00515   return;  // postpone exception 
00516       }
00517       else {
00518   
00520   // blockSize*recycledBlocks dependent
00521         int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
00522   //
00523   // Initialize the CG state storage
00524         // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
00525   // Generate CG state only if it does not exist, otherwise resize it.
00526         if (Z_ == Teuchos::null) {
00527           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00528           Z_ = MVT::Clone( *tmp, 1 );
00529         }
00530         if (P_ == Teuchos::null) {
00531           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00532           P_ = MVT::Clone( *tmp, 1 );
00533         }
00534         if (AP_ == Teuchos::null) {
00535           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00536           AP_ = MVT::Clone( *tmp, 1 );
00537         }
00538 
00539   if (C_ == Teuchos::null) {        
00540 
00541     // Get the multivector that is not null. 
00542     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00543     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00544            "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00545     TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00546            "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
00547     C_ = MVT::Clone( *tmp, savedBlocks_ );
00548   }
00549   else {
00550     // Generate C_ by cloning itself ONLY if more space is needed.
00551     if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
00552       Teuchos::RCP<const MV> tmp = C_;
00553       C_ = MVT::Clone( *tmp, savedBlocks_ );
00554     }
00555   }
00556   if (U_ == Teuchos::null) {        
00557     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00558     TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00559            "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
00560     U_ = MVT::Clone( *tmp, savedBlocks_ );
00561   }
00562   else {
00563     // Generate U_ by cloning itself ONLY if more space is needed.
00564     if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
00565       Teuchos::RCP<const MV> tmp = U_;
00566       U_ = MVT::Clone( *tmp, savedBlocks_ );
00567     }
00568   }
00569         if (keepDiagonal_) {
00570           if (D_ == Teuchos::null) {
00571             D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00572           }
00573           if (initDiagonal_) {
00574             D_->shape( newsd, newsd );
00575           }
00576           else {
00577             if (D_->numRows() < newsd || D_->numCols() < newsd) {
00578               D_->shapeUninitialized( newsd, newsd );
00579             }
00580           }
00581         }
00582   // State storage has now been initialized.
00583   stateStorageInitialized_ = true;
00584       } // if there is a vector to clone from
00585     } // if !stateStorageInitialized_ 
00586   } // end of setStateSize
00587 
00589   // Initialize the iteration object
00590   template <class ScalarType, class MV, class OP>
00591   void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV> newstate)
00592   {
00593 
00594     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00595            "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
00596     
00597     // Requirements: R_ and consistent multivectors widths and lengths
00598     //
00599     std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
00600 
00601     if (newstate.R != Teuchos::null){ 
00602 
00603       R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
00604       if (newstate.U == Teuchos::null){ 
00605         prevUdim_ = 0;
00606         newstate.U = U_;
00607         newstate.C = C_;
00608       }
00609       else {
00610         prevUdim_ =  newstate.curDim;
00611         if (newstate.C == Teuchos::null){  // Stub for new feature
00612           std::vector<int> index(prevUdim_);
00613           for (int i=0; i< prevUdim_; ++i)  
00614             index[i] = i; 
00615           Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
00616           newstate.C = MVT::Clone( *newstate.U, prevUdim_ ); 
00617           Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );    
00618           lp_->apply( *Ukeff, *Ckeff );
00619         }
00620         curDim_ = prevUdim_ ;
00621       }
00622 
00623       // Initialize the state storage if not already allocated in the constructor
00624       if (!stateStorageInitialized_) 
00625         setStateSize();
00626 
00627       //TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.V) != MVText::GetGlobalLength(*V_), std::invalid_argument, errstr );
00628       //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
00629 
00630       newstate.prevUdim =  prevUdim_; // big change in functionality from GCRODR 
00631       newstate.curDim =  curDim_; 
00632 
00633       //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
00634 
00635       std::vector<int> zero_index(1);
00636       zero_index[0] = 0;
00637       if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction 
00638         lp_->applyLeftPrec( *R_, *Z_ );
00639         MVT::SetBlock( *Z_,  zero_index , *P_ );  // P(:,zero_index) := Z
00640       } else {                                      
00641         Z_ = R_;
00642         MVT::SetBlock( *R_, zero_index, *P_ );
00643       }
00644 
00645       std::vector<int> nextind(1);
00646       nextind[0] = curDim_;
00647 
00648       MVT::SetBlock( *P_,  nextind, *newstate.U ); // U(:,curDim_ ) := P_
00649 
00650       ++curDim_;
00651       newstate.curDim = curDim_; 
00652 
00653       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
00654                           std::invalid_argument, errstr );
00655       if (newstate.U != U_) { // Why this is needed?
00656   U_ = newstate.U;
00657       }
00658 
00659       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
00660                           std::invalid_argument, errstr );
00661       if (newstate.C != C_) {
00662   C_ = newstate.C;
00663       }
00664     }
00665     else {
00666 
00667       TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00668                          "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
00669     }
00670 
00671     // the solver is initialized
00672     initialized_ = true;
00673   }
00674 
00675 
00677   // Iterate until the status test informs us we should stop.
00678   template <class ScalarType, class MV, class OP>
00679   void PCPGIter<ScalarType,MV,OP>::iterate()
00680   {
00681     //
00682     // Allocate/initialize data structures
00683     //
00684     if (initialized_ == false) {
00685       initialize();
00686     }
00687     const bool debug = false;
00688 
00689     // Allocate memory for scalars.
00690     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00691     Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
00692     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00693     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
00694 
00695     if( iter_ != 0 )
00696       std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl;  //DMD
00697 
00698     // GenOrtho Project Stubs
00699     std::vector<int> prevInd;
00700     Teuchos::RCP<const MV> Uprev;
00701     Teuchos::RCP<const MV> Cprev;
00702     Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
00703 
00704     if( prevUdim_ ){
00705       prevInd.resize( prevUdim_ );
00706       for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
00707       CZ.reshape( prevUdim_ , 1 );
00708       Uprev = MVT::CloneView(*U_, prevInd);
00709       Cprev = MVT::CloneView(*C_, prevInd);
00710     }
00711 
00712     // Get the current solution std::vector.
00713     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00714 
00715     // Check that the current solution std::vector only has one column.
00716     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
00717                         "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
00718 
00719     //Check that the input is correctly set up 
00720     TEUCHOS_TEST_FOR_EXCEPTION( curDim_  != prevUdim_ + 1, PCPGIterInitFailure,
00721                         "Belos::CGIter::iterate(): mistake in initialization !" );
00722 
00723 
00724     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00725     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00726 
00727 
00728     std::vector<int> curind(1);
00729     std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
00730     if (prevUdim_ > 0){                 // A-orthonalize P=Z to Uprev
00731       Teuchos::RCP<MV> P; 
00732       curind[0] = curDim_ - 1;          // column = dimension - 1 
00733       P = MVT::CloneViewNonConst(*U_,curind); 
00734       MVT::MvTransMv( one, *Cprev, *P, CZ );
00735       MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );       // P -= U*(C'Z)
00736 
00737       if( debug ){
00738         MVT::MvTransMv( one, *Cprev, *P, CZ );
00739         std::cout << " Input CZ post ortho " << std::endl;
00740         CZ.print( std::cout );
00741       }
00742       if( curDim_ == savedBlocks_ ){
00743         std::vector<int> zero_index(1);
00744         zero_index[0] = 0;
00745         MVT::SetBlock( *P, zero_index, *P_ );
00746       }
00747       P = Teuchos::null;
00748     }
00749 
00750     // Compute first <r,z> a.k.a. rHz
00751     MVT::MvTransMv( one, *R_, *Z_, rHz );
00752 
00754     // iterate until the status test is satisfied
00755     //
00756     while (stest_->checkStatus(this) != Passed ) {
00757       Teuchos::RCP<const MV> P; 
00758       Teuchos::RCP<MV> AP;
00759       iter_++;                          // The next iteration begins.
00760       //std::vector<int> curind(1);
00761       curind[0] = curDim_ - 1;          // column = dimension - 1 
00762       if( debug ){
00763         MVT::MvNorm(*R_, rnorm);
00764         std::cout << iter_ << "  " << curDim_ <<  "   " << rnorm[0] << std::endl;
00765       }
00766       if( prevUdim_ + iter_ < savedBlocks_ ){
00767         P = MVT::CloneView(*U_,curind); 
00768         AP = MVT::CloneViewNonConst(*C_,curind); 
00769         lp_->applyOp( *P, *AP );
00770         MVT::MvTransMv( one, *P, *AP, pAp );
00771       }else{
00772         if( prevUdim_ + iter_ == savedBlocks_ ){
00773           AP = MVT::CloneViewNonConst(*C_,curind); 
00774           lp_->applyOp( *P_, *AP );
00775           MVT::MvTransMv( one, *P_, *AP, pAp );
00776         }else{
00777           lp_->applyOp( *P_, *AP_ );
00778           MVT::MvTransMv( one, *P_, *AP_, pAp );
00779         }
00780       }
00781 
00782       if( keepDiagonal_  && prevUdim_ + iter_ <= savedBlocks_ )
00783         (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
00784 
00785       // positive pAp required 
00786       TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, PCPGIterateFailure,
00787                           "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00788 
00789       // alpha := <R_,Z_> / <P,AP>
00790       alpha(0,0) = rHz(0,0) / pAp(0,0);
00791 
00792       // positive alpha required 
00793       TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
00794                           "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
00795 
00796       // solution update  x += alpha * P
00797       if( curDim_ < savedBlocks_ ){
00798          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
00799       }else{
00800          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00801       }
00802       //lp_->updateSolution(); ... does nothing.
00803       //
00804       // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
00805       //
00806       rHz_old(0,0) = rHz(0,0);
00807       //
00808       // residual update R_ := R_ - alpha * AP
00809       //
00810       if( prevUdim_ + iter_ <= savedBlocks_ ){
00811          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
00812          AP = Teuchos::null;
00813       }else{
00814          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00815       }
00816       //
00817       // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
00818       //
00819       if ( lp_->getLeftPrec() != Teuchos::null ) {
00820         lp_->applyLeftPrec( *R_, *Z_ );
00821       } else {
00822         Z_ = R_;
00823       }
00824       //
00825       MVT::MvTransMv( one, *R_, *Z_, rHz );
00826       //
00827       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00828       //
00829       if( curDim_ < savedBlocks_ ){
00830          curDim_++;                                                         // update basis dim
00831          curind[0] = curDim_ - 1;
00832          Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
00833          MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
00834          if( prevUdim_ ){ // Deflate seed space 
00835              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00836              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
00837              if( debug ){
00838                std::cout << " Check CZ " << std::endl;
00839                MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
00840                CZ.print( std::cout );
00841              }
00842          }
00843          P = Teuchos::null;
00844          if( curDim_ == savedBlocks_ ){
00845            std::vector<int> zero_index(1);
00846            zero_index[0] = 0;
00847            MVT::SetBlock( *Pnext, zero_index, *P_ );
00848          }
00849          Pnext = Teuchos::null;
00850       }else{
00851          MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00852          if( prevUdim_ ){ // Deflate seed space
00853              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00854              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );       // P_ -= U*(C'Z)
00855 
00856              if( debug ){
00857                std::cout << " Check CZ " << std::endl;
00858                MVT::MvTransMv( one, *Cprev, *P_, CZ );
00859                CZ.print( std::cout );
00860              }
00861          }
00862       }
00863       // CGB: 5/26/2010
00864       // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
00865       // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
00866       // to ensure that this wasn't a bug, I verify below that we have set P == null, i.e., that we are not going to use it again
00867       // same for AP
00868       TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
00869     } // end coupled two-term recursion
00870     if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
00871   }
00872 
00873 } // end Belos namespace
00874 
00875 #endif /* BELOS_PCPG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines