BelosPseudoBlockGmresIter.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_PSEUDO_BLOCK_GMRES_ITER_HPP
00030 #define BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 #include "BelosIteration.hpp"
00039 
00040 #include "BelosLinearProblem.hpp"
00041 #include "BelosMatOrthoManager.hpp"
00042 #include "BelosOutputManager.hpp"
00043 #include "BelosStatusTest.hpp"
00044 #include "BelosOperatorTraits.hpp"
00045 #include "BelosMultiVecTraits.hpp"
00046 
00047 #include "Teuchos_BLAS.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 
00067 namespace Belos {
00068   
00070 
00071   
00076   template <class ScalarType, class MV>
00077   struct PseudoBlockGmresIterState {
00078 
00079     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00080     typedef typename SCT::magnitudeType MagnitudeType;
00081 
00086     int curDim;
00088     std::vector<Teuchos::RCP<const MV> > V;
00094     std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > H;
00096     std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > R;
00098     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > Z;
00100     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > sn;
00101     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > > cs;
00102 
00103     PseudoBlockGmresIterState() : curDim(0), V(0),
00104           H(0), R(0), Z(0),
00105                                   sn(0), cs(0)
00106     {}
00107   };
00108   
00110 
00111   
00124   class PseudoBlockGmresIterInitFailure : public BelosError {public:
00125     PseudoBlockGmresIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00126     {}};
00127   
00134   class PseudoBlockGmresIterOrthoFailure : public BelosError {public:
00135     PseudoBlockGmresIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00136     {}};
00137   
00139   
00140   
00141   template<class ScalarType, class MV, class OP>
00142   class PseudoBlockGmresIter : virtual public Iteration<ScalarType,MV,OP> {
00143     
00144   public:
00145     
00146     //
00147     // Convenience typedefs
00148     //
00149     typedef MultiVecTraits<ScalarType,MV> MVT;
00150     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00151     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00152     typedef typename SCT::magnitudeType MagnitudeType;
00153     
00155 
00156     
00165     PseudoBlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00166         const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00167         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00168         const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00169         Teuchos::ParameterList &params );
00170     
00172     virtual ~PseudoBlockGmresIter() {};
00174     
00175     
00177 
00178     
00200     void iterate();
00201     
00223     void initialize(PseudoBlockGmresIterState<ScalarType,MV> newstate);
00224     
00228     void initialize()
00229     {
00230       PseudoBlockGmresIterState<ScalarType,MV> empty;
00231       initialize(empty);
00232     }
00233     
00241     PseudoBlockGmresIterState<ScalarType,MV> getState() const {
00242       PseudoBlockGmresIterState<ScalarType,MV> state;
00243       state.curDim = curDim_;
00244       state.V.resize(numRHS_);
00245       state.H.resize(numRHS_);
00246       state.Z.resize(numRHS_);
00247       state.sn.resize(numRHS_);
00248       state.cs.resize(numRHS_);  
00249       for (int i=0; i<numRHS_; ++i) {
00250   state.V[i] = V_[i];
00251   state.H[i] = H_[i];
00252   state.Z[i] = Z_[i];
00253         state.sn[i] = sn_[i];
00254         state.cs[i] = cs_[i];
00255       }
00256       return state;
00257     }
00258     
00260     
00261     
00263 
00264     
00266     int getNumIters() const { return iter_; }
00267     
00269     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00270     
00273     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00274     
00276 
00281     Teuchos::RCP<MV> getCurrentUpdate() const;
00282     
00284 
00287     void updateLSQR( int dim = -1 );
00288     
00290     int getCurSubspaceDim() const { 
00291       if (!initialized_) return 0;
00292       return curDim_;
00293     };
00294     
00296     int getMaxSubspaceDim() const { return numBlocks_; }
00297     
00299     
00300     
00302 
00303     
00305     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00306     
00308     int getBlockSize() const { return 1; }
00309     
00311     void setBlockSize(int blockSize) { 
00312       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00313        "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
00314     }
00315     
00317     int getNumBlocks() const { return numBlocks_; }
00318     
00320     void setNumBlocks(int numBlocks);
00321     
00323     bool isInitialized() { return initialized_; }
00324     
00326     
00327   private:
00328     
00329     //
00330     // Classes inputed through constructor that define the linear problem to be solved.
00331     //
00332     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00333     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00334     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00335     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00336     
00337     //
00338     // Algorithmic parameters
00339     //  
00340     // numRHS_ is the current number of linear systems being solved.
00341     int numRHS_;
00342     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00343     int numBlocks_; 
00344     
00345     // Storage for QR factorization of the least squares system.
00346     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
00347     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
00348     
00349     // Pointers to a work std::vector used to improve aggregate performance.
00350     RCP<MV> U_vec_, AU_vec_;    
00351 
00352     // Pointers to the current right-hand side and solution multivecs being solved for.
00353     RCP<MV> cur_block_rhs_, cur_block_sol_;
00354 
00355     // 
00356     // Current solver state
00357     //
00358     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00359     // is capable of running; _initialize is controlled  by the initialize() member method
00360     // For the implications of the state of initialized_, please see documentation for initialize()
00361     bool initialized_;
00362     
00363     // Current subspace dimension, and number of iterations performed.
00364     int curDim_, iter_;
00365 
00366     // 
00367     // State Storage
00368     //
00369     std::vector<Teuchos::RCP<MV> > V_;
00370     //
00371     // Projected matrices
00372     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00373     //
00374     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
00375     // 
00376     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00377     // R_: Upper triangular reduction of H
00378     // Z_: Q applied to right-hand side of the least squares system
00379     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
00380     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;  
00381   };
00382   
00384   // Constructor.
00385   template<class ScalarType, class MV, class OP>
00386   PseudoBlockGmresIter<ScalarType,MV,OP>::PseudoBlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00387                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00388                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00389                      const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00390                      Teuchos::ParameterList &params ):
00391     lp_(problem),
00392     om_(printer),
00393     stest_(tester),
00394     ortho_(ortho),
00395     numRHS_(0),
00396     numBlocks_(0),
00397     initialized_(false),
00398     curDim_(0),
00399     iter_(0)
00400   {
00401     // Get the maximum number of blocks allowed for each Krylov subspace
00402     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00403                        "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00404     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00405 
00406     setNumBlocks( nb );
00407   }
00408   
00410   // Set the block size and make necessary adjustments.
00411   template <class ScalarType, class MV, class OP>
00412   void PseudoBlockGmresIter<ScalarType,MV,OP>::setNumBlocks (int numBlocks)
00413   {
00414     // This routine only allocates space; it doesn't not perform any computation
00415     // any change in size will invalidate the state of the solver.
00416     
00417     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
00418 
00419     numBlocks_ = numBlocks;
00420     curDim_ = 0;
00421 
00422     initialized_ = false;
00423   }
00424 
00426   // Get the current update from this subspace.
00427   template <class ScalarType, class MV, class OP>
00428   Teuchos::RCP<MV> PseudoBlockGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00429   {
00430     //
00431     // If this is the first iteration of the Arnoldi factorization, 
00432     // there is no update, so return Teuchos::null. 
00433     //
00434     RCP<MV> currentUpdate = Teuchos::null;
00435     if (curDim_==0) { 
00436       return currentUpdate; 
00437     } else {
00438       currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
00439       std::vector<int> index(1), index2(curDim_);
00440       for (int i=0; i<curDim_; ++i) {
00441         index2[i] = i;
00442       }
00443       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00444       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00445       Teuchos::BLAS<int,ScalarType> blas;
00446       
00447       for (int i=0; i<numRHS_; ++i) {
00448         index[0] = i;
00449         RCP<MV> cur_block_copy_vec = MVT::CloneView( *currentUpdate, index );
00450         //
00451         //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00452         //
00453         Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
00454         //
00455         //  Solve the least squares problem and compute current solutions.
00456         //
00457         blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00458              Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00459        H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
00460   
00461   RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
00462   MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
00463       }
00464     }
00465     return currentUpdate;
00466   }
00467   
00468 
00470   // Get the native residuals stored in this iteration.  
00471   // Note:  No residual std::vector will be returned by Gmres.
00472   template <class ScalarType, class MV, class OP>
00473   Teuchos::RCP<const MV> PseudoBlockGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00474   {
00475     //
00476     // NOTE: Make sure the incoming std::vector is the correct size!
00477     //
00478     if ( norms && (int)norms->size() < numRHS_ )                         
00479       norms->resize( numRHS_ );                                          
00480 
00481     if (norms) {
00482       Teuchos::BLAS<int,ScalarType> blas;
00483       for (int j=0; j<numRHS_; j++) {
00484         (*norms)[j] = Teuchos::ScalarTraits<ScalarType>::magnitude( (*Z_[j])(curDim_) );
00485       }
00486     }
00487     return Teuchos::null;
00488   }
00489 
00490   
00491 
00493   // Initialize this iteration object
00494   template <class ScalarType, class MV, class OP>
00495   void PseudoBlockGmresIter<ScalarType,MV,OP>::initialize(PseudoBlockGmresIterState<ScalarType,MV> newstate)
00496   {
00497     // Get the number of right-hand sides we're solving for now.
00498     int numRHS = MVT::GetNumberVecs(*(lp_->getCurrLHSVec()));
00499     numRHS_ = numRHS;
00500 
00501     // NOTE:  In PseudoBlockGmresIter, V and Z are required!!!  
00502     // inconsitent multivectors widths and lengths will not be tolerated, and
00503     // will be treated with exceptions.
00504     //
00505     std::string errstr("Belos::PseudoBlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
00506 
00507     // Check that we have V and Z.
00508     TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0, std::invalid_argument,
00509                        "Belos::PseudoBlockGmresIter::initialize(): V and/or Z is not specified.");
00510 
00511     // Get the multivector that is not null.
00512     Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00513     Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00514     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00515     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00516                        "Belos::PseudoBlockGmresIter::initialize(): linear problem does not specify multivectors to clone from.");
00517 
00518     // Check the new dimension is not more that the maximum number of allowable blocks.  
00519     TEST_FOR_EXCEPTION( newstate.curDim > numBlocks_+1,
00520                         std::invalid_argument, errstr );
00521     curDim_ = newstate.curDim;
00522 
00523     // Initialize the state storage
00524     // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00525     V_.resize(numRHS_);
00526     for (int i=0; i<numRHS_; ++i) {
00527       // Create a new std::vector if we need to.
00528       if (V_[i] == Teuchos::null || MVT::GetNumberVecs(*V_[i]) < numBlocks_+1 ) {
00529         V_[i] = MVT::Clone(*tmp,numBlocks_+1);
00530       }
00531       // Check that the newstate std::vector is consistent.
00532       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V[i]) != MVT::GetVecLength(*V_[i]),
00533                           std::invalid_argument, errstr );
00534       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
00535                           std::invalid_argument, errstr );
00536       
00537       int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
00538       if (newstate.V[i] != V_[i]) {
00539         // Cnly copy over the first block and print a warning.
00540         if (curDim_ == 0 && lclDim > 1) {
00541           om_->stream(Warnings) << "Belos::PseudoBlockGmresIter::initialize(): the solver was initialized with a kernel of " 
00542                                 << lclDim << std::endl << "The block size however is only " << 1 << std::endl
00543                                 << "The last " << lclDim - 1 << " vectors will be discarded." << std::endl;
00544         }
00545         std::vector<int> nevind(curDim_+1);
00546         for (int j=0; j<curDim_+1; j++) nevind[j] = j;
00547         Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V[i], nevind );
00548         Teuchos::RCP<MV> lclV = MVT::CloneView( *V_[i], nevind );
00549         const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00550         const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00551         MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
00552 
00553         // Done with local pointers
00554         lclV = Teuchos::null;
00555       }
00556     }
00557 
00558 
00559     // Check size of Z
00560     Z_.resize(numRHS_);
00561     for (int i=0; i<numRHS_; ++i) {
00562       // Create a std::vector if we need to.
00563       if (Z_[i] == Teuchos::null) {
00564   Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
00565       }
00566       if (Z_[i]->length() < numBlocks_+1) {
00567   Z_[i]->shapeUninitialized(numBlocks_+1, 1); 
00568       }
00569       
00570       // Check that the newstate std::vector is consistent.
00571       TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
00572       
00573       // Put data into Z_, make sure old information is not still hanging around.
00574       if (newstate.Z[i] != Z_[i]) {
00575   if (curDim_==0)
00576     Z_[i]->putScalar();
00577   
00578         Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
00579         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
00580         lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
00581         lclZ->assign(newZ);
00582   
00583         // Done with local pointers
00584   lclZ = Teuchos::null;
00585       }
00586     }
00587 
00588 
00589     // Check size of H
00590     H_.resize(numRHS_);
00591     for (int i=0; i<numRHS_; ++i) {
00592       // Create a matrix if we need to.
00593       if (H_[i] == Teuchos::null) {
00594   H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00595       }
00596       if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
00597   H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
00598       }
00599       
00600       // Put data into H_ if it exists, make sure old information is not still hanging around.
00601       if ((int)newstate.H.size() == numRHS_) {
00602   
00603   // Check that the newstate matrix is consistent.
00604   TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument, 
00605          "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
00606   
00607   if (newstate.H[i] != H_[i]) {
00608     // H_[i]->putScalar();
00609     
00610     Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
00611     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00612     lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
00613     lclH->assign(newH);
00614     
00615     // Done with local pointers
00616     lclH = Teuchos::null;
00617   }
00618       }
00619     }      
00620     
00622     // Reinitialize storage for least squares solve
00623     //
00624     cs_.resize(numRHS_);
00625     sn_.resize(numRHS_);
00626       
00627     // Copy over rotation angles if they exist
00628     if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
00629       for (int i=0; i<numRHS_; ++i) {
00630         if (cs_[i] != newstate.cs[i])
00631           cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
00632         if (sn_[i] != newstate.sn[i])
00633           sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
00634       }
00635     } 
00636       
00637     // Resize or create the vectors as necessary
00638     for (int i=0; i<numRHS_; ++i) {
00639       if (cs_[i] == Teuchos::null) 
00640         cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
00641       else
00642         cs_[i]->resize(numBlocks_+1);
00643       if (sn_[i] == Teuchos::null)
00644         sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
00645       else
00646         sn_[i]->resize(numBlocks_+1);
00647     }
00648 
00649     // the solver is initialized
00650     initialized_ = true;
00651       
00652       /*
00653   if (om_->isVerbosity( Debug ) ) {
00654       // Check almost everything here
00655       CheckList chk;
00656       chk.checkV = true;
00657       chk.checkArn = true;
00658       chk.checkAux = true;
00659       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00660     }
00661     */
00662 
00663   }
00664 
00665 
00667   // Iterate until the status test informs us we should stop.
00668   template <class ScalarType, class MV, class OP>
00669   void PseudoBlockGmresIter<ScalarType,MV,OP>::iterate()
00670   {
00671     //
00672     // Allocate/initialize data structures
00673     //
00674     if (initialized_ == false) {
00675       initialize();
00676     }
00677 
00678     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00679     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00680     
00681     // Compute the current search dimension. 
00682     int searchDim = numBlocks_;
00683     //
00684     // Associate each initial block of V_[i] with U_vec[i]
00685     // Reset the index std::vector (this might have been changed if there was a restart)
00686     //
00687     std::vector<int> index(1);
00688     std::vector<int> index2(1);
00689     index[0] = curDim_;
00690     Teuchos::RCP<MV> tmp_vec;
00691     Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
00692 
00693     // Create AU_vec to hold A*U_vec.
00694     Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
00695 
00696     for (int i=0; i<numRHS_; ++i) {
00697       index2[0] = i;
00698       tmp_vec = MVT::CloneView( *V_[i], index );
00699       MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *MVT::CloneView( *U_vec, index2 ) );
00700     }
00701     
00703     // iterate until the status test tells us to stop.
00704     //
00705     // also break if our basis is full
00706     //
00707     while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
00708 
00709       iter_++;
00710       //
00711       // Apply the operator to _work_vector
00712       //
00713       lp_->apply( *U_vec, *AU_vec );
00714       //
00715       //
00716       // Resize index.
00717       //
00718       int num_prev = curDim_+1;
00719       index.resize( num_prev );
00720       for (int i=0; i<num_prev; ++i) { 
00721   index[i] = i; 
00722       }
00723       //
00724       // Orthogonalize next Krylov std::vector for each right-hand side.
00725       //
00726       for (int i=0; i<numRHS_; ++i) {
00727   //
00728   // Get previous Krylov vectors.
00729   //
00730   RCP<MV> V_prev = MVT::CloneView( *V_[i], index );
00731   Teuchos::Array< RCP<const MV> > V_array( 1, V_prev );
00732   //
00733   // Get a view of the new candidate std::vector.
00734   //
00735   index2[0] = i;
00736   RCP<MV> V_new = MVT::CloneView( *AU_vec, index2 );
00737   //
00738   // Get a view of the current part of the upper-hessenberg matrix.
00739   //
00740   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new 
00741     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00742         ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
00743   Teuchos::Array< RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
00744   
00745   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
00746     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00747         ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
00748   //
00749   // Orthonormalize the new block of the Krylov expansion
00750   // NOTE:  Rank deficiencies are not checked because this is a single-std::vector Krylov method.
00751   //
00752   ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
00753   //
00754   // NOTE:  V_new is a copy of the iter+1 std::vector in V_[i], so the normalized std::vector has to be
00755   // be copied back in when V_new is changed.  
00756   //
00757   index2[0] = curDim_+1;
00758   tmp_vec = MVT::CloneView( *V_[i], index2 );
00759   MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
00760       }
00761       // 
00762       // Now _AU_vec is the new _U_vec, so swap these two vectors.
00763       // NOTE: This alleviates the need for allocating a std::vector for AU_vec each iteration.
00764       // 
00765       RCP<MV> tmp_AU_vec = U_vec;
00766       U_vec = AU_vec;
00767       AU_vec = tmp_AU_vec;
00768       //
00769       // V has been extended, and H has been extended. 
00770       //
00771       // Update the QR factorization of the upper Hessenberg matrix
00772       //
00773       updateLSQR();
00774       //
00775       // Update basis dim and release all pointers.
00776       //
00777       curDim_ += 1;
00778       //        
00779       /*      
00780       // When required, monitor some orthogonalities
00781       if (om_->isVerbosity( Debug ) ) {
00782       // Check almost everything here
00783       CheckList chk;
00784       chk.checkV = true;
00785       chk.checkArn = true;
00786       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00787       }
00788       else if (om_->isVerbosity( OrthoDetails ) ) {
00789         CheckList chk;
00790         chk.checkV = true;
00791         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00792       }
00793       */ 
00794       
00795     } // end while (statusTest == false)
00796    
00797   }
00798 
00800   // Update the least squares solution for each right-hand side.
00801   template<class ScalarType, class MV, class OP>
00802   void PseudoBlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00803   {
00804     // Get correct dimension based on input "dim"
00805     // Remember that ortho failures result in an exit before updateLSQR() is called.
00806     // Therefore, it is possible that dim == curDim_.
00807     int curDim = curDim_;
00808     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00809       curDim = dim;
00810     }
00811     
00812     int i, j;
00813     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00814 
00815     Teuchos::BLAS<int, ScalarType> blas;
00816     
00817     for (i=0; i<numRHS_; ++i) {
00818       //
00819       // Update the least-squares QR for each linear system.
00820       //
00821       // QR factorization of Least-Squares system with Givens rotations
00822       //
00823       for (j=0; j<curDim; j++) {
00824   //
00825   // Apply previous Givens rotations to new column of Hessenberg matrix
00826   //
00827   blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
00828       }
00829       //
00830       // Calculate new Givens rotation
00831       //
00832       blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00833       (*H_[i])(curDim+1,curDim) = zero;
00834       //
00835       // Update RHS w/ new transformation
00836       //
00837       blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00838     }
00839 
00840   } // end updateLSQR()
00841 
00842 } // end Belos namespace
00843 
00844 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7