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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_PCPG_ITER_HPP
00030 #define BELOS_PCPG_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosMatOrthoManager.hpp"
00041 #include "BelosOutputManager.hpp"
00042 #include "BelosStatusTest.hpp"
00043 #include "BelosOperatorTraits.hpp"
00044 #include "BelosMultiVecTraits.hpp"
00045 
00046 #include "Teuchos_BLAS.hpp"
00047 #include "Teuchos_LAPACK.hpp"
00048 #include "Teuchos_SerialDenseMatrix.hpp"
00049 #include "Teuchos_SerialDenseVector.hpp"
00050 #include "Teuchos_ScalarTraits.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00065 namespace Belos {
00066   
00068 
00069   
00074   template <class ScalarType, class MV>
00075   struct PCPGIterState {
00081     int curDim;
00083     int prevUdim;
00084 
00086     Teuchos::RCP<MV> R;
00087 
00089     Teuchos::RCP<MV> Z;
00090 
00092     Teuchos::RCP<MV> P;
00093 
00095     Teuchos::RCP<MV> AP;
00096 
00098     Teuchos::RCP<MV> U;
00099 
00101     Teuchos::RCP<MV> C;
00102 
00107     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
00108 
00109     PCPGIterState() : curDim(0), 
00110                       prevUdim(0), 
00111                       R(Teuchos::null), Z(Teuchos::null), 
00112                       P(Teuchos::null), AP(Teuchos::null),
00113                       U(Teuchos::null), C(Teuchos::null),
00114           D(Teuchos::null)
00115     {}
00116   };
00117   
00119   
00121 
00122   
00134   class PCPGIterInitFailure : public BelosError {public:
00135     PCPGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00136     {}};
00137 
00142   class PCPGIterateFailure : public BelosError {public:
00143     PCPGIterateFailure(const std::string& what_arg) : BelosError(what_arg)
00144     {}};
00145 
00146 
00147   
00154   class PCPGIterOrthoFailure : public BelosError {public:
00155     PCPGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00156     {}};
00157   
00164   class PCPGIterLAPACKFailure : public BelosError {public:
00165     PCPGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00166     {}};
00167   
00169   
00170   
00171   template<class ScalarType, class MV, class OP>
00172   class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
00173     
00174   public:
00175     
00176     //
00177     // Convenience typedefs
00178     //
00179     typedef MultiVecTraits<ScalarType,MV> MVT;
00180     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00181     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00182     typedef typename SCT::magnitudeType MagnitudeType;
00183     
00185 
00186     
00194     PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00195     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00196     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00197     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00198     Teuchos::ParameterList &params );
00199     
00201     virtual ~PCPGIter() {};
00203     
00204     
00206 
00207     
00226     void iterate();
00227     
00246     void initialize(PCPGIterState<ScalarType,MV> newstate);
00247     
00251     void initialize()
00252     {
00253       PCPGIterState<ScalarType,MV> empty;
00254       initialize(empty);
00255     }
00256     
00264     PCPGIterState<ScalarType,MV> getState() const {
00265       PCPGIterState<ScalarType,MV> state;
00266       state.Z = Z_;         // CG state
00267       state.P = P_;
00268       state.AP = AP_;
00269       state.R = R_;
00270       state.U = U_;         // seed state 
00271       state.C = C_; 
00272       state.D = D_;
00273       state.curDim = curDim_;
00274       state.prevUdim = prevUdim_;
00275       return state;
00276     }
00277     
00279     
00280     
00282 
00283     
00285     int getNumIters() const { return iter_; }
00286     
00288     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00289     
00292     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00293 
00295 
00298     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00299 
00301     int getCurSubspaceDim() const { 
00302       if (!initialized_) return 0;
00303       return curDim_;
00304     };
00305 
00307     int getPrevSubspaceDim() const { 
00308       if (!initialized_) return 0;
00309       return prevUdim_;
00310     };
00311     
00313     
00314     
00316 
00317     
00319     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00320     
00322     int getBlockSize() const { return 1; }
00323     
00325     int getNumRecycledBlocks() const { return savedBlocks_; }
00326 
00328     
00330     void setBlockSize(int blockSize) {
00331       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00332        "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
00333     }
00334 
00336     void setSize( int savedBlocks );
00337 
00339     bool isInitialized() { return initialized_; }
00340 
00342     void resetState(); 
00343     
00345     
00346   private:
00347     
00348     //
00349     // Internal methods
00350     //
00352     void setStateSize();
00353     
00354     //
00355     // Classes inputed through constructor that define the linear problem to be solved.
00356     //
00357     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00358     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00359     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00360     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00361 
00362     //
00363     // Algorithmic parameters
00364     // savedBlocks_ is the number of blocks allocated for the reused subspace
00365     int savedBlocks_; 
00366     //
00367     // 
00368     // Current solver state
00369     //
00370     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00371     // is capable of running; _initialize is controlled  by the initialize() member method
00372     // For the implications of the state of initialized_, please see documentation for initialize()
00373     bool initialized_;
00374     
00375     // stateStorageInitialized_ indicates that the state storage has be initialized to the current
00376     // savedBlocks_.  State storage initialization may be postponed if the linear problem was
00377     // generated without either the right-hand side or solution vectors.
00378     bool stateStorageInitialized_;
00379 
00380     // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots 
00381     bool keepDiagonal_;
00382 
00383     // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
00384     // out all entries before an iteration is started.
00385     bool initDiagonal_;
00386     
00387     // Current subspace dimension
00388     int curDim_;
00389 
00390     // Dimension of seed space used to solve current linear system
00391     int prevUdim_;
00392     
00393     // Number of iterations performed
00394     int iter_;
00395     // 
00396     // State Storage    ... of course this part is different for CG
00397     //
00398     // Residual
00399     Teuchos::RCP<MV> R_;
00400     //
00401     // Preconditioned residual
00402     Teuchos::RCP<MV> Z_;
00403     //
00404     // Direction std::vector
00405     Teuchos::RCP<MV> P_;
00406     //
00407     // Operator applied to direction std::vector
00408     Teuchos::RCP<MV> AP_;
00409     //
00410     // Recycled subspace vectors.
00411     Teuchos::RCP<MV> U_;
00412     // 
00413     // C = A * U,  linear system is Ax=b 
00414     Teuchos::RCP<MV> C_;
00415     //
00416     // Projected matrices
00417     // D_ : Diagonal matrix of pivots D = P'AP 
00418     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
00419   };
00420   
00422   // Constructor.
00423   template<class ScalarType, class MV, class OP>
00424   PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00425              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00426              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00427              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00428              Teuchos::ParameterList &params ):
00429     lp_(problem),
00430     om_(printer),
00431     stest_(tester),
00432     ortho_(ortho),
00433     savedBlocks_(0),
00434     initialized_(false),
00435     stateStorageInitialized_(false),
00436     keepDiagonal_(false), 
00437     initDiagonal_(false),
00438     curDim_(0),
00439     prevUdim_(0),
00440     iter_(0)
00441   {
00442     // Get the maximum number of blocks allowed for this Krylov subspace
00443 
00444     TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
00445                        "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
00446     int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
00447 
00448     // Find out whether we are saving the Diagonal matrix.
00449     keepDiagonal_ = params.get("Keep Diagonal", false);
00450 
00451     // Find out whether we are initializing the Diagonal matrix.
00452     initDiagonal_ = params.get("Initialize Diagonal", false);
00453 
00454     // Set the number of blocks and allocate data
00455     setSize( rb );
00456   }
00457   
00459   // Set the block size and adjust as necessary
00460   template <class ScalarType, class MV, class OP>
00461   void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
00462   {
00463     // allocate space only; perform no computation
00464     // Any change in size invalidates the state of the solver as implemented here.
00465 
00466     TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
00467 
00468     if ( savedBlocks_ != savedBlocks) {
00469       stateStorageInitialized_ = false;
00470       savedBlocks_ = savedBlocks;
00471       initialized_ = false;
00472       curDim_ = 0;
00473       prevUdim_ = 0;
00474       setStateSize(); // Use the current savedBlocks_ to initialize the state storage.    
00475     }
00476   }
00477 
00479   // Enable the reuse of a single solver object for completely different linear systems
00480   template <class ScalarType, class MV, class OP>
00481   void PCPGIter<ScalarType,MV,OP>::resetState()
00482   {
00483       stateStorageInitialized_ = false;
00484       initialized_ = false;
00485       curDim_ = 0;
00486       prevUdim_ = 0;
00487       setStateSize();
00488   }
00489 
00491   // Setup the state storage.  Called by either initialize or, if savedBlocks_ changes, setSize.
00492   template <class ScalarType, class MV, class OP>
00493   void PCPGIter<ScalarType,MV,OP>::setStateSize ()
00494   {
00495     if (!stateStorageInitialized_) {
00496 
00497       // Check if there is any multivector to clone from.
00498       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00499       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00500       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00501   return;  // postpone exception 
00502       }
00503       else {
00504   
00506   // blockSize*recycledBlocks dependent
00507         int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
00508   //
00509   // Initialize the CG state storage
00510         // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
00511   // Generate CG state only if it does not exist, otherwise resize it.
00512         if (Z_ == Teuchos::null) {
00513           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00514           Z_ = MVT::Clone( *tmp, 1 );
00515         }
00516         if (P_ == Teuchos::null) {
00517           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00518           P_ = MVT::Clone( *tmp, 1 );
00519         }
00520         if (AP_ == Teuchos::null) {
00521           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00522           AP_ = MVT::Clone( *tmp, 1 );
00523         }
00524 
00525   if (C_ == Teuchos::null) {        
00526 
00527     // Get the multivector that is not null. 
00528     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00529     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00530            "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00531     TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00532            "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
00533     C_ = MVT::Clone( *tmp, savedBlocks_ );
00534   }
00535   else {
00536     // Generate C_ by cloning itself ONLY if more space is needed.
00537     if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
00538       Teuchos::RCP<const MV> tmp = C_;
00539       C_ = MVT::Clone( *tmp, savedBlocks_ );
00540     }
00541   }
00542   if (U_ == Teuchos::null) {        
00543     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00544     TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
00545            "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
00546     U_ = MVT::Clone( *tmp, savedBlocks_ );
00547   }
00548   else {
00549     // Generate U_ by cloning itself ONLY if more space is needed.
00550     if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
00551       Teuchos::RCP<const MV> tmp = U_;
00552       U_ = MVT::Clone( *tmp, savedBlocks_ );
00553     }
00554   }
00555         if (keepDiagonal_) {
00556           if (D_ == Teuchos::null) {
00557             D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00558           }
00559           if (initDiagonal_) {
00560             D_->shape( newsd, newsd );
00561           }
00562           else {
00563             if (D_->numRows() < newsd || D_->numCols() < newsd) {
00564               D_->shapeUninitialized( newsd, newsd );
00565             }
00566           }
00567         }
00568   // State storage has now been initialized.
00569   stateStorageInitialized_ = true;
00570       } // if there is a vector to clone from
00571     } // if !stateStorageInitialized_ 
00572   } // end of setStateSize
00573 
00575   // Initialize the iteration object
00576   template <class ScalarType, class MV, class OP>
00577   void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV> newstate)
00578   {
00579 
00580     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00581            "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
00582     
00583     // Requirements: R_ and consistent multivectors widths and lengths
00584     //
00585     std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
00586 
00587     if (newstate.R != Teuchos::null){ 
00588 
00589       R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
00590       if (newstate.U == Teuchos::null){ 
00591         prevUdim_ = 0;
00592         newstate.U = U_;
00593         newstate.C = C_;
00594       }
00595       else {
00596         prevUdim_ =  newstate.curDim;
00597         if (newstate.C == Teuchos::null){  // Stub for new feature
00598           std::vector<int> index(prevUdim_);
00599           for (int i=0; i< prevUdim_; ++i)  
00600             index[i] = i; 
00601           Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
00602           newstate.C = MVT::Clone( *newstate.U, prevUdim_ ); 
00603           Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );    
00604           lp_->apply( *Ukeff, *Ckeff );
00605         }
00606         curDim_ = prevUdim_ ;
00607       }
00608 
00609       // Initialize the state storage if not already allocated in the constructor
00610       if (!stateStorageInitialized_) 
00611         setStateSize();
00612 
00613       //TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_), std::invalid_argument, errstr );
00614       //TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
00615 
00616       newstate.prevUdim =  prevUdim_; // big change in functionality from GCRODR 
00617       newstate.curDim =  curDim_; 
00618 
00619       //TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
00620 
00621       std::vector<int> zero_index(1);
00622       zero_index[0] = 0;
00623       if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction 
00624         lp_->applyLeftPrec( *R_, *Z_ );
00625         MVT::SetBlock( *Z_,  zero_index , *P_ );  // P(:,zero_index) := Z
00626       } else {                                      
00627         Z_ = R_;
00628         MVT::SetBlock( *R_, zero_index, *P_ );
00629       }
00630 
00631       std::vector<int> nextind(1);
00632       nextind[0] = curDim_;
00633 
00634       MVT::SetBlock( *P_,  nextind, *newstate.U ); // U(:,curDim_ ) := P_
00635 
00636       ++curDim_;
00637       newstate.curDim = curDim_; 
00638 
00639       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
00640                           std::invalid_argument, errstr );
00641       if (newstate.U != U_) { // Why this is needed?
00642   U_ = newstate.U;
00643       }
00644 
00645       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
00646                           std::invalid_argument, errstr );
00647       if (newstate.C != C_) {
00648   C_ = newstate.C;
00649       }
00650     }
00651     else {
00652 
00653       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00654                          "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
00655     }
00656 
00657     // the solver is initialized
00658     initialized_ = true;
00659   }
00660 
00661 
00663   // Iterate until the status test informs us we should stop.
00664   template <class ScalarType, class MV, class OP>
00665   void PCPGIter<ScalarType,MV,OP>::iterate()
00666   {
00667     //
00668     // Allocate/initialize data structures
00669     //
00670     if (initialized_ == false) {
00671       initialize();
00672     }
00673     bool debug = false;
00674 
00675     // Allocate memory for scalars.
00676     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00677     Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
00678     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00679     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
00680 
00681     if( iter_ != 0 )
00682       std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl;  //DMD
00683 
00684     // GenOrtho Project Stubs
00685     std::vector<int> prevInd;
00686     Teuchos::RCP<const MV> Uprev;
00687     Teuchos::RCP<const MV> Cprev;
00688     Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
00689 
00690     if( prevUdim_ ){
00691       prevInd.resize( prevUdim_ );
00692       for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
00693       CZ.reshape( prevUdim_ , 1 );
00694       Uprev = MVT::CloneView(*U_, prevInd);
00695       Cprev = MVT::CloneView(*C_, prevInd);
00696     }
00697 
00698     // Get the current solution std::vector.
00699     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00700 
00701     // Check that the current solution std::vector only has one column.
00702     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
00703                         "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
00704 
00705     //Check that the input is correctly set up 
00706     TEST_FOR_EXCEPTION( curDim_  != prevUdim_ + 1, PCPGIterInitFailure,
00707                         "Belos::CGIter::iterate(): mistake in initialization !" );
00708 
00709 
00710     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00711     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00712 
00713 
00714     std::vector<int> curind(1);
00715     std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
00716     if (prevUdim_ > 0){                 // A-orthonalize P=Z to Uprev
00717       Teuchos::RCP<MV> P; 
00718       curind[0] = curDim_ - 1;          // column = dimension - 1 
00719       P = MVT::CloneViewNonConst(*U_,curind); 
00720       MVT::MvTransMv( one, *Cprev, *P, CZ );
00721       MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );       // P -= U*(C'Z)
00722 
00723       if( debug ){
00724         MVT::MvTransMv( one, *Cprev, *P, CZ );
00725         std::cout << " Input CZ post ortho " << std::endl;
00726         CZ.print( std::cout );
00727       }
00728       if( curDim_ == savedBlocks_ ){
00729         std::vector<int> zero_index(1);
00730         zero_index[0] = 0;
00731         MVT::SetBlock( *P, zero_index, *P_ );
00732       }
00733       P = Teuchos::null;
00734     }
00735 
00736     // Compute first <r,z> a.k.a. rHz
00737     MVT::MvTransMv( one, *R_, *Z_, rHz );
00738 
00740     // iterate until the status test is satisfied
00741     //
00742     while (stest_->checkStatus(this) != Passed ) {
00743       Teuchos::RCP<const MV> P; 
00744       Teuchos::RCP<MV> AP;
00745       iter_++;                          // The next iteration begins.
00746       //std::vector<int> curind(1);
00747       curind[0] = curDim_ - 1;          // column = dimension - 1 
00748       if( debug ){
00749         MVT::MvNorm(*R_, rnorm);
00750         std::cout << iter_ << "  " << curDim_ <<  "   " << rnorm[0] << std::endl;
00751       }
00752       if( prevUdim_ + iter_ < savedBlocks_ ){
00753         P = MVT::CloneView(*U_,curind); 
00754         AP = MVT::CloneViewNonConst(*C_,curind); 
00755         lp_->applyOp( *P, *AP );
00756         MVT::MvTransMv( one, *P, *AP, pAp );
00757       }else{
00758         if( prevUdim_ + iter_ == savedBlocks_ ){
00759           AP = MVT::CloneViewNonConst(*C_,curind); 
00760           lp_->applyOp( *P_, *AP );
00761           MVT::MvTransMv( one, *P_, *AP, pAp );
00762         }else{
00763           lp_->applyOp( *P_, *AP_ );
00764           MVT::MvTransMv( one, *P_, *AP_, pAp );
00765         }
00766       }
00767 
00768       if( keepDiagonal_  && prevUdim_ + iter_ <= savedBlocks_ )
00769         (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
00770 
00771       // positive pAp required 
00772       TEST_FOR_EXCEPTION( pAp(0,0) <= zero, PCPGIterateFailure,
00773                           "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00774 
00775       // alpha := <R_,Z_> / <P,AP>
00776       alpha(0,0) = rHz(0,0) / pAp(0,0);
00777 
00778       // positive alpha required 
00779       TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
00780                           "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
00781 
00782       // solution update  x += alpha * P
00783       if( curDim_ < savedBlocks_ ){
00784          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
00785       }else{
00786          MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00787       }
00788       //lp_->updateSolution(); ... does nothing.
00789       //
00790       // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
00791       //
00792       rHz_old(0,0) = rHz(0,0);
00793       //
00794       // residual update R_ := R_ - alpha * AP
00795       //
00796       if( prevUdim_ + iter_ <= savedBlocks_ ){
00797          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
00798          AP = Teuchos::null;
00799       }else{
00800          MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00801       }
00802       //
00803       // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
00804       //
00805       if ( lp_->getLeftPrec() != Teuchos::null ) {
00806         lp_->applyLeftPrec( *R_, *Z_ );
00807       } else {
00808         Z_ = R_;
00809       }
00810       //
00811       MVT::MvTransMv( one, *R_, *Z_, rHz );
00812       //
00813       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00814       //
00815       if( curDim_ < savedBlocks_ ){
00816          curDim_++;                                                         // update basis dim
00817          curind[0] = curDim_ - 1;
00818          Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
00819          MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
00820          if( prevUdim_ ){ // Deflate seed space 
00821              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00822              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
00823              if( debug ){
00824                std::cout << " Check CZ " << std::endl;
00825                MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
00826                CZ.print( std::cout );
00827              }
00828          }
00829          P = Teuchos::null;
00830          if( curDim_ == savedBlocks_ ){
00831            std::vector<int> zero_index(1);
00832            zero_index[0] = 0;
00833            MVT::SetBlock( *Pnext, zero_index, *P_ );
00834          }
00835          Pnext = Teuchos::null;
00836       }else{
00837          MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00838          if( prevUdim_ ){ // Deflate seed space
00839              MVT::MvTransMv( one, *Cprev, *Z_, CZ );
00840              MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );       // P_ -= U*(C'Z)
00841 
00842              if( debug ){
00843                std::cout << " Check CZ " << std::endl;
00844                MVT::MvTransMv( one, *Cprev, *P_, CZ );
00845                CZ.print( std::cout );
00846              }
00847          }
00848       }
00849       // CGB: 5/26/2010
00850       // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
00851       // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
00852       // 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
00853       // same for AP
00854       TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
00855     } // end coupled two-term recursion
00856     if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
00857   }
00858 
00859 } // end Belos namespace
00860 
00861 #endif /* BELOS_PCPG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:14 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3