Belos Package Browser (Single Doxygen Collection) Development
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_BLOCK_GMRES_ITER_HPP
00043 #define BELOS_BLOCK_GMRES_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosGmresIteration.hpp"
00052 
00053 #include "BelosLinearProblem.hpp"
00054 #include "BelosMatOrthoManager.hpp"
00055 #include "BelosOutputManager.hpp"
00056 #include "BelosStatusTest.hpp"
00057 #include "BelosOperatorTraits.hpp"
00058 #include "BelosMultiVecTraits.hpp"
00059 
00060 #include "Teuchos_BLAS.hpp"
00061 #include "Teuchos_SerialDenseMatrix.hpp"
00062 #include "Teuchos_SerialDenseVector.hpp"
00063 #include "Teuchos_ScalarTraits.hpp"
00064 #include "Teuchos_ParameterList.hpp"
00065 #include "Teuchos_TimeMonitor.hpp"
00066 
00080 namespace Belos {
00081   
00082 template<class ScalarType, class MV, class OP>
00083 class BlockGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
00084 
00085   public:
00086     
00087   //
00088   // Convenience typedefs
00089   //
00090   typedef MultiVecTraits<ScalarType,MV> MVT;
00091   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00092   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00093   typedef typename SCT::magnitudeType MagnitudeType;
00094 
00096 
00097 
00107   BlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00108       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00109       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00110       const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00111       Teuchos::ParameterList &params );
00112 
00114   virtual ~BlockGmresIter() {};
00116 
00117 
00119 
00120   
00142   void iterate();
00143 
00165   void initializeGmres(GmresIterationState<ScalarType,MV> newstate);
00166 
00170   void initialize()
00171   {
00172     GmresIterationState<ScalarType,MV> empty;
00173     initializeGmres(empty);
00174   }
00175   
00183   GmresIterationState<ScalarType,MV> getState() const {
00184     GmresIterationState<ScalarType,MV> state;
00185     state.curDim = curDim_;
00186     state.V = V_;
00187     state.H = H_;
00188     state.R = R_;
00189     state.z = z_;
00190     return state;
00191   }
00192 
00194 
00195   
00197 
00198 
00200   int getNumIters() const { return iter_; }
00201   
00203   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00204 
00207   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00208 
00210 
00215   Teuchos::RCP<MV> getCurrentUpdate() const;
00216 
00218 
00221   void updateLSQR( int dim = -1 );
00222 
00224   int getCurSubspaceDim() const { 
00225     if (!initialized_) return 0;
00226     return curDim_;
00227   };
00228 
00230   int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }
00231 
00233 
00234   
00236 
00237 
00239   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00240 
00242   int getBlockSize() const { return blockSize_; }
00243   
00245   void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }
00246   
00248   int getNumBlocks() const { return numBlocks_; }
00249   
00251   void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }
00252   
00259   void setSize(int blockSize, int numBlocks);
00260 
00262   bool isInitialized() { return initialized_; }
00263 
00265 
00266   private:
00267 
00268   //
00269   // Internal methods
00270   //
00272   void setStateSize();
00273   
00274   //
00275   // Classes inputed through constructor that define the linear problem to be solved.
00276   //
00277   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00278   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00279   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00280   const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00281 
00282   //
00283   // Algorithmic parameters
00284   //  
00285   // blockSize_ is the solver block size.
00286   // It controls the number of vectors added to the basis on each iteration.
00287   int blockSize_;
00288   // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00289   int numBlocks_; 
00290   
00291   // Storage for QR factorization of the least squares system.
00292   Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
00293   Teuchos::SerialDenseVector<int,MagnitudeType> cs;
00294   
00295   // 
00296   // Current solver state
00297   //
00298   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00299   // is capable of running; _initialize is controlled  by the initialize() member method
00300   // For the implications of the state of initialized_, please see documentation for initialize()
00301   bool initialized_;
00302 
00303   // stateStorageInitialized_ specifies that the state storage has be initialized to the current
00304   // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00305   // generated without the right-hand side or solution vectors.
00306   bool stateStorageInitialized_;
00307 
00308   // keepHessenberg_ specifies that the iteration must keep the Hessenberg matrix formed via the 
00309   // Arnoldi factorization and the upper triangular matrix that is the Hessenberg matrix reduced via
00310   // QR factorization separate.
00311   bool keepHessenberg_;
00312 
00313   // initHessenberg_ specifies that the iteration should reinitialize the Hessenberg matrix by zeroing
00314   // out all entries before an iteration is started.
00315   bool initHessenberg_;
00316  
00317   // Current subspace dimension, and number of iterations performed.
00318   int curDim_, iter_;
00319   
00320   // 
00321   // State Storage
00322   //
00323   Teuchos::RCP<MV> V_;
00324   //
00325   // Projected matrices
00326   // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00327   //
00328   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00329   // 
00330   // QR decomposition of Projected matrices for solving the least squares system HY = B.
00331   // R_: Upper triangular reduction of H
00332   // z_: Q applied to right-hand side of the least squares system
00333   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00334   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;  
00335 };
00336 
00338   // Constructor.
00339   template<class ScalarType, class MV, class OP>
00340   BlockGmresIter<ScalarType,MV,OP>::BlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00341                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00342                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00343                const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00344                Teuchos::ParameterList &params ):
00345     lp_(problem),
00346     om_(printer),
00347     stest_(tester),
00348     ortho_(ortho),
00349     blockSize_(0),
00350     numBlocks_(0),
00351     initialized_(false),
00352     stateStorageInitialized_(false),
00353     keepHessenberg_(false),
00354     initHessenberg_(false),
00355     curDim_(0),
00356     iter_(0)
00357   {
00358     // Find out whether we are saving the Hessenberg matrix.
00359     keepHessenberg_ = params.get("Keep Hessenberg", false);
00360 
00361     // Find out whether we are initializing the Hessenberg matrix.
00362     initHessenberg_ = params.get("Initialize Hessenberg", false);
00363 
00364     // Get the maximum number of blocks allowed for this Krylov subspace
00365     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00366                        "Belos::BlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00367     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00368 
00369     // Set the block size and allocate data
00370     int bs = params.get("Block Size", 1);
00371     setSize( bs, nb );
00372   }
00373 
00375   // Set the block size and make necessary adjustments.
00376   template <class ScalarType, class MV, class OP>
00377   void BlockGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
00378   {
00379     // This routine only allocates space; it doesn't not perform any computation
00380     // any change in size will invalidate the state of the solver.
00381 
00382     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setSize was passed a non-positive argument.");
00383     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00384       // do nothing
00385       return;
00386     }
00387 
00388     if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
00389       stateStorageInitialized_ = false;
00390 
00391     blockSize_ = blockSize;
00392     numBlocks_ = numBlocks;
00393 
00394     initialized_ = false;
00395     curDim_ = 0;
00396 
00397     // Use the current blockSize_ and numBlocks_ to initialize the state storage.    
00398     setStateSize();
00399     
00400   }
00401 
00403   // Setup the state storage.
00404   template <class ScalarType, class MV, class OP>
00405   void BlockGmresIter<ScalarType,MV,OP>::setStateSize ()
00406   {
00407     if (!stateStorageInitialized_) {
00408 
00409       // Check if there is any multivector to clone from.
00410       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00411       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00412       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00413   stateStorageInitialized_ = false;
00414   return;
00415       }
00416       else {
00417   
00419   // blockSize*numBlocks dependent
00420   //
00421   int newsd = blockSize_*(numBlocks_+1);
00422   
00423   if (blockSize_==1) {
00424     cs.resize( newsd );
00425     sn.resize( newsd );
00426   }
00427   else {
00428     beta.resize( newsd );
00429   }
00430   
00431   // Initialize the state storage
00432         TEUCHOS_TEST_FOR_EXCEPTION(blockSize_*numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
00433                            "Belos::BlockGmresIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
00434 
00435   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00436   if (V_ == Teuchos::null) {
00437     // Get the multivector that is not null.
00438     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00439     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00440            "Belos::BlockGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00441     V_ = MVT::Clone( *tmp, newsd );
00442   }
00443   else {
00444     // Generate V_ by cloning itself ONLY if more space is needed.
00445     if (MVT::GetNumberVecs(*V_) < newsd) {
00446       Teuchos::RCP<const MV> tmp = V_;
00447       V_ = MVT::Clone( *tmp, newsd );
00448     }
00449   }
00450   
00451   // Generate R_ only if it doesn't exist, otherwise resize it.
00452   if (R_ == Teuchos::null) {
00453     R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00454         }
00455         if (initHessenberg_) {
00456           R_->shape( newsd, newsd-blockSize_ );
00457         }
00458         else {
00459           if (R_->numRows() < newsd || R_->numCols() < newsd-blockSize_) {
00460       R_->shapeUninitialized( newsd, newsd-blockSize_ );
00461     }
00462         }
00463   
00464   // Generate H_ only if it doesn't exist, and we are keeping the upper Hessenberg matrix.
00465   if (keepHessenberg_) {
00466     if (H_ == Teuchos::null) {
00467       H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00468           }
00469           if (initHessenberg_) {
00470             H_->shape( newsd, newsd-blockSize_ );
00471           }
00472           else {
00473             if (H_->numRows() < newsd || H_->numCols() < newsd-blockSize_) {
00474         H_->shapeUninitialized( newsd, newsd-blockSize_ );
00475       }
00476           }
00477   }
00478   else {
00479     // Point H_ and R_ at the same object.
00480     H_ = R_;
00481   }
00482   
00483   // Generate z_ only if it doesn't exist, otherwise resize it.
00484   if (z_ == Teuchos::null) {
00485     z_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00486         }
00487         if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
00488     z_->shapeUninitialized( newsd, blockSize_ );
00489   }
00490   
00491   // State storage has now been initialized.
00492   stateStorageInitialized_ = true;
00493       }
00494     }
00495   }
00496 
00498   // Get the current update from this subspace.
00499   template <class ScalarType, class MV, class OP>
00500   Teuchos::RCP<MV> BlockGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00501   {
00502     //
00503     // If this is the first iteration of the Arnoldi factorization, 
00504     // there is no update, so return Teuchos::null. 
00505     //
00506     Teuchos::RCP<MV> currentUpdate = Teuchos::null;
00507     if (curDim_==0) { 
00508       return currentUpdate; 
00509     } else {
00510       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00511       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00512       Teuchos::BLAS<int,ScalarType> blas;
00513       currentUpdate = MVT::Clone( *V_, blockSize_ );
00514       //
00515       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00516       //
00517       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
00518       //
00519       //  Solve the least squares problem.
00520       //
00521       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00522      Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,  
00523      R_->values(), R_->stride(), y.values(), y.stride() );
00524       //
00525       //  Compute the current update.
00526       //
00527       std::vector<int> index(curDim_);
00528       for ( int i=0; i<curDim_; i++ ) {
00529         index[i] = i;
00530       }
00531       Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
00532       MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
00533     }
00534     return currentUpdate;
00535   }
00536 
00537 
00539   // Get the native residuals stored in this iteration.  
00540   // Note:  No residual std::vector will be returned by Gmres.
00541   template <class ScalarType, class MV, class OP>
00542   Teuchos::RCP<const MV> BlockGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00543   {
00544     //
00545     // NOTE: Make sure the incoming std::vector is the correct size!
00546     //
00547     if ( norms && (int)norms->size() < blockSize_ )                         
00548       norms->resize( blockSize_ );                                          
00549     
00550     if (norms) {
00551       Teuchos::BLAS<int,ScalarType> blas;
00552       for (int j=0; j<blockSize_; j++) {
00553         (*norms)[j] = blas.NRM2( blockSize_, &(*z_)(curDim_, j), 1);
00554       }
00555     }
00556     return Teuchos::null;
00557   }
00558 
00559   
00560 
00562   // Initialize this iteration object
00563   template <class ScalarType, class MV, class OP>
00564   void BlockGmresIter<ScalarType,MV,OP>::initializeGmres(GmresIterationState<ScalarType,MV> newstate)
00565   {
00566     // Initialize the state storage if it isn't already.
00567     if (!stateStorageInitialized_) 
00568       setStateSize();
00569 
00570     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00571            "Belos::BlockGmresIter::initialize(): Cannot initialize state storage!");
00572     
00573     // NOTE:  In BlockGmresIter, V and Z are required!!!  
00574     // inconsitent multivectors widths and lengths will not be tolerated, and
00575     // will be treated with exceptions.
00576     //
00577     std::string errstr("Belos::BlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
00578 
00579     if (newstate.V != Teuchos::null && newstate.z != Teuchos::null) {
00580 
00581       // initialize V_,z_, and curDim_
00582 
00583       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00584                           std::invalid_argument, errstr );
00585       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00586                           std::invalid_argument, errstr );
00587       TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*(numBlocks_+1),
00588                           std::invalid_argument, errstr );
00589 
00590       curDim_ = newstate.curDim;
00591       int lclDim = MVT::GetNumberVecs(*newstate.V);
00592 
00593       // check size of Z
00594       TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_, std::invalid_argument, errstr);
00595       
00596 
00597       // copy basis vectors from newstate into V
00598       if (newstate.V != V_) {
00599         // only copy over the first block and print a warning.
00600         if (curDim_ == 0 && lclDim > blockSize_) {
00601     om_->stream(Warnings) << "Belos::BlockGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00602         << "The block size however is only " << blockSize_ << std::endl
00603         << "The last " << lclDim - blockSize_ << " vectors will be discarded." << std::endl;
00604   }
00605         std::vector<int> nevind(curDim_+blockSize_);
00606         for (int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
00607   Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
00608   Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_, nevind );
00609         MVT::MvAddMv( 1.0, *newV, 0.0, *newV, *lclV );
00610 
00611         // done with local pointers
00612         lclV = Teuchos::null;
00613       }
00614 
00615       // put data into z_, make sure old information is not still hanging around.
00616       if (newstate.z != z_) {
00617         z_->putScalar();
00618         Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+blockSize_,blockSize_);
00619         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclZ;
00620         lclZ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
00621         lclZ->assign(newZ);
00622 
00623         // done with local pointers
00624         lclZ = Teuchos::null;
00625       }
00626 
00627     }
00628     else {
00629 
00630       TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
00631                          "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
00632 
00633       TEUCHOS_TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00634                          "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial norms z_0.");
00635     }
00636 
00637     // the solver is initialized
00638     initialized_ = true;
00639 
00640     /*
00641     if (om_->isVerbosity( Debug ) ) {
00642       // Check almost everything here
00643       CheckList chk;
00644       chk.checkV = true;
00645       chk.checkArn = true;
00646       chk.checkAux = true;
00647       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00648     }
00649     */
00650 
00651   }
00652 
00653 
00655   // Iterate until the status test informs us we should stop.
00656   template <class ScalarType, class MV, class OP>
00657   void BlockGmresIter<ScalarType,MV,OP>::iterate()
00658   {
00659     //
00660     // Allocate/initialize data structures
00661     //
00662     if (initialized_ == false) {
00663       initialize();
00664     }
00665     
00666     // Compute the current search dimension. 
00667     int searchDim = blockSize_*numBlocks_;
00668 
00670     // iterate until the status test tells us to stop.
00671     //
00672     // also break if our basis is full
00673     //
00674     while (stest_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
00675 
00676       iter_++;
00677 
00678       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
00679       int lclDim = curDim_ + blockSize_; 
00680 
00681       // Get the current part of the basis.
00682       std::vector<int> curind(blockSize_);
00683       for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
00684       Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
00685 
00686       // Get a view of the previous vectors.
00687       // This is used for orthogonalization and for computing V^H K H
00688       for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
00689       Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
00690 
00691       // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
00692       lp_->apply(*Vprev,*Vnext);
00693       Vprev = Teuchos::null;
00694       
00695       // Remove all previous Krylov basis vectors from Vnext      
00696       // Get a view of all the previous vectors
00697       std::vector<int> prevind(lclDim);
00698       for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00699       Vprev = MVT::CloneView(*V_,prevind);
00700       Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00701       
00702       // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00703       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00704   subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00705            ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
00706       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00707       AsubH.append( subH );
00708       
00709       // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00710       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00711   subH2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00712             ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
00713       subH2->putScalar();  // Initialize subdiagonal to zero
00714       int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
00715 
00716       // Copy over the coefficients if we are saving the upper Hessenberg matrix,
00717       // just in case we run into an error.
00718       if (keepHessenberg_) {
00719   // Copy over the orthogonalization coefficients.
00720   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00721     subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00722              ( Teuchos::View,*R_,lclDim,blockSize_,0,curDim_ ) );
00723   subR->assign(*subH);
00724   
00725   // Copy over the lower diagonal block of the Hessenberg matrix.
00726   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00727     subR2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00728         ( Teuchos::View,*R_,blockSize_,blockSize_,lclDim,curDim_ ) );
00729   subR2->assign(*subH2);
00730       }
00731 
00732       TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,GmresIterationOrthoFailure,
00733        "Belos::BlockGmresIter::iterate(): couldn't generate basis of full rank.");
00734       //
00735       // V has been extended, and H has been extended. 
00736       //
00737       // Update the QR factorization of the upper Hessenberg matrix
00738       //
00739       updateLSQR();
00740       //
00741       // Update basis dim and release all pointers.
00742       //
00743       Vnext = Teuchos::null;
00744       curDim_ += blockSize_;
00745       //        
00746       /*      
00747       // When required, monitor some orthogonalities
00748       if (om_->isVerbosity( Debug ) ) {
00749       // Check almost everything here
00750       CheckList chk;
00751       chk.checkV = true;
00752       chk.checkArn = true;
00753       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00754       }
00755       else if (om_->isVerbosity( OrthoDetails ) ) {
00756         CheckList chk;
00757         chk.checkV = true;
00758         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00759       }
00760       */ 
00761       
00762     } // end while (statusTest == false)
00763    
00764   }
00765 
00766   
00767   template<class ScalarType, class MV, class OP>
00768   void BlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00769   {
00770     int i, j, maxidx;
00771     ScalarType sigma, mu, vscale, maxelem;
00772     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00773     
00774     // Get correct dimension based on input "dim"
00775     // Remember that ortho failures result in an exit before updateLSQR() is called.
00776     // Therefore, it is possible that dim == curDim_.
00777     int curDim = curDim_;
00778     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00779       curDim = dim;
00780     }
00781     
00782     Teuchos::BLAS<int, ScalarType> blas;
00783     //
00784     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00785     // system to upper-triangular form.
00786     //
00787     if (blockSize_ == 1) {
00788       //
00789       // QR factorization of Least-Squares system with Givens rotations
00790       //
00791       for (i=0; i<curDim; i++) {
00792   //
00793   // Apply previous Givens rotations to new column of Hessenberg matrix
00794   //
00795   blas.ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
00796       }
00797       //
00798       // Calculate new Givens rotation
00799       //
00800       blas.ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
00801       (*R_)(curDim+1,curDim) = zero;
00802       //
00803       // Update RHS w/ new transformation
00804       //
00805       blas.ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
00806     }
00807     else {
00808       //
00809       // QR factorization of Least-Squares system with Householder reflectors
00810       //
00811       for (j=0; j<blockSize_; j++) {
00812   //
00813   // Apply previous Householder reflectors to new block of Hessenberg matrix
00814   //
00815   for (i=0; i<curDim+j; i++) {
00816     sigma = blas.DOT( blockSize_, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
00817     sigma += (*R_)(i,curDim+j);
00818     sigma *= beta[i];
00819     blas.AXPY(blockSize_, -sigma, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
00820     (*R_)(i,curDim+j) -= sigma;
00821   }
00822   //
00823   // Compute new Householder reflector
00824   //
00825   maxidx = blas.IAMAX( blockSize_+1, &(*R_)(curDim+j,curDim+j), 1 );
00826   maxelem = (*R_)(curDim+j+maxidx-1,curDim+j);
00827   for (i=0; i<blockSize_+1; i++)
00828     (*R_)(curDim+j+i,curDim+j) /= maxelem;
00829   sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j), 1,
00830         &(*R_)(curDim+j+1,curDim+j), 1 );
00831   if (sigma == zero) {
00832     beta[curDim + j] = zero;
00833   } else {
00834     mu = Teuchos::ScalarTraits<ScalarType>::squareroot((*R_)(curDim+j,curDim+j)*(*R_)(curDim+j,curDim+j)+sigma);
00835     if ( Teuchos::ScalarTraits<ScalarType>::real((*R_)(curDim+j,curDim+j)) 
00836          < Teuchos::ScalarTraits<MagnitudeType>::zero() ) {
00837       vscale = (*R_)(curDim+j,curDim+j) - mu;
00838     } else {
00839       vscale = -sigma / ((*R_)(curDim+j,curDim+j) + mu);
00840     }
00841     beta[curDim+j] = 2.0*vscale*vscale/(sigma + vscale*vscale);
00842     (*R_)(curDim+j,curDim+j) = maxelem*mu;
00843     for (i=0; i<blockSize_; i++)
00844       (*R_)(curDim+j+1+i,curDim+j) /= vscale;
00845   }
00846   //
00847   // Apply new Householder reflector to rhs
00848   //
00849   for (i=0; i<blockSize_; i++) {
00850     sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j),
00851           1, &(*z_)(curDim+j+1,i), 1);
00852     sigma += (*z_)(curDim+j,i);
00853     sigma *= beta[curDim+j];
00854     blas.AXPY(blockSize_, -sigma, &(*R_)(curDim+j+1,curDim+j),
00855         1, &(*z_)(curDim+j+1,i), 1);
00856     (*z_)(curDim+j,i) -= sigma;
00857   }
00858       }
00859     } // end if (blockSize_ == 1)
00860 
00861     // If the least-squares problem is updated wrt "dim" then update the curDim_.
00862     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00863       curDim_ = dim + blockSize_;
00864     }
00865 
00866   } // end updateLSQR()
00867 
00868 } // end Belos namespace
00869 
00870 #endif /* BELOS_BLOCK_GMRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines