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