Belos Version of the Day
BelosBlockFGmresIter.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_FGMRES_ITER_HPP
00043 #define BELOS_BLOCK_FGMRES_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 BlockFGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
00084 
00085   public:
00086     
00087   //
00088   // Convenience typedefs
00089   //
00090   typedef MultiVecTraits<ScalarType,MV> MVT;
00091   typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00092   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00093   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00094   typedef typename SCT::magnitudeType MagnitudeType;
00095 
00097 
00098 
00108   BlockFGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00109        const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00110        const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00111        const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00112        Teuchos::ParameterList &params );
00113   
00115   virtual ~BlockFGmresIter() {};
00117   
00118   
00120 
00121   
00143   void iterate();
00144 
00166   void initializeGmres(GmresIterationState<ScalarType,MV> newstate);
00167 
00171   void initialize()
00172   {
00173     GmresIterationState<ScalarType,MV> empty;
00174     initializeGmres(empty);
00175   }
00176   
00184   GmresIterationState<ScalarType,MV> getState() const {
00185     GmresIterationState<ScalarType,MV> state;
00186     state.curDim = curDim_;
00187     state.V = V_;
00188     state.Z = Z_;
00189     state.H = H_;
00190     state.R = R_;
00191     state.z = z_;
00192     return state;
00193   }
00194 
00196 
00197   
00199 
00200 
00202   int getNumIters() const { return iter_; }
00203   
00205   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00206 
00209   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00210 
00212 
00217   Teuchos::RCP<MV> getCurrentUpdate() const;
00218 
00220 
00223   void updateLSQR( int dim = -1 );
00224 
00226   int getCurSubspaceDim() const { 
00227     if (!initialized_) return 0;
00228     return curDim_;
00229   };
00230 
00232   int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }
00233 
00235 
00236   
00238 
00239 
00241   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00242 
00244   int getBlockSize() const { return blockSize_; }
00245   
00247   void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }
00248   
00250   int getNumBlocks() const { return numBlocks_; }
00251   
00253   void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }
00254   
00261   void setSize(int blockSize, int numBlocks);
00262 
00264   bool isInitialized() { return initialized_; }
00265 
00267 
00268   private:
00269 
00270   //
00271   // Internal methods
00272   //
00274   void setStateSize();
00275   
00276   //
00277   // Classes inputed through constructor that define the linear problem to be solved.
00278   //
00279   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00280   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00281   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00282   const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00283 
00284   //
00285   // Algorithmic parameters
00286   //  
00287   // blockSize_ is the solver block size.
00288   // It controls the number of vectors added to the basis on each iteration.
00289   int blockSize_;
00290   // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00291   int numBlocks_; 
00292   
00293   // Storage for QR factorization of the least squares system.
00294   Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
00295   Teuchos::SerialDenseVector<int,MagnitudeType> cs;
00296   
00297   // 
00298   // Current solver state
00299   //
00300   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00301   // is capable of running; _initialize is controlled  by the initialize() member method
00302   // For the implications of the state of initialized_, please see documentation for initialize()
00303   bool initialized_;
00304 
00305   // stateStorageInitialized_ specified that the state storage has be initialized to the current
00306   // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00307   // generated without the right-hand side or solution vectors.
00308   bool stateStorageInitialized_;
00309 
00310   // Current subspace dimension, and number of iterations performed.
00311   int curDim_, iter_;
00312   
00313   // 
00314   // State Storage
00315   //
00316   Teuchos::RCP<MV> V_;
00317   Teuchos::RCP<MV> Z_;
00318   //
00319   // Projected matrices
00320   // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00321   //
00322   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00323   // 
00324   // QR decomposition of Projected matrices for solving the least squares system HY = B.
00325   // R_: Upper triangular reduction of H
00326   // z_: Q applied to right-hand side of the least squares system
00327   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00328   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;  
00329 };
00330 
00332   // Constructor.
00333   template<class ScalarType, class MV, class OP>
00334   BlockFGmresIter<ScalarType,MV,OP>::BlockFGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00335                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00336                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00337                const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00338                Teuchos::ParameterList &params ):
00339     lp_(problem),
00340     om_(printer),
00341     stest_(tester),
00342     ortho_(ortho),
00343     blockSize_(0),
00344     numBlocks_(0),
00345     initialized_(false),
00346     stateStorageInitialized_(false),
00347     curDim_(0),
00348     iter_(0)
00349   {
00350     // Get the maximum number of blocks allowed for this Krylov subspace
00351     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00352                        "Belos::BlockFGmresIter::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 BlockFGmresIter<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     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockFGmresIter::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 BlockFGmresIter<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         TEUCHOS_TEST_FOR_EXCEPTION( blockSize_*static_cast<ptrdiff_t>(numBlocks_) > MVText::GetGlobalLength(*rhsMV),std::invalid_argument,
00419                            "Belos::BlockFGmresIter::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     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00426            "Belos::BlockFGmresIter::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   if (Z_ == Teuchos::null) {
00438     // Get the multivector that is not null.
00439     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00440     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00441            "Belos::BlockFGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00442     Z_ = MVT::Clone( *tmp, newsd );
00443   }
00444   else {
00445     // Generate Z_ by cloning itself ONLY if more space is needed.
00446     if (MVT::GetNumberVecs(*Z_) < newsd) {
00447       Teuchos::RCP<const MV> tmp = Z_;
00448       Z_ = MVT::Clone( *tmp, newsd );
00449     }
00450   } 
00451 
00452   // Generate H_ only if it doesn't exist, otherwise resize it.
00453   if (H_ == Teuchos::null)
00454     H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-blockSize_ ) ); 
00455   else
00456     H_->shapeUninitialized( newsd, newsd-blockSize_ );
00457   
00458   // TODO:  Insert logic so that Hessenberg matrix can be saved and reduced matrix is stored in R_
00459   //R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-blockSize_ ) );
00460   // Generate z_ only if it doesn't exist, otherwise resize it.
00461   if (z_ == Teuchos::null)
00462     z_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,blockSize_) );
00463   else
00464     z_->shapeUninitialized( newsd, blockSize_ );
00465   
00466   // State storage has now been initialized.
00467   stateStorageInitialized_ = true;
00468       }
00469     }
00470   }
00471 
00473   // Get the current update from this subspace.
00474   template <class ScalarType, class MV, class OP>
00475   Teuchos::RCP<MV> BlockFGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00476   {
00477     //
00478     // If this is the first iteration of the Arnoldi factorization, 
00479     // there is no update, so return Teuchos::null. 
00480     //
00481     Teuchos::RCP<MV> currentUpdate = Teuchos::null;
00482     if (curDim_==0) { 
00483       return currentUpdate; 
00484     } else {
00485       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00486       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00487       Teuchos::BLAS<int,ScalarType> blas;
00488       currentUpdate = MVT::Clone( *Z_, blockSize_ );
00489       //
00490       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00491       //
00492       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
00493       //
00494       //  Solve the least squares problem.
00495       //
00496       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00497      Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,  
00498      H_->values(), H_->stride(), y.values(), y.stride() );
00499       //
00500       //  Compute the current update.
00501       //
00502       std::vector<int> index(curDim_);
00503       for ( int i=0; i<curDim_; i++ ) {   
00504         index[i] = i;
00505       }
00506       Teuchos::RCP<const MV> Zjp1 = MVT::CloneView( *Z_, index );
00507       MVT::MvTimesMatAddMv( one, *Zjp1, y, zero, *currentUpdate );
00508     }
00509     return currentUpdate;
00510   }
00511 
00512 
00514   // Get the native residuals stored in this iteration.  
00515   // Note:  No residual std::vector will be returned by FGmres.
00516   template <class ScalarType, class MV, class OP>
00517   Teuchos::RCP<const MV> BlockFGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00518   {
00519     //
00520     // NOTE: Make sure the incoming std::vector is the correct size!
00521     //
00522     if ( norms && (int)norms->size() < blockSize_ )                         
00523       norms->resize( blockSize_ );                                          
00524     
00525     if (norms) {
00526       Teuchos::BLAS<int,ScalarType> blas;
00527       for (int j=0; j<blockSize_; j++) {
00528         (*norms)[j] = blas.NRM2( blockSize_, &(*z_)(curDim_, j), 1);
00529       }
00530     }
00531     return Teuchos::null;
00532   }
00533 
00534   
00535 
00537   // Initialize this iteration object
00538   template <class ScalarType, class MV, class OP>
00539   void BlockFGmresIter<ScalarType,MV,OP>::initializeGmres(GmresIterationState<ScalarType,MV> newstate)
00540   {
00541     // Initialize the state storage if it isn't already.
00542     if (!stateStorageInitialized_) 
00543       setStateSize();
00544 
00545     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00546            "Belos::BlockFGmresIter::initialize(): Cannot initialize state storage!");
00547     
00548     // NOTE:  In BlockFGmresIter, V and Z are required!!!  
00549     // inconsitent multivectors widths and lengths will not be tolerated, and
00550     // will be treated with exceptions.
00551     //
00552     std::string errstr("Belos::BlockFGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
00553 
00554     if (newstate.V != Teuchos::null && newstate.z != Teuchos::null) {
00555 
00556       // initialize V_,z_, and curDim_
00557 
00558       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.V) != MVText::GetGlobalLength(*V_),
00559                           std::invalid_argument, errstr );
00560       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00561                           std::invalid_argument, errstr );
00562       TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*(numBlocks_+1),
00563                           std::invalid_argument, errstr );
00564 
00565       curDim_ = newstate.curDim;
00566       int lclDim = MVT::GetNumberVecs(*newstate.V);
00567 
00568       // check size of Z
00569       TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_, std::invalid_argument, errstr);
00570       
00571 
00572       // copy basis vectors from newstate into V
00573       if (newstate.V != V_) {
00574         // only copy over the first block and print a warning.
00575         if (curDim_ == 0 && lclDim > blockSize_) {
00576     om_->stream(Warnings) << "Belos::BlockFGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00577         << "The block size however is only " << blockSize_ << std::endl
00578         << "The last " << lclDim - blockSize_ << " vectors will be discarded." << std::endl;
00579   }
00580         std::vector<int> nevind(curDim_+blockSize_);
00581         for (int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
00582   Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
00583   Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_, nevind );
00584         MVT::MvAddMv( 1.0, *newV, 0.0, *newV, *lclV );
00585 
00586         // done with local pointers
00587         lclV = Teuchos::null;
00588       }
00589 
00590       // put data into z_, make sure old information is not still hanging around.
00591       if (newstate.z != z_) {
00592         z_->putScalar();
00593         Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+blockSize_,blockSize_);
00594         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclz;
00595         lclz = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
00596         lclz->assign(newZ);
00597 
00598         // done with local pointers
00599         lclz = Teuchos::null;
00600       }
00601 
00602     }
00603     else {
00604 
00605       TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
00606                          "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial kernel V_0.");
00607 
00608       TEUCHOS_TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00609                          "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial norms z_0.");
00610     }
00611 
00612     // the solver is initialized
00613     initialized_ = true;
00614 
00615     /*
00616     if (om_->isVerbosity( Debug ) ) {
00617       // Check almost everything here
00618       CheckList chk;
00619       chk.checkV = true;
00620       chk.checkArn = true;
00621       chk.checkAux = true;
00622       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00623     }
00624     */
00625 
00626   }
00627 
00628 
00630   // Iterate until the status test informs us we should stop.
00631   template <class ScalarType, class MV, class OP>
00632   void BlockFGmresIter<ScalarType,MV,OP>::iterate()
00633   {
00634     //
00635     // Allocate/initialize data structures
00636     //
00637     if (initialized_ == false) {
00638       initialize();
00639     }
00640     
00641     // Compute the current search dimension. 
00642     int searchDim = blockSize_*numBlocks_;
00643 
00645     // iterate until the status test tells us to stop.
00646     //
00647     // also break if our basis is full
00648     //
00649     while (stest_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
00650 
00651       iter_++;
00652 
00653       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
00654       int lclDim = curDim_ + blockSize_; 
00655 
00656       // Get the current part of the basis.
00657       std::vector<int> curind(blockSize_);
00658       for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
00659       Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
00660 
00661       // Get a view of the previous vectors.
00662       // This is used for orthogonalization and for computing V^H K H
00663       for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
00664       Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
00665       Teuchos::RCP<MV> Znext = MVT::CloneViewNonConst(*Z_,curind);
00666 
00667       // Compute the next std::vector in the Krylov basis:  Znext = M*Vprev
00668       lp_->applyRightPrec(*Vprev,*Znext);
00669       Vprev = Teuchos::null;
00670 
00671       // Compute the next std::vector in the Krylov basis:  Vnext = A*Znext
00672       lp_->applyOp(*Znext,*Vnext);
00673       Znext = Teuchos::null;
00674       
00675       // Remove all previous Krylov basis vectors from Vnext      
00676       // Get a view of all the previous vectors
00677       std::vector<int> prevind(lclDim);
00678       for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00679       Vprev = MVT::CloneView(*V_,prevind);
00680       Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00681       
00682       // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00683       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00684   subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00685            ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
00686       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00687       AsubH.append( subH );
00688       
00689       // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00690       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00691   subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00692            ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
00693       int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00694       TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,GmresIterationOrthoFailure,
00695        "Belos::BlockFGmresIter::iterate(): couldn't generate basis of full rank.");
00696       //
00697       // V has been extended, and H has been extended. 
00698       //
00699       // Update the QR factorization of the upper Hessenberg matrix
00700       //
00701       updateLSQR();
00702       //
00703       // Update basis dim and release all pointers.
00704       //
00705       Vnext = Teuchos::null;
00706       curDim_ += blockSize_;
00707       //        
00708       /*      
00709       // When required, monitor some orthogonalities
00710       if (om_->isVerbosity( Debug ) ) {
00711       // Check almost everything here
00712       CheckList chk;
00713       chk.checkV = true;
00714       chk.checkArn = true;
00715       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00716       }
00717       else if (om_->isVerbosity( OrthoDetails ) ) {
00718         CheckList chk;
00719         chk.checkV = true;
00720         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00721       }
00722       */ 
00723       
00724     } // end while (statusTest == false)
00725    
00726   }
00727 
00728   
00729   template<class ScalarType, class MV, class OP>
00730   void BlockFGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00731   {
00732     int i, j, maxidx;
00733     ScalarType sigma, mu, vscale, maxelem;
00734     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00735     
00736     // Get correct dimension based on input "dim"
00737     // Remember that ortho failures result in an exit before updateLSQR() is called.
00738     // Therefore, it is possible that dim == curDim_.
00739     int curDim = curDim_;
00740     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00741       curDim = dim;
00742     }
00743     
00744     Teuchos::BLAS<int, ScalarType> blas;
00745     //
00746     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00747     // system to upper-triangular form.
00748     //
00749     if (blockSize_ == 1) {
00750       //
00751       // QR factorization of Least-Squares system with Givens rotations
00752       //
00753       for (i=0; i<curDim; i++) {
00754   //
00755   // Apply previous Givens rotations to new column of Hessenberg matrix
00756   //
00757   blas.ROT( 1, &(*H_)(i,curDim), 1, &(*H_)(i+1, curDim), 1, &cs[i], &sn[i] );
00758       }
00759       //
00760       // Calculate new Givens rotation
00761       //
00762       blas.ROTG( &(*H_)(curDim,curDim), &(*H_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
00763       (*H_)(curDim+1,curDim) = zero;
00764       //
00765       // Update RHS w/ new transformation
00766       //
00767       blas.ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
00768     }
00769     else {
00770       //
00771       // QR factorization of Least-Squares system with Householder reflectors
00772       //
00773       for (j=0; j<blockSize_; j++) {
00774   //
00775   // Apply previous Householder reflectors to new block of Hessenberg matrix
00776   //
00777   for (i=0; i<curDim+j; i++) {
00778     sigma = blas.DOT( blockSize_, &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
00779     sigma += (*H_)(i,curDim+j);
00780     sigma *= beta[i];
00781     blas.AXPY(blockSize_, -sigma, &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
00782     (*H_)(i,curDim+j) -= sigma;
00783   }
00784   //
00785   // Compute new Householder reflector
00786   //
00787   maxidx = blas.IAMAX( blockSize_+1, &(*H_)(curDim+j,curDim+j), 1 );
00788   maxelem = (*H_)(curDim+j+maxidx-1,curDim+j);
00789   for (i=0; i<blockSize_+1; i++)
00790     (*H_)(curDim+j+i,curDim+j) /= maxelem;
00791   sigma = blas.DOT( blockSize_, &(*H_)(curDim+j+1,curDim+j), 1,
00792         &(*H_)(curDim+j+1,curDim+j), 1 );
00793   if (sigma == zero) {
00794     beta[curDim + j] = zero;
00795   } else {
00796     mu = Teuchos::ScalarTraits<ScalarType>::squareroot((*H_)(curDim+j,curDim+j)*(*H_)(curDim+j,curDim+j)+sigma);
00797     if ( Teuchos::ScalarTraits<ScalarType>::real((*H_)(curDim+j,curDim+j)) 
00798          < Teuchos::ScalarTraits<MagnitudeType>::zero() ) {
00799       vscale = (*H_)(curDim+j,curDim+j) - mu;
00800     } else {
00801       vscale = -sigma / ((*H_)(curDim+j,curDim+j) + mu);
00802     }
00803     beta[curDim+j] = 2.0*vscale*vscale/(sigma + vscale*vscale);
00804     (*H_)(curDim+j,curDim+j) = maxelem*mu;
00805     for (i=0; i<blockSize_; i++)
00806       (*H_)(curDim+j+1+i,curDim+j) /= vscale;
00807   }
00808   //
00809   // Apply new Householder reflector to rhs
00810   //
00811   for (i=0; i<blockSize_; i++) {
00812     sigma = blas.DOT( blockSize_, &(*H_)(curDim+j+1,curDim+j),
00813           1, &(*z_)(curDim+j+1,i), 1);
00814     sigma += (*z_)(curDim+j,i);
00815     sigma *= beta[curDim+j];
00816     blas.AXPY(blockSize_, -sigma, &(*H_)(curDim+j+1,curDim+j),
00817         1, &(*z_)(curDim+j+1,i), 1);
00818     (*z_)(curDim+j,i) -= sigma;
00819   }
00820       }
00821     } // end if (blockSize_ == 1)
00822 
00823     // If the least-squares problem is updated wrt "dim" then update the curDim_.
00824     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00825       curDim_ = dim + blockSize_;
00826     }
00827 
00828   } // end updateLSQR()
00829 
00830 } // end Belos namespace
00831 
00832 #endif /* BELOS_BLOCK_FGMRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines