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