BelosGCRODRIter.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_GCRODR_ITER_HPP
00030 #define BELOS_GCRODR_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosMatOrthoManager.hpp"
00041 #include "BelosOutputManager.hpp"
00042 #include "BelosStatusTest.hpp"
00043 #include "BelosOperatorTraits.hpp"
00044 #include "BelosMultiVecTraits.hpp"
00045 
00046 #include "Teuchos_BLAS.hpp"
00047 #include "Teuchos_LAPACK.hpp"
00048 #include "Teuchos_SerialDenseMatrix.hpp"
00049 #include "Teuchos_SerialDenseVector.hpp"
00050 #include "Teuchos_ScalarTraits.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00065 namespace Belos {
00066   
00068 
00069   
00074   template <class ScalarType, class MV>
00075   struct GCRODRIterState {
00080     int curDim;
00081     
00083     Teuchos::RCP<const MV> V;
00084    
00086     Teuchos::RCP<const MV> U, C;
00087 
00093     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00094 
00097     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > B;
00098 
00100     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > R;
00101 
00103     Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > z;
00104 
00105     GCRODRIterState() : curDim(0), V(Teuchos::null), 
00106       U(Teuchos::null), C(Teuchos::null),
00107       H(Teuchos::null), R(Teuchos::null),
00108       z(Teuchos::null)
00109     {}
00110   };
00111   
00113   
00115 
00116   
00128   class GCRODRIterInitFailure : public BelosError {public:
00129     GCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00130     {}};
00131   
00138   class GCRODRIterOrthoFailure : public BelosError {public:
00139     GCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00140     {}};
00141   
00148   class GCRODRIterLAPACKFailure : public BelosError {public:
00149     GCRODRIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00150     {}};
00151   
00153   
00154   
00155   template<class ScalarType, class MV, class OP>
00156   class GCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
00157     
00158   public:
00159     
00160     //
00161     // Convenience typedefs
00162     //
00163     typedef MultiVecTraits<ScalarType,MV> MVT;
00164     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00165     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00166     typedef typename SCT::magnitudeType MagnitudeType;
00167     
00169 
00170     
00179     GCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00180     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00181     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00182     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00183     Teuchos::ParameterList &params );
00184     
00186     virtual ~GCRODRIter() {};
00188     
00189     
00191 
00192     
00214     void iterate();
00215     
00237     void initialize(GCRODRIterState<ScalarType,MV> newstate);
00238     
00242     void initialize()
00243     {
00244       GCRODRIterState<ScalarType,MV> empty;
00245       initialize(empty);
00246     }
00247     
00255     GCRODRIterState<ScalarType,MV> getState() const {
00256       GCRODRIterState<ScalarType,MV> state;
00257       state.curDim = curDim_;
00258       state.V = V_;
00259       state.U = U_;
00260       state.C = C_;
00261       state.H = H_;
00262       state.B = B_;
00263       state.R = R_;
00264       state.z = z_;
00265       return state;
00266     }
00267     
00269     
00270     
00272 
00273     
00275     int getNumIters() const { return iter_; }
00276     
00278     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00279     
00282     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00283     
00285 
00290     Teuchos::RCP<MV> getCurrentUpdate() const;
00291     
00293 
00296     void updateLSQR( int dim = -1 );
00297     
00299     int getCurSubspaceDim() const { 
00300       if (!initialized_) return 0;
00301       return curDim_;
00302     };
00303     
00305     int getMaxSubspaceDim() const { return numBlocks_-recycledBlocks_; }
00306     
00308     
00309     
00311 
00312     
00314     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00315     
00317     int getNumBlocks() const { return numBlocks_; }
00318     
00320     void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
00321     
00323     int getRecycledBlocks() const { return recycledBlocks_; }
00324 
00326     void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
00327     
00329     int getBlockSize() const { return 1; }
00330     
00332     void setBlockSize(int blockSize) {
00333       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00334        "Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
00335     }
00336 
00338     void setSize( int recycledBlocks, int numBlocks );
00339 
00341     bool isInitialized() { return initialized_; }
00342     
00344     
00345   private:
00346     
00347     //
00348     // Internal methods
00349     //
00351     void setStateSize();
00352     
00353     //
00354     // Classes inputed through constructor that define the linear problem to be solved.
00355     //
00356     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00357     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00358     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00359     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00360     
00361     //
00362     // Algorithmic parameters
00363     //  
00364     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00365     int numBlocks_; 
00366 
00367     // recycledBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
00368     int recycledBlocks_; 
00369     
00370     // Storage for QR factorization of the least squares system.
00371     Teuchos::SerialDenseVector<int,ScalarType> sn;
00372     Teuchos::SerialDenseVector<int,MagnitudeType> cs;
00373     
00374     // 
00375     // Current solver state
00376     //
00377     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00378     // is capable of running; _initialize is controlled  by the initialize() member method
00379     // For the implications of the state of initialized_, please see documentation for initialize()
00380     bool initialized_;
00381     
00382     // stateStorageInitialized_ specified that the state storage has be initialized to the current
00383     // numBlocks_ and numRecycledBlocks_.  This initialization may be postponed if the linear problem was
00384     // generated without the right-hand side or solution vectors.
00385     bool stateStorageInitialized_;
00386     
00387     // Current subspace dimension, and number of iterations performed.
00388     int curDim_, iter_;
00389     
00390     // 
00391     // State Storage
00392     //
00393     // Krylov vectors.
00394     Teuchos::RCP<MV> V_;
00395     //
00396     // Recycled subspace vectors.
00397     Teuchos::RCP<const MV> U_, C_;
00398     //
00399     // Projected matrices
00400     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00401     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00402     //
00403     // B_ : Projected matrix from the recycled subspace B = C^H*A*V
00404     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00405     //
00406     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00407     // R_: Upper triangular reduction of H
00408     // z_: Q applied to right-hand side of the least squares system
00409     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00410     Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > z_;  
00411   };
00412   
00414   // Constructor.
00415   template<class ScalarType, class MV, class OP>
00416   GCRODRIter<ScalarType,MV,OP>::GCRODRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00417              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00418              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00419              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00420              Teuchos::ParameterList &params ):
00421     lp_(problem),
00422     om_(printer),
00423     stest_(tester),
00424     ortho_(ortho),
00425     numBlocks_(0),
00426     recycledBlocks_(0),
00427     initialized_(false),
00428     stateStorageInitialized_(false),
00429     curDim_(0),
00430     iter_(0)
00431   {
00432     // Get the maximum number of blocks allowed for this Krylov subspace
00433     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00434                        "Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00435     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00436 
00437     TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,
00438                        "Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00439     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00440 
00441     // Set the number of blocks and allocate data
00442     setSize( rb, nb );
00443   }
00444   
00446   // Set the block size and make necessary adjustments.
00447   template <class ScalarType, class MV, class OP>
00448   void GCRODRIter<ScalarType,MV,OP>::setSize( int recycledBlocks, int numBlocks )
00449   {
00450     // This routine only allocates space; it doesn't not perform any computation
00451     // any change in size will invalidate the state of the solver.
00452 
00453     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::GCRODRIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
00454     TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument, "Belos::GCRODRIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
00455     TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks, std::invalid_argument, "Belos::GCRODRIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
00456 
00457     if (numBlocks == numBlocks_ && recycledBlocks == recycledBlocks_) {
00458       // do nothing
00459       return;
00460     }
00461     else {
00462       stateStorageInitialized_ = false;
00463     }
00464 
00465     numBlocks_ = numBlocks;
00466     recycledBlocks_ = recycledBlocks;
00467 
00468     initialized_ = false;
00469     curDim_ = 0;
00470 
00471     // Use the current numBlocks_ to initialize the state storage.    
00472     setStateSize();
00473     
00474   }
00475 
00477   // Setup the state storage.
00478   template <class ScalarType, class MV, class OP>
00479   void GCRODRIter<ScalarType,MV,OP>::setStateSize ()
00480   {
00481     if (!stateStorageInitialized_) {
00482 
00483       // Check if there is any multivector to clone from.
00484       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00485       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00486       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00487   stateStorageInitialized_ = false;
00488   return;
00489       }
00490       else {
00491   
00493   // blockSize*numBlocks dependent
00494   //
00495   int newsd = numBlocks_ - recycledBlocks_ + 1;
00496   
00497         cs.resize( newsd );
00498   sn.resize( newsd );
00499   
00500   // Initialize the state storage
00501         TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
00502                            "Belos::GCRODRIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
00503 
00504   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00505   if (V_ == Teuchos::null) {
00506     // Get the multivector that is not null.
00507     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00508     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00509            "Belos::GCRODRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00510     V_ = MVT::Clone( *tmp, newsd );
00511   }
00512   else {
00513     // Generate V_ by cloning itself ONLY if more space is needed.
00514     if (MVT::GetNumberVecs(*V_) < newsd) {
00515       Teuchos::RCP<const MV> tmp = V_;
00516       V_ = MVT::Clone( *tmp, newsd );
00517     }
00518   }
00519 
00520   // Generate B_ only if it doesn't exist, otherwise resize it.
00521   if (B_ == Teuchos::null)
00522     B_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_, newsd-1 ) );
00523   else
00524     B_->shapeUninitialized( recycledBlocks_, newsd-1 );
00525   
00526   // Generate H_ only if it doesn't exist, otherwise resize it.
00527   if (H_ == Teuchos::null)
00528     H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-1 ) );    
00529         else
00530     H_->shape( newsd, newsd-1 );
00531   
00532   // Generate R_ only if it doesn't exist, otherwise resize it.
00533   if (R_ == Teuchos::null)
00534     R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-1 ) );
00535   else
00536     R_->shape( newsd, newsd-1 );
00537 
00538   if (z_ == Teuchos::null)
00539     z_ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(newsd) );
00540   else
00541     z_->shapeUninitialized( newsd, 1 );
00542   
00543   // State storage has now been initialized.
00544   stateStorageInitialized_ = true;
00545       }
00546     }
00547   }
00548 
00550   // Get the current update from this subspace.
00551   template <class ScalarType, class MV, class OP>
00552   Teuchos::RCP<MV> GCRODRIter<ScalarType,MV,OP>::getCurrentUpdate() const
00553   {
00554     //
00555     // If this is the first iteration of the Arnoldi factorization, 
00556     // there is no update, so return Teuchos::null. 
00557     //
00558     RCP<MV> currentUpdate = Teuchos::null;
00559     if (curDim_==0) { 
00560       return currentUpdate; 
00561     } else {
00562       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00563       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00564       Teuchos::BLAS<int,ScalarType> blas;
00565       currentUpdate = MVT::Clone( *V_, 1 );
00566       //
00567       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00568       //
00569       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, 1 );
00570       //
00571       //  Solve the least squares problem.
00572       //
00573       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00574      Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00575      R_->values(), R_->stride(), y.values(), y.stride() );
00576       //
00577       //  Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
00578       //
00579       std::vector<int> index(curDim_);
00580       for ( int i=0; i<curDim_; i++ ) {   
00581         index[i] = i;
00582       }
00583       RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
00584       MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
00585       //
00586       //  Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
00587       //
00588       Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
00589       Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
00590       z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
00591       MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
00592     }
00593     return currentUpdate;
00594   }
00595 
00596 
00598   // Get the native residuals stored in this iteration.  
00599   // Note:  No residual std::vector will be returned by Gmres.
00600   template <class ScalarType, class MV, class OP>
00601   Teuchos::RCP<const MV> GCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00602   {
00603     //
00604     // NOTE: Make sure the incoming std::vector is the correct size!
00605     //
00606     if ( norms && (int)norms->size()==0 )                         
00607       norms->resize( 1 );                                          
00608     
00609     if (norms) {
00610       Teuchos::BLAS<int,ScalarType> blas;
00611       (*norms)[0] = blas.NRM2( 1, &(*z_)(curDim_), 1);
00612     }
00613     return Teuchos::null;
00614   }
00615   
00616   
00617 
00619   // Initialize this iteration object
00620   template <class ScalarType, class MV, class OP>
00621   void GCRODRIter<ScalarType,MV,OP>::initialize(GCRODRIterState<ScalarType,MV> newstate)
00622   {
00623     // Initialize the state storage if it isn't already.
00624     if (!stateStorageInitialized_) 
00625       setStateSize();
00626 
00627     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00628            "Belos::GCRODRIter::initialize(): Cannot initialize state storage!");
00629     
00630     // NOTE:  In GCRODRIter, V and z are required!!!  
00631     // inconsitent multivectors widths and lengths will not be tolerated, and
00632     // will be treated with exceptions.
00633     //
00634     std::string errstr("Belos::GCRODRIter::initialize(): Specified multivectors must have a consistent length and width.");
00635 
00636     if (newstate.V != Teuchos::null && newstate.U != Teuchos::null && newstate.C != Teuchos::null && newstate.z != Teuchos::null) {
00637 
00638       // initialize V_,z_, and curDim_
00639 
00640       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00641                           std::invalid_argument, errstr );
00642       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1,
00643                           std::invalid_argument, errstr );
00644       TEST_FOR_EXCEPTION( newstate.curDim > 1*(numBlocks_+1),
00645                           std::invalid_argument, errstr );
00646 
00647       curDim_ = newstate.curDim;
00648       int lclDim = MVT::GetNumberVecs(*newstate.V);
00649 
00650       // check size of Z
00651       TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
00652       
00653 
00654       // copy basis vectors from newstate into V
00655       if (newstate.V != V_) {
00656         // only copy over the first block and print a warning.
00657         if (curDim_ == 0 && lclDim > 1) {
00658     om_->stream(Warnings) << "Belos::GCRODRIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00659         << "The last " << lclDim - 1 << " vectors will be discarded." << std::endl;
00660   }
00661         std::vector<int> nevind(curDim_+1);
00662         for (int i=0; i<curDim_+1; i++) nevind[i] = i;
00663   Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
00664   Teuchos::RCP<MV> lclV = MVT::CloneView( *V_, nevind );
00665         MVT::MvAddMv( 1.0, *newV, 0.0, *newV, *lclV );
00666 
00667         // done with local pointers
00668         lclV = Teuchos::null;
00669       }
00670 
00671       // put data into z_, make sure old information is not still hanging around.
00672       if (newstate.z != z_) {
00673         z_->putScalar();
00674         Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+1,1);
00675         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclZ;
00676         lclZ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+1,1) );
00677         lclZ->assign(newZ);
00678 
00679         // done with local pointers
00680         lclZ = Teuchos::null;
00681       }
00682 
00683       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != recycledBlocks_,
00684                           std::invalid_argument, errstr );
00685       if (newstate.U != U_) {
00686   U_ = newstate.U;
00687       }
00688 
00689       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != recycledBlocks_,
00690                           std::invalid_argument, errstr );
00691       if (newstate.C != C_) {
00692   C_ = newstate.C;
00693       }
00694     }
00695     else {
00696 
00697       TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
00698                          "Belos::GCRODRIter::initialize(): GCRODRStateIterState does not have initial kernel V_0.");
00699 
00700       TEST_FOR_EXCEPTION(newstate.U == Teuchos::null,std::invalid_argument,
00701                          "Belos::GCRODRIter::initialize(): GCRODRStateIterState does not have recycled basis U.");
00702 
00703       TEST_FOR_EXCEPTION(newstate.C == Teuchos::null,std::invalid_argument,
00704                          "Belos::GCRODRIter::initialize(): GCRODRStateIterState does not have recycled basis C = A*U.");
00705 
00706       TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00707                          "Belos::GCRODRIter::initialize(): GCRODRStateIterState does not have initial norms z_0.");
00708     }
00709 
00710     // the solver is initialized
00711     initialized_ = true;
00712 
00713     /*
00714     if (om_->isVerbosity( Debug ) ) {
00715       // Check almost everything here
00716       CheckList chk;
00717       chk.checkV = true;
00718       chk.checkArn = true;
00719       chk.checkAux = true;
00720       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00721     }
00722     */
00723 
00724   }
00725 
00726 
00728   // Iterate until the status test informs us we should stop.
00729   template <class ScalarType, class MV, class OP>
00730   void GCRODRIter<ScalarType,MV,OP>::iterate()
00731   {
00732     //
00733     // Allocate/initialize data structures
00734     //
00735     if (initialized_ == false) {
00736       initialize();
00737     }
00738     
00739     // Compute the current search dimension. 
00740     int searchDim = numBlocks_ - recycledBlocks_;
00741 
00743     // iterate until the status test tells us to stop.
00744     //
00745     // also break if our basis is full
00746     //
00747     while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {
00748 
00749       iter_++;
00750 
00751       // F can be found at the curDim_ block, but the next block is at curDim_ + 1.
00752       int lclDim = curDim_ + 1; 
00753 
00754       // Get the current part of the basis.
00755       std::vector<int> curind(1);
00756       curind[0] = lclDim;
00757       Teuchos::RCP<MV> Vnext = MVT::CloneView(*V_,curind);
00758 
00759       // Get a view of the previous vectors.
00760       // This is used for orthogonalization and for computing V^H K H
00761       curind[0] = curDim_;
00762       Teuchos::RCP<MV> Vprev = MVT::CloneView(*V_,curind);
00763 
00764       // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
00765       lp_->apply(*Vprev,*Vnext);
00766       Vprev = Teuchos::null;
00767 
00768       // First, remove the recycled subspace (C) from Vnext and put coefficients in B.
00769       Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
00770       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00771   subB = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00772            ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
00773       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
00774       AsubB.append( subB );
00775 
00776       // Project out the recycled subspace.
00777       ortho_->project( *Vnext, AsubB, C );
00778 
00779       // Now, remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.          
00780       // Get a view of all the previous vectors
00781       std::vector<int> prevind(lclDim);
00782       for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00783       Vprev = MVT::CloneView(*V_,prevind);
00784       Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00785   
00786       // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00787       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00788   subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00789            ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
00790       Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00791       AsubH.append( subH );
00792       
00793       // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00794       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00795   subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00796            ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
00797 
00798       // Project out the previous Krylov vectors and normalize the next vector.
00799       int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00800 
00801       // Copy over the coefficients to R just in case we run into an error.
00802       Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,*R_,lclDim+1,1,0,curDim_ );
00803       Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
00804       subR2.assign(subH2);
00805       
00806       TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure,
00807        "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00808       //
00809       // V has been extended, and H has been extended. 
00810       //
00811       // Update the QR factorization of the upper Hessenberg matrix
00812       //
00813       updateLSQR();
00814       //
00815       // Update basis dim and release all pointers.
00816       //
00817       Vnext = Teuchos::null;
00818       curDim_++;
00819       //        
00820       /*      
00821       // When required, monitor some orthogonalities
00822       if (om_->isVerbosity( Debug ) ) {
00823       // Check almost everything here
00824       CheckList chk;
00825       chk.checkV = true;
00826       chk.checkArn = true;
00827       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00828       }
00829       else if (om_->isVerbosity( OrthoDetails ) ) {
00830         CheckList chk;
00831         chk.checkV = true;
00832         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00833       }
00834       */ 
00835       
00836     } // end while (statusTest == false)
00837    
00838   }
00839 
00840   
00841   template<class ScalarType, class MV, class OP>
00842   void GCRODRIter<ScalarType,MV,OP>::updateLSQR( int dim )
00843   {
00844     int i;
00845     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00846     
00847     // Get correct dimension based on input "dim"
00848     // Remember that ortho failures result in an exit before updateLSQR() is called.
00849     // Therefore, it is possible that dim == curDim_.
00850     int curDim = curDim_;
00851     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00852       curDim = dim;
00853     }
00854     
00855     Teuchos::LAPACK<int, ScalarType> lapack;
00856     Teuchos::BLAS<int, ScalarType> blas;
00857     //
00858     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00859     // system to upper-triangular form.
00860     //
00861     // QR factorization of Least-Squares system with Givens rotations
00862     //
00863     for (i=0; i<curDim; i++) {
00864       //
00865       // Apply previous Givens rotations to new column of Hessenberg matrix
00866       //
00867       blas.ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
00868     }
00869     //
00870     // Calculate new Givens rotation
00871     //
00872     blas.ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
00873     (*R_)(curDim+1,curDim) = zero;
00874     //
00875     // Update RHS w/ new transformation
00876     //
00877     blas.ROT( 1, &(*z_)(curDim), 1, &(*z_)(curDim+1), 1, &cs[curDim], &sn[curDim] );
00878 
00879   } // end updateLSQR()
00880 
00881 } // end Belos namespace
00882 
00883 #endif /* BELOS_GCRODR_ITER_HPP */

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7