BelosBlockGmresIter.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_BLOCK_GMRES_ITER_HPP
00030 #define BELOS_BLOCK_GMRES_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 #include "BelosGmresIteration.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   
00068 template<class ScalarType, class MV, class OP>
00069 class BlockGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
00070 
00071   public:
00072     
00073   //
00074   // Convenience typedefs
00075   //
00076   typedef MultiVecTraits<ScalarType,MV> MVT;
00077   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00078   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00079   typedef typename SCT::magnitudeType MagnitudeType;
00080 
00082 
00083 
00093   BlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00094       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00095       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00096       const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00097       Teuchos::ParameterList &params );
00098 
00100   virtual ~BlockGmresIter() {};
00102 
00103 
00105 
00106   
00128   void iterate();
00129 
00151   void initializeGmres(GmresIterationState<ScalarType,MV> newstate);
00152 
00156   void initialize()
00157   {
00158     GmresIterationState<ScalarType,MV> empty;
00159     initializeGmres(empty);
00160   }
00161   
00169   GmresIterationState<ScalarType,MV> getState() const {
00170     GmresIterationState<ScalarType,MV> state;
00171     state.curDim = curDim_;
00172     state.V = V_;
00173     state.H = H_;
00174     state.R = R_;
00175     state.z = z_;
00176     return state;
00177   }
00178 
00180 
00181   
00183 
00184 
00186   int getNumIters() const { return iter_; }
00187   
00189   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00190 
00193   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00194 
00196 
00201   Teuchos::RCP<MV> getCurrentUpdate() const;
00202 
00204 
00207   void updateLSQR( int dim = -1 );
00208 
00210   int getCurSubspaceDim() const { 
00211     if (!initialized_) return 0;
00212     return curDim_;
00213   };
00214 
00216   int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }
00217 
00219 
00220   
00222 
00223 
00225   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00226 
00228   int getBlockSize() const { return blockSize_; }
00229   
00231   void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }
00232   
00234   int getNumBlocks() const { return numBlocks_; }
00235   
00237   void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }
00238   
00245   void setSize(int blockSize, int numBlocks);
00246 
00248   bool isInitialized() { return initialized_; }
00249 
00251 
00252   private:
00253 
00254   //
00255   // Internal methods
00256   //
00258   void setStateSize();
00259   
00260   //
00261   // Classes inputed through constructor that define the linear problem to be solved.
00262   //
00263   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00264   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00265   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00266   const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00267 
00268   //
00269   // Algorithmic parameters
00270   //  
00271   // blockSize_ is the solver block size.
00272   // It controls the number of vectors added to the basis on each iteration.
00273   int blockSize_;
00274   // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00275   int numBlocks_; 
00276   
00277   // Storage for QR factorization of the least squares system.
00278   Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
00279   Teuchos::SerialDenseVector<int,MagnitudeType> cs;
00280   
00281   // 
00282   // Current solver state
00283   //
00284   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00285   // is capable of running; _initialize is controlled  by the initialize() member method
00286   // For the implications of the state of initialized_, please see documentation for initialize()
00287   bool initialized_;
00288 
00289   // stateStorageInitialized_ specifies that the state storage has be initialized to the current
00290   // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00291   // generated without the right-hand side or solution vectors.
00292   bool stateStorageInitialized_;
00293 
00294   // keepHessenberg_ specifies that the iteration must keep the Hessenberg matrix formed via the 
00295   // Arnoldi factorization and the upper triangular matrix that is the Hessenberg matrix reduced via
00296   // QR factorization separate.
00297   bool keepHessenberg_;
00298 
00299   // initHessenberg_ specifies that the iteration should reinitialize the Hessenberg matrix by zeroing
00300   // out all entries before an iteration is started.
00301   bool initHessenberg_;
00302  
00303   // Current subspace dimension, and number of iterations performed.
00304   int curDim_, iter_;
00305   
00306   // 
00307   // State Storage
00308   //
00309   Teuchos::RCP<MV> V_;
00310   //
00311   // Projected matrices
00312   // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00313   //
00314   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00315   // 
00316   // QR decomposition of Projected matrices for solving the least squares system HY = B.
00317   // R_: Upper triangular reduction of H
00318   // z_: Q applied to right-hand side of the least squares system
00319   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00320   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;  
00321 };
00322 
00324   // Constructor.
00325   template<class ScalarType, class MV, class OP>
00326   BlockGmresIter<ScalarType,MV,OP>::BlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00327                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00328                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00329                const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00330                Teuchos::ParameterList &params ):
00331     lp_(problem),
00332     om_(printer),
00333     stest_(tester),
00334     ortho_(ortho),
00335     blockSize_(0),
00336     numBlocks_(0),
00337     initialized_(false),
00338     stateStorageInitialized_(false),
00339     keepHessenberg_(false),
00340     initHessenberg_(false),
00341     curDim_(0),
00342     iter_(0)
00343   {
00344     // Find out whether we are saving the Hessenberg matrix.
00345     keepHessenberg_ = params.get("Keep Hessenberg", false);
00346 
00347     // Find out whether we are initializing the Hessenberg matrix.
00348     initHessenberg_ = params.get("Initialize Hessenberg", false);
00349 
00350     // Get the maximum number of blocks allowed for this Krylov subspace
00351     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00352                        "Belos::BlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00353     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00354 
00355     // Set the block size and allocate data
00356     int bs = params.get("Block Size", 1);
00357     setSize( bs, nb );
00358   }
00359 
00361   // Set the block size and make necessary adjustments.
00362   template <class ScalarType, class MV, class OP>
00363   void BlockGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
00364   {
00365     // This routine only allocates space; it doesn't not perform any computation
00366     // any change in size will invalidate the state of the solver.
00367 
00368     TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setSize was passed a non-positive argument.");
00369     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00370       // do nothing
00371       return;
00372     }
00373 
00374     if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
00375       stateStorageInitialized_ = false;
00376 
00377     blockSize_ = blockSize;
00378     numBlocks_ = numBlocks;
00379 
00380     initialized_ = false;
00381     curDim_ = 0;
00382 
00383     // Use the current blockSize_ and numBlocks_ to initialize the state storage.    
00384     setStateSize();
00385     
00386   }
00387 
00389   // Setup the state storage.
00390   template <class ScalarType, class MV, class OP>
00391   void BlockGmresIter<ScalarType,MV,OP>::setStateSize ()
00392   {
00393     if (!stateStorageInitialized_) {
00394 
00395       // Check if there is any multivector to clone from.
00396       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00397       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00398       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00399   stateStorageInitialized_ = false;
00400   return;
00401       }
00402       else {
00403   
00405   // blockSize*numBlocks dependent
00406   //
00407   int newsd = blockSize_*(numBlocks_+1);
00408   
00409   if (blockSize_==1) {
00410     cs.resize( newsd );
00411     sn.resize( newsd );
00412   }
00413   else {
00414     beta.resize( newsd );
00415   }
00416   
00417   // Initialize the state storage
00418         TEST_FOR_EXCEPTION(blockSize_*numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
00419                            "Belos::BlockGmresIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
00420 
00421   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00422   if (V_ == Teuchos::null) {
00423     // Get the multivector that is not null.
00424     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00425     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00426            "Belos::BlockGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00427     V_ = MVT::Clone( *tmp, newsd );
00428   }
00429   else {
00430     // Generate V_ by cloning itself ONLY if more space is needed.
00431     if (MVT::GetNumberVecs(*V_) < newsd) {
00432       Teuchos::RCP<const MV> tmp = V_;
00433       V_ = MVT::Clone( *tmp, newsd );
00434     }
00435   }
00436   
00437   // Generate R_ only if it doesn't exist, otherwise resize it.
00438   if (R_ == Teuchos::null) {
00439     R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00440         }
00441         if (initHessenberg_) {
00442           R_->shape( newsd, newsd-blockSize_ );
00443         }
00444         else {
00445           if (R_->numRows() < newsd || R_->numCols() < newsd-blockSize_) {
00446       R_->shapeUninitialized( newsd, newsd-blockSize_ );
00447     }
00448         }
00449   
00450   // Generate H_ only if it doesn't exist, and we are keeping the upper Hessenberg matrix.
00451   if (keepHessenberg_) {
00452     if (H_ == Teuchos::null) {
00453       H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00454           }
00455           if (initHessenberg_) {
00456             H_->shape( newsd, newsd-blockSize_ );
00457           }
00458           else {
00459             if (H_->numRows() < newsd || H_->numCols() < newsd-blockSize_) {
00460         H_->shapeUninitialized( newsd, newsd-blockSize_ );
00461       }
00462           }
00463   }
00464   else {
00465     // Point H_ and R_ at the same object.
00466     H_ = R_;
00467   }
00468   
00469   // Generate z_ only if it doesn't exist, otherwise resize it.
00470   if (z_ == Teuchos::null) {
00471     z_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00472         }
00473         if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
00474     z_->shapeUninitialized( newsd, blockSize_ );
00475   }
00476   
00477   // State storage has now been initialized.
00478   stateStorageInitialized_ = true;
00479       }
00480     }
00481   }
00482 
00484   // Get the current update from this subspace.
00485   template <class ScalarType, class MV, class OP>
00486   Teuchos::RCP<MV> BlockGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00487   {
00488     //
00489     // If this is the first iteration of the Arnoldi factorization, 
00490     // there is no update, so return Teuchos::null. 
00491     //
00492     RCP<MV> currentUpdate = Teuchos::null;
00493     if (curDim_==0) { 
00494       return currentUpdate; 
00495     } else {
00496       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00497       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00498       Teuchos::BLAS<int,ScalarType> blas;
00499       currentUpdate = MVT::Clone( *V_, blockSize_ );
00500       //
00501       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00502       //
00503       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
00504       //
00505       //  Solve the least squares problem.
00506       //
00507       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00508      Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,  
00509      R_->values(), R_->stride(), y.values(), y.stride() );
00510       //
00511       //  Compute the current update.
00512       //
00513       std::vector<int> index(curDim_);
00514       for ( int i=0; i<curDim_; i++ ) {   
00515         index[i] = i;
00516       }
00517       RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
00518       MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
00519     }
00520     return currentUpdate;
00521   }
00522 
00523 
00525   // Get the native residuals stored in this iteration.  
00526   // Note:  No residual std::vector will be returned by Gmres.
00527   template <class ScalarType, class MV, class OP>
00528   Teuchos::RCP<const MV> BlockGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00529   {
00530     //
00531     // NOTE: Make sure the incoming std::vector is the correct size!
00532     //
00533     if ( norms && (int)norms->size() < blockSize_ )                         
00534       norms->resize( blockSize_ );                                          
00535     
00536     if (norms) {
00537       Teuchos::BLAS<int,ScalarType> blas;
00538       for (int j=0; j<blockSize_; j++) {
00539         (*norms)[j] = blas.NRM2( blockSize_, &(*z_)(curDim_, j), 1);
00540       }
00541     }
00542     return Teuchos::null;
00543   }
00544 
00545   
00546 
00548   // Initialize this iteration object
00549   template <class ScalarType, class MV, class OP>
00550   void BlockGmresIter<ScalarType,MV,OP>::initializeGmres(GmresIterationState<ScalarType,MV> newstate)
00551   {
00552     // Initialize the state storage if it isn't already.
00553     if (!stateStorageInitialized_) 
00554       setStateSize();
00555 
00556     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00557            "Belos::BlockGmresIter::initialize(): Cannot initialize state storage!");
00558     
00559     // NOTE:  In BlockGmresIter, V and Z are required!!!  
00560     // inconsitent multivectors widths and lengths will not be tolerated, and
00561     // will be treated with exceptions.
00562     //
00563     std::string errstr("Belos::BlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
00564 
00565     if (newstate.V != Teuchos::null && newstate.z != Teuchos::null) {
00566 
00567       // initialize V_,z_, and curDim_
00568 
00569       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00570                           std::invalid_argument, errstr );
00571       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00572                           std::invalid_argument, errstr );
00573       TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*(numBlocks_+1),
00574                           std::invalid_argument, errstr );
00575 
00576       curDim_ = newstate.curDim;
00577       int lclDim = MVT::GetNumberVecs(*newstate.V);
00578 
00579       // check size of Z
00580       TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_, std::invalid_argument, errstr);
00581       
00582 
00583       // copy basis vectors from newstate into V
00584       if (newstate.V != V_) {
00585         // only copy over the first block and print a warning.
00586         if (curDim_ == 0 && lclDim > blockSize_) {
00587     om_->stream(Warnings) << "Belos::BlockGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00588         << "The block size however is only " << blockSize_ << std::endl
00589         << "The last " << lclDim - blockSize_ << " vectors will be discarded." << std::endl;
00590   }
00591         std::vector<int> nevind(curDim_+blockSize_);
00592         for (int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
00593   Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
00594   Teuchos::RCP<MV> lclV = MVT::CloneView( *V_, nevind );
00595         MVT::MvAddMv( 1.0, *newV, 0.0, *newV, *lclV );
00596 
00597         // done with local pointers
00598         lclV = Teuchos::null;
00599       }
00600 
00601       // put data into z_, make sure old information is not still hanging around.
00602       if (newstate.z != z_) {
00603         z_->putScalar();
00604         Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+blockSize_,blockSize_);
00605         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclZ;
00606         lclZ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
00607         lclZ->assign(newZ);
00608 
00609         // done with local pointers
00610         lclZ = Teuchos::null;
00611       }
00612 
00613     }
00614     else {
00615 
00616       TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
00617                          "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
00618 
00619       TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00620                          "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial norms z_0.");
00621     }
00622 
00623     // the solver is initialized
00624     initialized_ = true;
00625 
00626     /*
00627     if (om_->isVerbosity( Debug ) ) {
00628       // Check almost everything here
00629       CheckList chk;
00630       chk.checkV = true;
00631       chk.checkArn = true;
00632       chk.checkAux = true;
00633       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00634     }
00635     */
00636 
00637   }
00638 
00639 
00641   // Iterate until the status test informs us we should stop.
00642   template <class ScalarType, class MV, class OP>
00643   void BlockGmresIter<ScalarType,MV,OP>::iterate()
00644   {
00645     //
00646     // Allocate/initialize data structures
00647     //
00648     if (initialized_ == false) {
00649       initialize();
00650     }
00651     
00652     // Compute the current search dimension. 
00653     int searchDim = blockSize_*numBlocks_;
00654 
00656     // iterate until the status test tells us to stop.
00657     //
00658     // also break if our basis is full
00659     //
00660     while (stest_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
00661 
00662       iter_++;
00663 
00664       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
00665       int lclDim = curDim_ + blockSize_; 
00666 
00667       // Get the current part of the basis.
00668       std::vector<int> curind(blockSize_);
00669       for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
00670       Teuchos::RCP<MV> Vnext = MVT::CloneView(*V_,curind);
00671 
00672       // Get a view of the previous vectors.
00673       // This is used for orthogonalization and for computing V^H K H
00674       for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
00675       Teuchos::RCP<MV> Vprev = MVT::CloneView(*V_,curind);
00676 
00677       // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
00678       lp_->apply(*Vprev,*Vnext);
00679       Vprev = Teuchos::null;
00680       
00681       // Remove all previous Krylov basis vectors from Vnext      
00682       // Get a view of all the previous vectors
00683       std::vector<int> prevind(lclDim);
00684       for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00685       Vprev = MVT::CloneView(*V_,prevind);
00686       Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00687       
00688       // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00689       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00690   subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00691            ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
00692       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00693       AsubH.append( subH );
00694       
00695       // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00696       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00697   subH2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00698             ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
00699       subH2->putScalar();  // Initialize subdiagonal to zero
00700       int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
00701 
00702       // Copy over the coefficients if we are saving the upper Hessenberg matrix,
00703       // just in case we run into an error.
00704       if (keepHessenberg_) {
00705   // Copy over the orthogonalization coefficients.
00706   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00707     subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00708              ( Teuchos::View,*R_,lclDim,blockSize_,0,curDim_ ) );
00709   subR->assign(*subH);
00710   
00711   // Copy over the lower diagonal block of the Hessenberg matrix.
00712   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00713     subR2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00714         ( Teuchos::View,*R_,blockSize_,blockSize_,lclDim,curDim_ ) );
00715   subR2->assign(*subH2);
00716       }
00717 
00718       TEST_FOR_EXCEPTION(rank != blockSize_,GmresIterationOrthoFailure,
00719        "Belos::BlockGmresIter::iterate(): couldn't generate basis of full rank.");
00720       //
00721       // V has been extended, and H has been extended. 
00722       //
00723       // Update the QR factorization of the upper Hessenberg matrix
00724       //
00725       updateLSQR();
00726       //
00727       // Update basis dim and release all pointers.
00728       //
00729       Vnext = Teuchos::null;
00730       curDim_ += blockSize_;
00731       //        
00732       /*      
00733       // When required, monitor some orthogonalities
00734       if (om_->isVerbosity( Debug ) ) {
00735       // Check almost everything here
00736       CheckList chk;
00737       chk.checkV = true;
00738       chk.checkArn = true;
00739       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00740       }
00741       else if (om_->isVerbosity( OrthoDetails ) ) {
00742         CheckList chk;
00743         chk.checkV = true;
00744         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00745       }
00746       */ 
00747       
00748     } // end while (statusTest == false)
00749    
00750   }
00751 
00752   
00753   template<class ScalarType, class MV, class OP>
00754   void BlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00755   {
00756     int i, j, maxidx;
00757     ScalarType sigma, mu, vscale, maxelem;
00758     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00759     
00760     // Get correct dimension based on input "dim"
00761     // Remember that ortho failures result in an exit before updateLSQR() is called.
00762     // Therefore, it is possible that dim == curDim_.
00763     int curDim = curDim_;
00764     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00765       curDim = dim;
00766     }
00767     
00768     Teuchos::LAPACK<int, ScalarType> lapack;
00769     Teuchos::BLAS<int, ScalarType> blas;
00770     //
00771     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00772     // system to upper-triangular form.
00773     //
00774     if (blockSize_ == 1) {
00775       //
00776       // QR factorization of Least-Squares system with Givens rotations
00777       //
00778       for (i=0; i<curDim; i++) {
00779   //
00780   // Apply previous Givens rotations to new column of Hessenberg matrix
00781   //
00782   blas.ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
00783       }
00784       //
00785       // Calculate new Givens rotation
00786       //
00787       blas.ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
00788       (*R_)(curDim+1,curDim) = zero;
00789       //
00790       // Update RHS w/ new transformation
00791       //
00792       blas.ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
00793     }
00794     else {
00795       //
00796       // QR factorization of Least-Squares system with Householder reflectors
00797       //
00798       for (j=0; j<blockSize_; j++) {
00799   //
00800   // Apply previous Householder reflectors to new block of Hessenberg matrix
00801   //
00802   for (i=0; i<curDim+j; i++) {
00803     sigma = blas.DOT( blockSize_, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
00804     sigma += (*R_)(i,curDim+j);
00805     sigma *= beta[i];
00806     blas.AXPY(blockSize_, -sigma, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
00807     (*R_)(i,curDim+j) -= sigma;
00808   }
00809   //
00810   // Compute new Householder reflector
00811   //
00812   maxidx = blas.IAMAX( blockSize_+1, &(*R_)(curDim+j,curDim+j), 1 );
00813   maxelem = (*R_)(curDim+j+maxidx-1,curDim+j);
00814   for (i=0; i<blockSize_+1; i++)
00815     (*R_)(curDim+j+i,curDim+j) /= maxelem;
00816   sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j), 1,
00817         &(*R_)(curDim+j+1,curDim+j), 1 );
00818   if (sigma == zero) {
00819     beta[curDim + j] = zero;
00820   } else {
00821     mu = std::sqrt((*R_)(curDim+j,curDim+j)*(*R_)(curDim+j,curDim+j)+sigma);
00822     if ( Teuchos::ScalarTraits<ScalarType>::real((*R_)(curDim+j,curDim+j)) 
00823          < Teuchos::ScalarTraits<MagnitudeType>::zero() ) {
00824       vscale = (*R_)(curDim+j,curDim+j) - mu;
00825     } else {
00826       vscale = -sigma / ((*R_)(curDim+j,curDim+j) + mu);
00827     }
00828     beta[curDim+j] = 2.0*vscale*vscale/(sigma + vscale*vscale);
00829     (*R_)(curDim+j,curDim+j) = maxelem*mu;
00830     for (i=0; i<blockSize_; i++)
00831       (*R_)(curDim+j+1+i,curDim+j) /= vscale;
00832   }
00833   //
00834   // Apply new Householder reflector to rhs
00835   //
00836   for (i=0; i<blockSize_; i++) {
00837     sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j),
00838           1, &(*z_)(curDim+j+1,i), 1);
00839     sigma += (*z_)(curDim+j,i);
00840     sigma *= beta[curDim+j];
00841     blas.AXPY(blockSize_, -sigma, &(*R_)(curDim+j+1,curDim+j),
00842         1, &(*z_)(curDim+j+1,i), 1);
00843     (*z_)(curDim+j,i) -= sigma;
00844   }
00845       }
00846     } // end if (blockSize_ == 1)
00847 
00848     // If the least-squares problem is updated wrt "dim" then update the curDim_.
00849     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00850       curDim_ = dim + blockSize_;
00851     }
00852 
00853   } // end updateLSQR()
00854 
00855 } // end Belos namespace
00856 
00857 #endif /* BELOS_BLOCK_GMRES_ITER_HPP */

Generated on Wed May 12 21:45:50 2010 for Belos by  doxygen 1.4.7