Belos Package Browser (Single Doxygen Collection) Development
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 OperatorTraits<ScalarType,MV,OP> OPT;
00194     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00195     typedef typename SCT::magnitudeType MagnitudeType;
00196     
00198 
00199     
00207     PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00208     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00209     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00210     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00211     Teuchos::ParameterList &params );
00212     
00214     virtual ~PCPGIter() {};
00216     
00217     
00219 
00220     
00239     void iterate();
00240     
00259     void initialize(PCPGIterState<ScalarType,MV> newstate);
00260     
00264     void initialize()
00265     {
00266       PCPGIterState<ScalarType,MV> empty;
00267       initialize(empty);
00268     }
00269     
00277     PCPGIterState<ScalarType,MV> getState() const {
00278       PCPGIterState<ScalarType,MV> state;
00279       state.Z = Z_;         // CG state
00280       state.P = P_;
00281       state.AP = AP_;
00282       state.R = R_;
00283       state.U = U_;         // seed state 
00284       state.C = C_; 
00285       state.D = D_;
00286       state.curDim = curDim_;
00287       state.prevUdim = prevUdim_;
00288       return state;
00289     }
00290     
00292     
00293     
00295 
00296     
00298     int getNumIters() const { return iter_; }
00299     
00301     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00302     
00305     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00306 
00308 
00311     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00312 
00314     int getCurSubspaceDim() const { 
00315       if (!initialized_) return 0;
00316       return curDim_;
00317     };
00318 
00320     int getPrevSubspaceDim() const { 
00321       if (!initialized_) return 0;
00322       return prevUdim_;
00323     };
00324     
00326     
00327     
00329 
00330     
00332     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00333     
00335     int getBlockSize() const { return 1; }
00336     
00338     int getNumRecycledBlocks() const { return savedBlocks_; }
00339 
00341     
00343     void setBlockSize(int blockSize) {
00344       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00345        "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
00346     }
00347 
00349     void setSize( int savedBlocks );
00350 
00352     bool isInitialized() { return initialized_; }
00353 
00355     void resetState(); 
00356     
00358     
00359   private:
00360     
00361     //
00362     // Internal methods
00363     //
00365     void setStateSize();
00366     
00367     //
00368     // Classes inputed through constructor that define the linear problem to be solved.
00369     //
00370     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00371     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00372     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00373     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00374 
00375     //
00376     // Algorithmic parameters
00377     // savedBlocks_ is the number of blocks allocated for the reused subspace
00378     int savedBlocks_; 
00379     //
00380     // 
00381     // Current solver state
00382     //
00383     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00384     // is capable of running; _initialize is controlled  by the initialize() member method
00385     // For the implications of the state of initialized_, please see documentation for initialize()
00386     bool initialized_;
00387     
00388     // stateStorageInitialized_ indicates that the state storage has be initialized to the current
00389     // savedBlocks_.  State storage initialization may be postponed if the linear problem was
00390     // generated without either the right-hand side or solution vectors.
00391     bool stateStorageInitialized_;
00392 
00393     // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots 
00394     bool keepDiagonal_;
00395 
00396     // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
00397     // out all entries before an iteration is started.
00398     bool initDiagonal_;
00399     
00400     // Current subspace dimension
00401     int curDim_;
00402 
00403     // Dimension of seed space used to solve current linear system
00404     int prevUdim_;
00405     
00406     // Number of iterations performed
00407     int iter_;
00408     // 
00409     // State Storage    ... of course this part is different for CG
00410     //
00411     // Residual
00412     Teuchos::RCP<MV> R_;
00413     //
00414     // Preconditioned residual
00415     Teuchos::RCP<MV> Z_;
00416     //
00417     // Direction std::vector
00418     Teuchos::RCP<MV> P_;
00419     //
00420     // Operator applied to direction std::vector
00421     Teuchos::RCP<MV> AP_;
00422     //
00423     // Recycled subspace vectors.
00424     Teuchos::RCP<MV> U_;
00425     // 
00426     // C = A * U,  linear system is Ax=b 
00427     Teuchos::RCP<MV> C_;
00428     //
00429     // Projected matrices
00430     // D_ : Diagonal matrix of pivots D = P'AP 
00431     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
00432   };
00433   
00435   // Constructor.
00436   template<class ScalarType, class MV, class OP>
00437   PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00438              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00439              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00440              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00441              Teuchos::ParameterList &params ):
00442     lp_(problem),
00443     om_(printer),
00444     stest_(tester),
00445     ortho_(ortho),
00446     savedBlocks_(0),
00447     initialized_(false),
00448     stateStorageInitialized_(false),
00449     keepDiagonal_(false), 
00450     initDiagonal_(false),
00451     curDim_(0),
00452     prevUdim_(0),
00453     iter_(0)
00454   {
00455     // Get the maximum number of blocks allowed for this Krylov subspace
00456 
00457     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
00458                        "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
00459     int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
00460 
00461     // Find out whether we are saving the Diagonal matrix.
00462     keepDiagonal_ = params.get("Keep Diagonal", false);
00463 
00464     // Find out whether we are initializing the Diagonal matrix.
00465     initDiagonal_ = params.get("Initialize Diagonal", false);
00466 
00467     // Set the number of blocks and allocate data
00468     setSize( rb );
00469   }
00470   
00472   // Set the block size and adjust as necessary
00473   template <class ScalarType, class MV, class OP>
00474   void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
00475   {
00476     // allocate space only; perform no computation
00477     // Any change in size invalidates the state of the solver as implemented here.
00478 
00479     TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
00480 
00481     if ( savedBlocks_ != savedBlocks) {
00482       stateStorageInitialized_ = false;
00483       savedBlocks_ = savedBlocks;
00484       initialized_ = false;
00485       curDim_ = 0;
00486       prevUdim_ = 0;
00487       setStateSize(); // Use the current savedBlocks_ to initialize the state storage.    
00488     }
00489   }
00490 
00492   // Enable the reuse of a single solver object for completely different linear systems
00493   template <class ScalarType, class MV, class OP>
00494   void PCPGIter<ScalarType,MV,OP>::resetState()
00495   {
00496       stateStorageInitialized_ = false;
00497       initialized_ = false;
00498       curDim_ = 0;
00499       prevUdim_ = 0;
00500       setStateSize();
00501   }
00502 
00504   // Setup the state storage.  Called by either initialize or, if savedBlocks_ changes, setSize.
00505   template <class ScalarType, class MV, class OP>
00506   void PCPGIter<ScalarType,MV,OP>::setStateSize ()
00507   {
00508     if (!stateStorageInitialized_) {
00509 
00510       // Check if there is any multivector to clone from.
00511       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00512       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00513       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00514   return;  // postpone exception 
00515       }
00516       else {
00517   
00519   // blockSize*recycledBlocks dependent
00520         int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
00521   //
00522   // Initialize the CG state storage
00523         // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
00524   // Generate CG state only if it does not exist, otherwise resize it.
00525         if (Z_ == Teuchos::null) {
00526           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00527           Z_ = MVT::Clone( *tmp, 1 );
00528         }
00529         if (P_ == Teuchos::null) {
00530           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00531           P_ = MVT::Clone( *tmp, 1 );
00532         }
00533         if (AP_ == Teuchos::null) {
00534           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00535           AP_ = MVT::Clone( *tmp, 1 );
00536         }
00537 
00538   if (C_ == Teuchos::null) {        
00539 
00540     // Get the multivector that is not null. 
00541     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00542     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00543            "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00544     TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00545            "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
00546     C_ = MVT::Clone( *tmp, savedBlocks_ );
00547   }
00548   else {
00549     // Generate C_ by cloning itself ONLY if more space is needed.
00550     if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
00551       Teuchos::RCP<const MV> tmp = C_;
00552       C_ = MVT::Clone( *tmp, savedBlocks_ );
00553     }
00554   }
00555   if (U_ == Teuchos::null) {        
00556     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00557     TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00558            "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
00559     U_ = MVT::Clone( *tmp, savedBlocks_ );
00560   }
00561   else {
00562     // Generate U_ by cloning itself ONLY if more space is needed.
00563     if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
00564       Teuchos::RCP<const MV> tmp = U_;
00565       U_ = MVT::Clone( *tmp, savedBlocks_ );
00566     }
00567   }
00568         if (keepDiagonal_) {
00569           if (D_ == Teuchos::null) {
00570             D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00571           }
00572           if (initDiagonal_) {
00573             D_->shape( newsd, newsd );
00574           }
00575           else {
00576             if (D_->numRows() < newsd || D_->numCols() < newsd) {
00577               D_->shapeUninitialized( newsd, newsd );
00578             }
00579           }
00580         }
00581   // State storage has now been initialized.
00582   stateStorageInitialized_ = true;
00583       } // if there is a vector to clone from
00584     } // if !stateStorageInitialized_ 
00585   } // end of setStateSize
00586 
00588   // Initialize the iteration object
00589   template <class ScalarType, class MV, class OP>
00590   void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV> newstate)
00591   {
00592 
00593     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00594            "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
00595     
00596     // Requirements: R_ and consistent multivectors widths and lengths
00597     //
00598     std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
00599 
00600     if (newstate.R != Teuchos::null){ 
00601 
00602       R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
00603       if (newstate.U == Teuchos::null){ 
00604         prevUdim_ = 0;
00605         newstate.U = U_;
00606         newstate.C = C_;
00607       }
00608       else {
00609         prevUdim_ =  newstate.curDim;
00610         if (newstate.C == Teuchos::null){  // Stub for new feature
00611           std::vector<int> index(prevUdim_);
00612           for (int i=0; i< prevUdim_; ++i)  
00613             index[i] = i; 
00614           Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
00615           newstate.C = MVT::Clone( *newstate.U, prevUdim_ ); 
00616           Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );    
00617           lp_->apply( *Ukeff, *Ckeff );
00618         }
00619         curDim_ = prevUdim_ ;
00620       }
00621 
00622       // Initialize the state storage if not already allocated in the constructor
00623       if (!stateStorageInitialized_) 
00624         setStateSize();
00625 
00626       //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_), std::invalid_argument, errstr );
00627       //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
00628 
00629       newstate.prevUdim =  prevUdim_; // big change in functionality from GCRODR 
00630       newstate.curDim =  curDim_; 
00631 
00632       //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
00633 
00634       std::vector<int> zero_index(1);
00635       zero_index[0] = 0;
00636       if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction 
00637         lp_->applyLeftPrec( *R_, *Z_ );
00638         MVT::SetBlock( *Z_,  zero_index , *P_ );  // P(:,zero_index) := Z
00639       } else {                                      
00640         Z_ = R_;
00641         MVT::SetBlock( *R_, zero_index, *P_ );
00642       }
00643 
00644       std::vector<int> nextind(1);
00645       nextind[0] = curDim_;
00646 
00647       MVT::SetBlock( *P_,  nextind, *newstate.U ); // U(:,curDim_ ) := P_
00648 
00649       ++curDim_;
00650       newstate.curDim = curDim_; 
00651 
00652       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
00653                           std::invalid_argument, errstr );
00654       if (newstate.U != U_) { // Why this is needed?
00655   U_ = newstate.U;
00656       }
00657 
00658       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
00659                           std::invalid_argument, errstr );
00660       if (newstate.C != C_) {
00661   C_ = newstate.C;
00662       }
00663     }
00664     else {
00665 
00666       TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00667                          "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
00668     }
00669 
00670     // the solver is initialized
00671     initialized_ = true;
00672   }
00673 
00674 
00676   // Iterate until the status test informs us we should stop.
00677   template <class ScalarType, class MV, class OP>
00678   void PCPGIter<ScalarType,MV,OP>::iterate()
00679   {
00680     //
00681     // Allocate/initialize data structures
00682     //
00683     if (initialized_ == false) {
00684       initialize();
00685     }
00686     bool debug = false;
00687 
00688     // Allocate memory for scalars.
00689     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00690     Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
00691     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00692     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
00693 
00694     if( iter_ != 0 )
00695       std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl;  //DMD
00696 
00697     // GenOrtho Project Stubs
00698     std::vector<int> prevInd;
00699     Teuchos::RCP<const MV> Uprev;
00700     Teuchos::RCP<const MV> Cprev;
00701     Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
00702 
00703     if( prevUdim_ ){
00704       prevInd.resize( prevUdim_ );
00705       for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
00706       CZ.reshape( prevUdim_ , 1 );
00707       Uprev = MVT::CloneView(*U_, prevInd);
00708       Cprev = MVT::CloneView(*C_, prevInd);
00709     }
00710 
00711     // Get the current solution std::vector.
00712     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00713 
00714     // Check that the current solution std::vector only has one column.
00715     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
00716                         "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
00717 
00718     //Check that the input is correctly set up 
00719     TEUCHOS_TEST_FOR_EXCEPTION( curDim_  != prevUdim_ + 1, PCPGIterInitFailure,
00720                         "Belos::CGIter::iterate(): mistake in initialization !" );
00721 
00722 
00723     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00724     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00725 
00726 
00727     std::vector<int> curind(1);
00728     std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
00729     if (prevUdim_ > 0){                 // A-orthonalize P=Z to Uprev
00730       Teuchos::RCP<MV> P; 
00731       curind[0] = curDim_ - 1;          // column = dimension - 1 
00732       P = MVT::CloneViewNonConst(*U_,curind); 
00733       MVT::MvTransMv( one, *Cprev, *P, CZ );
00734       MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );       // P -= U*(C'Z)
00735 
00736       if( debug ){
00737         MVT::MvTransMv( one, *Cprev, *P, CZ );
00738         std::cout << " Input CZ post ortho " << std::endl;
00739         CZ.print( std::cout );
00740       }
00741       if( curDim_ == savedBlocks_ ){
00742         std::vector<int> zero_index(1);
00743         zero_index[0] = 0;
00744         MVT::SetBlock( *P, zero_index, *P_ );
00745       }
00746       P = Teuchos::null;
00747     }
00748 
00749     // Compute first <r,z> a.k.a. rHz
00750     MVT::MvTransMv( one, *R_, *Z_, rHz );
00751 
00753     // iterate until the status test is satisfied
00754     //
00755     while (stest_->checkStatus(this) != Passed ) {
00756       Teuchos::RCP<const MV> P; 
00757       Teuchos::RCP<MV> AP;
00758       iter_++;                          // The next iteration begins.
00759       //std::vector<int> curind(1);
00760       curind[0] = curDim_ - 1;          // column = dimension - 1 
00761       if( debug ){
00762         MVT::MvNorm(*R_, rnorm);
00763         std::cout << iter_ << "  " << curDim_ <<  "   " << rnorm[0] << std::endl;
00764       }
00765       if( prevUdim_ + iter_ < savedBlocks_ ){
00766         P = MVT::CloneView(*U_,curind); 
00767         AP = MVT::CloneViewNonConst(*C_,curind); 
00768         lp_->applyOp( *P, *AP );
00769         MVT::MvTransMv( one, *P, *AP, pAp );
00770       }else{
00771         if( prevUdim_ + iter_ == savedBlocks_ ){
00772           AP = MVT::CloneViewNonConst(*C_,curind); 
00773           lp_->applyOp( *P_, *AP );
00774           MVT::MvTransMv( one, *P_, *AP, pAp );
00775         }else{
00776           lp_->applyOp( *P_, *AP_ );
00777           MVT::MvTransMv( one, *P_, *AP_, pAp );
00778         }
00779       }
00780 
00781       if( keepDiagonal_  && prevUdim_ + iter_ <= savedBlocks_ )
00782         (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
00783 
00784       // positive pAp required 
00785       TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, PCPGIterateFailure,
00786                           "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00787 
00788       // alpha := <R_,Z_> / <P,AP>
00789       alpha(0,0) = rHz(0,0) / pAp(0,0);
00790 
00791       // positive alpha required 
00792       TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
00793                           "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
00794 
00795       // solution update  x += alpha * P
00796       if( curDim_ < savedBlocks_ ){
00797          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
00798       }else{
00799          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00800       }
00801       //lp_->updateSolution(); ... does nothing.
00802       //
00803       // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
00804       //
00805       rHz_old(0,0) = rHz(0,0);
00806       //
00807       // residual update R_ := R_ - alpha * AP
00808       //
00809       if( prevUdim_ + iter_ <= savedBlocks_ ){
00810          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
00811          AP = Teuchos::null;
00812       }else{
00813          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00814       }
00815       //
00816       // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
00817       //
00818       if ( lp_->getLeftPrec() != Teuchos::null ) {
00819         lp_->applyLeftPrec( *R_, *Z_ );
00820       } else {
00821         Z_ = R_;
00822       }
00823       //
00824       MVT::MvTransMv( one, *R_, *Z_, rHz );
00825       //
00826       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00827       //
00828       if( curDim_ < savedBlocks_ ){
00829          curDim_++;                                                         // update basis dim
00830          curind[0] = curDim_ - 1;
00831          Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
00832          MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
00833          if( prevUdim_ ){ // Deflate seed space 
00834              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00835              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
00836              if( debug ){
00837                std::cout << " Check CZ " << std::endl;
00838                MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
00839                CZ.print( std::cout );
00840              }
00841          }
00842          P = Teuchos::null;
00843          if( curDim_ == savedBlocks_ ){
00844            std::vector<int> zero_index(1);
00845            zero_index[0] = 0;
00846            MVT::SetBlock( *Pnext, zero_index, *P_ );
00847          }
00848          Pnext = Teuchos::null;
00849       }else{
00850          MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00851          if( prevUdim_ ){ // Deflate seed space
00852              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00853              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );       // P_ -= U*(C'Z)
00854 
00855              if( debug ){
00856                std::cout << " Check CZ " << std::endl;
00857                MVT::MvTransMv( one, *Cprev, *P_, CZ );
00858                CZ.print( std::cout );
00859              }
00860          }
00861       }
00862       // CGB: 5/26/2010
00863       // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
00864       // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
00865       // 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
00866       // same for AP
00867       TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
00868     } // end coupled two-term recursion
00869     if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
00870   }
00871 
00872 } // end Belos namespace
00873 
00874 #endif /* BELOS_PCPG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines