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   MVT::MvAddMv( 1.0, *newV, 0.0, *newV, *lclV );
00550   
00551   // Done with local pointers
00552   lclV = Teuchos::null;
00553       }
00554     }
00555 
00556 
00557     // Check size of Z
00558     Z_.resize(numRHS_);
00559     for (int i=0; i<numRHS_; ++i) {
00560       // Create a std::vector if we need to.
00561       if (Z_[i] == Teuchos::null) {
00562   Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
00563       }
00564       if (Z_[i]->length() < numBlocks_+1) {
00565   Z_[i]->shapeUninitialized(numBlocks_+1, 1); 
00566       }
00567       
00568       // Check that the newstate std::vector is consistent.
00569       TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
00570       
00571       // Put data into Z_, make sure old information is not still hanging around.
00572       if (newstate.Z[i] != Z_[i]) {
00573   if (curDim_==0)
00574     Z_[i]->putScalar();
00575   
00576         Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
00577         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
00578         lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
00579         lclZ->assign(newZ);
00580   
00581         // Done with local pointers
00582   lclZ = Teuchos::null;
00583       }
00584     }
00585 
00586 
00587     // Check size of H
00588     H_.resize(numRHS_);
00589     for (int i=0; i<numRHS_; ++i) {
00590       // Create a matrix if we need to.
00591       if (H_[i] == Teuchos::null) {
00592   H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00593       }
00594       if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
00595   H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
00596       }
00597       
00598       // Put data into H_ if it exists, make sure old information is not still hanging around.
00599       if ((int)newstate.H.size() == numRHS_) {
00600   
00601   // Check that the newstate matrix is consistent.
00602   TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument, 
00603          "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
00604   
00605   if (newstate.H[i] != H_[i]) {
00606     // H_[i]->putScalar();
00607     
00608     Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
00609     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00610     lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
00611     lclH->assign(newH);
00612     
00613     // Done with local pointers
00614     lclH = Teuchos::null;
00615   }
00616       }
00617     }      
00618     
00620     // Reinitialize storage for least squares solve
00621     //
00622     cs_.resize(numRHS_);
00623     sn_.resize(numRHS_);
00624       
00625     // Copy over rotation angles if they exist
00626     if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
00627       for (int i=0; i<numRHS_; ++i) {
00628         if (cs_[i] != newstate.cs[i])
00629           cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
00630         if (sn_[i] != newstate.sn[i])
00631           sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
00632       }
00633     } 
00634       
00635     // Resize or create the vectors as necessary
00636     for (int i=0; i<numRHS_; ++i) {
00637       if (cs_[i] == Teuchos::null) 
00638         cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
00639       else
00640         cs_[i]->resize(numBlocks_+1);
00641       if (sn_[i] == Teuchos::null)
00642         sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
00643       else
00644         sn_[i]->resize(numBlocks_+1);
00645     }
00646 
00647     // the solver is initialized
00648     initialized_ = true;
00649       
00650       /*
00651   if (om_->isVerbosity( Debug ) ) {
00652       // Check almost everything here
00653       CheckList chk;
00654       chk.checkV = true;
00655       chk.checkArn = true;
00656       chk.checkAux = true;
00657       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00658     }
00659     */
00660 
00661   }
00662 
00663 
00665   // Iterate until the status test informs us we should stop.
00666   template <class ScalarType, class MV, class OP>
00667   void PseudoBlockGmresIter<ScalarType,MV,OP>::iterate()
00668   {
00669     //
00670     // Allocate/initialize data structures
00671     //
00672     if (initialized_ == false) {
00673       initialize();
00674     }
00675 
00676     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00677     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00678     
00679     // Compute the current search dimension. 
00680     int searchDim = numBlocks_;
00681     //
00682     // Associate each initial block of V_[i] with U_vec[i]
00683     // Reset the index std::vector (this might have been changed if there was a restart)
00684     //
00685     std::vector<int> index(1);
00686     std::vector<int> index2(1);
00687     index[0] = curDim_;
00688     Teuchos::RCP<MV> tmp_vec;
00689     Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
00690 
00691     // Create AU_vec to hold A*U_vec.
00692     Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
00693 
00694     for (int i=0; i<numRHS_; ++i) {
00695       index2[0] = i;
00696       tmp_vec = MVT::CloneView( *V_[i], index );
00697       MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *MVT::CloneView( *U_vec, index2 ) );
00698     }
00699     
00701     // iterate until the status test tells us to stop.
00702     //
00703     // also break if our basis is full
00704     //
00705     while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
00706 
00707       iter_++;
00708       //
00709       // Apply the operator to _work_vector
00710       //
00711       lp_->apply( *U_vec, *AU_vec );
00712       //
00713       //
00714       // Resize index.
00715       //
00716       int num_prev = curDim_+1;
00717       index.resize( num_prev );
00718       for (int i=0; i<num_prev; ++i) { 
00719   index[i] = i; 
00720       }
00721       //
00722       // Orthogonalize next Krylov std::vector for each right-hand side.
00723       //
00724       for (int i=0; i<numRHS_; ++i) {
00725   //
00726   // Get previous Krylov vectors.
00727   //
00728   RCP<MV> V_prev = MVT::CloneView( *V_[i], index );
00729   Teuchos::Array< RCP<const MV> > V_array( 1, V_prev );
00730   //
00731   // Get a view of the new candidate std::vector.
00732   //
00733   index2[0] = i;
00734   RCP<MV> V_new = MVT::CloneView( *AU_vec, index2 );
00735   //
00736   // Get a view of the current part of the upper-hessenberg matrix.
00737   //
00738   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new 
00739     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00740         ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
00741   Teuchos::Array< RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
00742   
00743   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
00744     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00745         ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
00746   //
00747   // Orthonormalize the new block of the Krylov expansion
00748   // NOTE:  Rank deficiencies are not checked because this is a single-std::vector Krylov method.
00749   //
00750   ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
00751   //
00752   // NOTE:  V_new is a copy of the iter+1 std::vector in V_[i], so the normalized std::vector has to be
00753   // be copied back in when V_new is changed.  
00754   //
00755   index2[0] = curDim_+1;
00756   tmp_vec = MVT::CloneView( *V_[i], index2 );
00757   MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
00758       }
00759       // 
00760       // Now _AU_vec is the new _U_vec, so swap these two vectors.
00761       // NOTE: This alleviates the need for allocating a std::vector for AU_vec each iteration.
00762       // 
00763       RCP<MV> tmp_AU_vec = U_vec;
00764       U_vec = AU_vec;
00765       AU_vec = tmp_AU_vec;
00766       //
00767       // V has been extended, and H has been extended. 
00768       //
00769       // Update the QR factorization of the upper Hessenberg matrix
00770       //
00771       updateLSQR();
00772       //
00773       // Update basis dim and release all pointers.
00774       //
00775       curDim_ += 1;
00776       //        
00777       /*      
00778       // When required, monitor some orthogonalities
00779       if (om_->isVerbosity( Debug ) ) {
00780       // Check almost everything here
00781       CheckList chk;
00782       chk.checkV = true;
00783       chk.checkArn = true;
00784       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00785       }
00786       else if (om_->isVerbosity( OrthoDetails ) ) {
00787         CheckList chk;
00788         chk.checkV = true;
00789         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00790       }
00791       */ 
00792       
00793     } // end while (statusTest == false)
00794    
00795   }
00796 
00798   // Update the least squares solution for each right-hand side.
00799   template<class ScalarType, class MV, class OP>
00800   void PseudoBlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00801   {
00802     // Get correct dimension based on input "dim"
00803     // Remember that ortho failures result in an exit before updateLSQR() is called.
00804     // Therefore, it is possible that dim == curDim_.
00805     int curDim = curDim_;
00806     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00807       curDim = dim;
00808     }
00809     
00810     int i, j;
00811     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00812 
00813     Teuchos::BLAS<int, ScalarType> blas;
00814     
00815     for (i=0; i<numRHS_; ++i) {
00816       //
00817       // Update the least-squares QR for each linear system.
00818       //
00819       // QR factorization of Least-Squares system with Givens rotations
00820       //
00821       for (j=0; j<curDim; j++) {
00822   //
00823   // Apply previous Givens rotations to new column of Hessenberg matrix
00824   //
00825   blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
00826       }
00827       //
00828       // Calculate new Givens rotation
00829       //
00830       blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00831       (*H_[i])(curDim+1,curDim) = zero;
00832       //
00833       // Update RHS w/ new transformation
00834       //
00835       blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00836     }
00837 
00838   } // end updateLSQR()
00839 
00840 } // end Belos namespace
00841 
00842 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */

Generated on Wed May 12 21:30:08 2010 for Belos by  doxygen 1.4.7