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>::
00335   BlockFGmresIter (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00336                    const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00337                    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00338                    const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00339                    Teuchos::ParameterList &params ):
00340     lp_(problem),
00341     om_(printer),
00342     stest_(tester),
00343     ortho_(ortho),
00344     blockSize_(0),
00345     numBlocks_(0),
00346     initialized_(false),
00347     stateStorageInitialized_(false),
00348     curDim_(0),
00349     iter_(0)
00350   {
00351     // Get the maximum number of blocks allowed for this Krylov subspace
00352     TEUCHOS_TEST_FOR_EXCEPTION(
00353       ! params.isParameter ("Num Blocks"), std::invalid_argument,
00354       "Belos::BlockFGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00355     const int nb = params.get<int> ("Num Blocks");
00356 
00357     // Set the block size and allocate data.
00358     const int bs = params.get ("Block Size", 1);
00359     setSize (bs, nb);
00360   }
00361 
00363   // Set the block size and make necessary adjustments.
00364   template <class ScalarType, class MV, class OP>
00365   void BlockFGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
00366   {
00367     // This routine only allocates space; it doesn't not perform any computation
00368     // any change in size will invalidate the state of the solver.
00369 
00370     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockFGmresIter::setSize was passed a non-positive argument.");
00371     if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00372       // do nothing
00373       return;
00374     }
00375 
00376     if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
00377       stateStorageInitialized_ = false;
00378 
00379     blockSize_ = blockSize;
00380     numBlocks_ = numBlocks;
00381 
00382     initialized_ = false;
00383     curDim_ = 0;
00384 
00385     // Use the current blockSize_ and numBlocks_ to initialize the state storage.
00386     setStateSize();
00387 
00388   }
00389 
00391   // Setup the state storage.
00392   template <class ScalarType, class MV, class OP>
00393   void BlockFGmresIter<ScalarType,MV,OP>::setStateSize ()
00394   {
00395     using Teuchos::RCP;
00396     using Teuchos::rcp;
00397     typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
00398 
00399     if (! stateStorageInitialized_) {
00400       // Check if there is any multivector to clone from.
00401       RCP<const MV> lhsMV = lp_->getLHS();
00402       RCP<const MV> rhsMV = lp_->getRHS();
00403       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00404         stateStorageInitialized_ = false;
00405         return;
00406       }
00407       else {
00409         // blockSize*numBlocks dependent
00410         //
00411         int newsd = blockSize_*(numBlocks_+1);
00412 
00413         if (blockSize_==1) {
00414           cs.resize (newsd);
00415           sn.resize (newsd);
00416         }
00417         else {
00418           beta.resize (newsd);
00419         }
00420 
00421         // Initialize the state storage
00422         TEUCHOS_TEST_FOR_EXCEPTION(
00423           blockSize_ * static_cast<ptrdiff_t> (numBlocks_) > MVText::GetGlobalLength (*rhsMV),
00424           std::invalid_argument, "Belos::BlockFGmresIter::setStateSize(): "
00425           "Cannot generate a Krylov basis with dimension larger the operator!");
00426 
00427         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00428         if (V_ == Teuchos::null) {
00429           // Get the multivector that is not null.
00430           RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
00431           TEUCHOS_TEST_FOR_EXCEPTION(
00432             tmp == Teuchos::null, std::invalid_argument,
00433             "Belos::BlockFGmresIter::setStateSize(): "
00434             "linear problem does not specify multivectors to clone from.");
00435           V_ = MVT::Clone (*tmp, newsd);
00436         }
00437         else {
00438           // Generate V_ by cloning itself ONLY if more space is needed.
00439           if (MVT::GetNumberVecs (*V_) < newsd) {
00440             RCP<const MV> tmp = V_;
00441             V_ = MVT::Clone (*tmp, newsd);
00442           }
00443         }
00444 
00445         if (Z_ == Teuchos::null) {
00446           // Get the multivector that is not null.
00447           RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
00448           TEUCHOS_TEST_FOR_EXCEPTION(
00449             tmp == Teuchos::null, std::invalid_argument,
00450             "Belos::BlockFGmresIter::setStateSize(): "
00451             "linear problem does not specify multivectors to clone from.");
00452           Z_ = MVT::Clone (*tmp, newsd);
00453         }
00454         else {
00455           // Generate Z_ by cloning itself ONLY if more space is needed.
00456           if (MVT::GetNumberVecs (*Z_) < newsd) {
00457             RCP<const MV> tmp = Z_;
00458             Z_ = MVT::Clone (*tmp, newsd);
00459           }
00460         }
00461 
00462         // Generate H_ only if it doesn't exist, otherwise resize it.
00463         if (H_ == Teuchos::null) {
00464           H_ = rcp (new SDM (newsd, newsd-blockSize_));
00465         }
00466         else {
00467           H_->shapeUninitialized (newsd, newsd - blockSize_);
00468         }
00469 
00470         // TODO:  Insert logic so that Hessenberg matrix can be saved and reduced matrix is stored in R_
00471         //R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-blockSize_ ) );
00472         // Generate z_ only if it doesn't exist, otherwise resize it.
00473         if (z_ == Teuchos::null) {
00474           z_ = rcp (new SDM (newsd, blockSize_));
00475         }
00476         else {
00477           z_->shapeUninitialized (newsd, blockSize_);
00478         }
00479 
00480         // State storage has now been initialized.
00481         stateStorageInitialized_ = true;
00482       }
00483     }
00484   }
00485 
00486 
00487   template <class ScalarType, class MV, class OP>
00488   Teuchos::RCP<MV>
00489   BlockFGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00490   {
00491     typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
00492 
00493     Teuchos::RCP<MV> currentUpdate = Teuchos::null;
00494     if (curDim_ == 0) {
00495       // If this is the first iteration of the Arnoldi factorization,
00496       // then there is no update, so return Teuchos::null.
00497       return currentUpdate;
00498     }
00499     else {
00500       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero ();
00501       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one ();
00502       Teuchos::BLAS<int,ScalarType> blas;
00503 
00504       currentUpdate = MVT::Clone (*Z_, blockSize_);
00505 
00506       // Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00507       SDM y (Teuchos::Copy, *z_, curDim_, blockSize_);
00508 
00509       // Solve the least squares problem.
00510       blas.TRSM (Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00511                  Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
00512                  H_->values (), H_->stride (), y.values (), y.stride ());
00513 
00514       // Compute the current update.
00515       std::vector<int> index (curDim_);
00516       for (int i = 0; i < curDim_; ++i) {
00517         index[i] = i;
00518       }
00519       Teuchos::RCP<const MV> Zjp1 = MVT::CloneView (*Z_, index);
00520       MVT::MvTimesMatAddMv (one, *Zjp1, y, zero, *currentUpdate);
00521     }
00522     return currentUpdate;
00523   }
00524 
00525 
00526   template <class ScalarType, class MV, class OP>
00527   Teuchos::RCP<const MV>
00528   BlockFGmresIter<ScalarType,MV,OP>::
00529   getNativeResiduals (std::vector<MagnitudeType> *norms) const
00530   {
00531     // NOTE: Make sure the incoming std::vector is the correct size!
00532     if (norms != NULL && (int)norms->size() < blockSize_) {
00533       norms->resize (blockSize_);
00534     }
00535 
00536     if (norms != NULL) {
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 
00543     // FGmres does not return a residual (multi)vector.
00544     return Teuchos::null;
00545   }
00546 
00547 
00548   template <class ScalarType, class MV, class OP>
00549   void BlockFGmresIter<ScalarType,MV,OP>::
00550   initializeGmres (GmresIterationState<ScalarType,MV> newstate)
00551   {
00552     using Teuchos::RCP;
00553     using Teuchos::rcp;
00554     using std::endl;
00555     typedef Teuchos::ScalarTraits<ScalarType> STS;
00556     typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
00557     const ScalarType ZERO = STS::zero ();
00558     const ScalarType ONE = STS::one ();
00559 
00560     // Initialize the state storage if it isn't already.
00561     if (! stateStorageInitialized_) {
00562       setStateSize ();
00563     }
00564 
00565     TEUCHOS_TEST_FOR_EXCEPTION(
00566       ! stateStorageInitialized_, std::invalid_argument,
00567       "Belos::BlockFGmresIter::initialize(): Cannot initialize state storage!");
00568 
00569     // NOTE: In BlockFGmresIter, V and Z are required!!!  Inconsistent
00570     // multivectors widths and lengths will not be tolerated, and will
00571     // be treated with exceptions.
00572     const char errstr[] = "Belos::BlockFGmresIter::initialize(): The given "
00573       "multivectors must have a consistent length and width.";
00574 
00575     if (! newstate.V.is_null () && ! newstate.z.is_null ()) {
00576 
00577       // initialize V_,z_, and curDim_
00578 
00579       TEUCHOS_TEST_FOR_EXCEPTION(
00580         MVText::GetGlobalLength(*newstate.V) != MVText::GetGlobalLength(*V_),
00581         std::invalid_argument, errstr );
00582       TEUCHOS_TEST_FOR_EXCEPTION(
00583         MVT::GetNumberVecs(*newstate.V) < blockSize_,
00584         std::invalid_argument, errstr );
00585       TEUCHOS_TEST_FOR_EXCEPTION(
00586         newstate.curDim > blockSize_*(numBlocks_+1),
00587         std::invalid_argument, errstr );
00588 
00589       curDim_ = newstate.curDim;
00590       const int lclDim = MVT::GetNumberVecs(*newstate.V);
00591 
00592       // check size of Z
00593       TEUCHOS_TEST_FOR_EXCEPTION(
00594         newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_,
00595         std::invalid_argument, errstr);
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           std::ostream& warn = om_->stream (Warnings);
00602           warn << "Belos::BlockFGmresIter::initialize(): the solver was "
00603                << "initialized with a kernel of " << lclDim << endl
00604                << "The block size however is only " << blockSize_ << endl
00605                << "The last " << lclDim - blockSize_
00606                << " vectors will be discarded." << endl;
00607         }
00608         std::vector<int> nevind (curDim_ + blockSize_);
00609         for (int i = 0; i < curDim_ + blockSize_; ++i) {
00610           nevind[i] = i;
00611         }
00612         RCP<const MV> newV = MVT::CloneView (*newstate.V, nevind);
00613         RCP<MV> lclV = MVT::CloneViewNonConst (*V_, nevind);
00614         MVT::MvAddMv (ONE, *newV, ZERO, *newV, *lclV);
00615 
00616         // done with local pointers
00617         lclV = Teuchos::null;
00618       }
00619 
00620       // put data into z_, make sure old information is not still hanging around.
00621       if (newstate.z != z_) {
00622         z_->putScalar();
00623         SDM newZ (Teuchos::View, *newstate.z, curDim_ + blockSize_, blockSize_);
00624         RCP<SDM> lclz;
00625         lclz = rcp (new SDM (Teuchos::View, *z_, curDim_ + blockSize_, blockSize_));
00626         lclz->assign (newZ);
00627         lclz = Teuchos::null; // done with local pointers
00628       }
00629     }
00630     else {
00631       TEUCHOS_TEST_FOR_EXCEPTION(
00632         newstate.V == Teuchos::null,std::invalid_argument,
00633         "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial kernel V_0.");
00634 
00635       TEUCHOS_TEST_FOR_EXCEPTION(
00636         newstate.z == Teuchos::null,std::invalid_argument,
00637         "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial norms z_0.");
00638     }
00639 
00640     // the solver is initialized
00641     initialized_ = true;
00642   }
00643 
00644 
00645   template <class ScalarType, class MV, class OP>
00646   void BlockFGmresIter<ScalarType,MV,OP>::iterate()
00647   {
00648     using Teuchos::Array;
00649     using Teuchos::null;
00650     using Teuchos::RCP;
00651     using Teuchos::rcp;
00652     using Teuchos::View;
00653     typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
00654 
00655     // Allocate/initialize data structures
00656     if (initialized_ == false) {
00657       initialize();
00658     }
00659 
00660     // Compute the current search dimension.
00661     const int searchDim = blockSize_ * numBlocks_;
00662 
00663     // Iterate until the status test tells us to stop.
00664     // Raise an exception if a computed block is not full rank.
00665     while (stest_->checkStatus (this) != Passed && curDim_+blockSize_ <= searchDim) {
00666       ++iter_;
00667 
00668       // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
00669       const int lclDim = curDim_ + blockSize_;
00670 
00671       // Get the current part of the basis.
00672       std::vector<int> curind (blockSize_);
00673       for (int i = 0; i < blockSize_; ++i) {
00674         curind[i] = lclDim + i;
00675       }
00676       RCP<MV> Vnext = MVT::CloneViewNonConst (*V_, curind);
00677 
00678       // Get a view of the previous vectors.
00679       // This is used for orthogonalization and for computing V^H K H.
00680       for (int i = 0; i < blockSize_; ++i) {
00681         curind[i] = curDim_ + i;
00682       }
00683       RCP<const MV> Vprev = MVT::CloneView (*V_, curind);
00684       RCP<MV> Znext = MVT::CloneViewNonConst (*Z_, curind);
00685 
00686       // Compute the next (multi)vector in the Krylov basis:  Znext = M*Vprev
00687       lp_->applyRightPrec (*Vprev, *Znext);
00688       Vprev = null;
00689 
00690       // Compute the next (multi)vector in the Krylov basis:  Vnext = A*Znext
00691       lp_->applyOp (*Znext, *Vnext);
00692       Znext = null;
00693 
00694       // Remove all previous Krylov basis vectors from Vnext
00695       // Get a view of all the previous vectors
00696       std::vector<int> prevind (lclDim);
00697       for (int i = 0; i < lclDim; ++i) {
00698         prevind[i] = i;
00699       }
00700       Vprev = MVT::CloneView (*V_, prevind);
00701       Array<RCP<const MV> > AVprev (1, Vprev);
00702 
00703       // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00704       RCP<SDM> subH = rcp (new SDM (View, *H_, lclDim, blockSize_, 0, curDim_));
00705       Array<RCP<SDM> > AsubH;
00706       AsubH.append (subH);
00707 
00708       // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00709       RCP<SDM> subR = rcp (new SDM (View, *H_, blockSize_, blockSize_, lclDim, curDim_));
00710       const int rank = ortho_->projectAndNormalize (*Vnext, AsubH, subR, AVprev);
00711       TEUCHOS_TEST_FOR_EXCEPTION(
00712         rank != blockSize_, GmresIterationOrthoFailure,
00713         "Belos::BlockFGmresIter::iterate(): After orthogonalization, the new "
00714         "basis block does not have full rank.  It contains " << blockSize_
00715         << " vector" << (blockSize_ != 1 ? "s" : "")
00716         << ", but its rank is " << rank << ".");
00717 
00718       //
00719       // V has been extended, and H has been extended.
00720       //
00721       // Update the QR factorization of the upper Hessenberg matrix
00722       //
00723       updateLSQR ();
00724       //
00725       // Update basis dim and release all pointers.
00726       //
00727       Vnext = null;
00728       curDim_ += blockSize_;
00729     } // end while (statusTest == false)
00730   }
00731 
00732 
00733   template<class ScalarType, class MV, class OP>
00734   void BlockFGmresIter<ScalarType,MV,OP>::updateLSQR (int dim)
00735   {
00736     typedef Teuchos::ScalarTraits<ScalarType> STS;
00737     typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00738 
00739     const ScalarType zero = STS::zero ();
00740     const ScalarType two = (STS::one () + STS::one());
00741     ScalarType sigma, mu, vscale, maxelem;
00742     Teuchos::BLAS<int, ScalarType> blas;
00743 
00744     // Get correct dimension based on input 'dim'.  Remember that
00745     // orthogonalization failures result in an exit before
00746     // updateLSQR() is called.  Therefore, it is possible that dim ==
00747     // curDim_.
00748     int curDim = curDim_;
00749     if (dim >= curDim_ && dim < getMaxSubspaceDim ()) {
00750       curDim = dim;
00751     }
00752 
00753     // Apply previous transformations, and compute new transformation
00754     // to reduce upper Hessenberg system to upper triangular form.
00755     // The type of transformation we use depends the block size.  We
00756     // use Givens rotations for a block size of 1, and Householder
00757     // reflectors otherwise.
00758     if (blockSize_ == 1) {
00759       // QR factorization of upper Hessenberg matrix using Givens rotations
00760       for (int i = 0; i < curDim; ++i) {
00761         // Apply previous Givens rotations to new column of Hessenberg matrix
00762         blas.ROT (1, &(*H_)(i, curDim), 1, &(*H_)(i+1, curDim), 1, &cs[i], &sn[i]);
00763       }
00764       // Calculate new Givens rotation
00765       blas.ROTG (&(*H_)(curDim, curDim), &(*H_)(curDim+1, curDim), &cs[curDim], &sn[curDim]);
00766       (*H_)(curDim+1, curDim) = zero;
00767 
00768       // Update RHS w/ new transformation
00769       blas.ROT (1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim]);
00770     }
00771     else {
00772       // QR factorization of least-squares system using Householder reflectors.
00773       for (int j = 0; j < blockSize_; ++j) {
00774         // Apply previous Householder reflectors to new block of Hessenberg matrix
00775         for (int i = 0; i < curDim + j; ++i) {
00776           sigma = blas.DOT (blockSize_, &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
00777           sigma += (*H_)(i,curDim+j);
00778           sigma *= beta[i];
00779           blas.AXPY (blockSize_, -sigma, &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
00780           (*H_)(i,curDim+j) -= sigma;
00781         }
00782 
00783         // Compute new Householder reflector
00784         const int maxidx = blas.IAMAX (blockSize_+1, &(*H_)(curDim+j,curDim+j), 1);
00785         maxelem = (*H_)(curDim + j + maxidx - 1, curDim + j);
00786         for (int i = 0; i < blockSize_ + 1; ++i) {
00787           (*H_)(curDim+j+i,curDim+j) /= maxelem;
00788         }
00789         sigma = blas.DOT (blockSize_, &(*H_)(curDim + j + 1, curDim + j), 1,
00790                           &(*H_)(curDim + j + 1, curDim + j), 1);
00791         if (sigma == zero) {
00792           beta[curDim + j] = zero;
00793         } else {
00794           mu = STS::squareroot ((*H_)(curDim+j,curDim+j)*(*H_)(curDim+j,curDim+j)+sigma);
00795           if (STS::real ((*H_)(curDim + j, curDim + j)) < STM::zero ()) {
00796             vscale = (*H_)(curDim+j,curDim+j) - mu;
00797           } else {
00798             vscale = -sigma / ((*H_)(curDim+j, curDim+j) + mu);
00799           }
00800           beta[curDim+j] = two * vscale * vscale / (sigma + vscale*vscale);
00801           (*H_)(curDim+j, curDim+j) = maxelem*mu;
00802           for (int i = 0; i < blockSize_; ++i) {
00803             (*H_)(curDim+j+1+i,curDim+j) /= vscale;
00804           }
00805         }
00806 
00807         // Apply new Householder reflector to the right-hand side.
00808         for (int i = 0; i < blockSize_; ++i) {
00809           sigma = blas.DOT (blockSize_, &(*H_)(curDim+j+1,curDim+j),
00810                             1, &(*z_)(curDim+j+1,i), 1);
00811           sigma += (*z_)(curDim+j,i);
00812           sigma *= beta[curDim+j];
00813           blas.AXPY (blockSize_, -sigma, &(*H_)(curDim+j+1,curDim+j),
00814                      1, &(*z_)(curDim+j+1,i), 1);
00815           (*z_)(curDim+j,i) -= sigma;
00816         }
00817       }
00818     } // end if (blockSize_ == 1)
00819 
00820     // If the least-squares problem is updated wrt "dim" then update curDim_.
00821     if (dim >= curDim_ && dim < getMaxSubspaceDim ()) {
00822       curDim_ = dim + blockSize_;
00823     }
00824   } // end updateLSQR()
00825 
00826 } // namespace Belos
00827 
00828 #endif /* BELOS_BLOCK_FGMRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines