Belos Version of the Day
BelosPseudoBlockGmresIter.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_PSEUDO_BLOCK_GMRES_ITER_HPP
00043 #define BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosIteration.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   
00083 
00084   
00089   template <class ScalarType, class MV>
00090   struct PseudoBlockGmresIterState {
00091 
00092     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00093     typedef typename SCT::magnitudeType MagnitudeType;
00094 
00099     int curDim;
00101     std::vector<Teuchos::RCP<const MV> > V;
00107     std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > H;
00109     std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > R;
00111     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > Z;
00113     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > sn;
00114     std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > > cs;
00115 
00116     PseudoBlockGmresIterState() : curDim(0), V(0),
00117           H(0), R(0), Z(0),
00118                                   sn(0), cs(0)
00119     {}
00120   };
00121   
00123 
00124   
00137   class PseudoBlockGmresIterInitFailure : public BelosError {public:
00138     PseudoBlockGmresIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00139     {}};
00140   
00147   class PseudoBlockGmresIterOrthoFailure : public BelosError {public:
00148     PseudoBlockGmresIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00149     {}};
00150   
00152   
00153   
00154   template<class ScalarType, class MV, class OP>
00155   class PseudoBlockGmresIter : virtual public Iteration<ScalarType,MV,OP> {
00156     
00157   public:
00158     
00159     //
00160     // Convenience typedefs
00161     //
00162     typedef MultiVecTraits<ScalarType,MV> MVT;
00163     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00164     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00165     typedef typename SCT::magnitudeType MagnitudeType;
00166     
00168 
00169     
00178     PseudoBlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00179         const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00180         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00181         const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00182         Teuchos::ParameterList &params );
00183     
00185     virtual ~PseudoBlockGmresIter() {};
00187     
00188     
00190 
00191     
00213     void iterate();
00214     
00236     void initialize(PseudoBlockGmresIterState<ScalarType,MV> newstate);
00237     
00241     void initialize()
00242     {
00243       PseudoBlockGmresIterState<ScalarType,MV> empty;
00244       initialize(empty);
00245     }
00246     
00254     PseudoBlockGmresIterState<ScalarType,MV> getState() const {
00255       PseudoBlockGmresIterState<ScalarType,MV> state;
00256       state.curDim = curDim_;
00257       state.V.resize(numRHS_);
00258       state.H.resize(numRHS_);
00259       state.Z.resize(numRHS_);
00260       state.sn.resize(numRHS_);
00261       state.cs.resize(numRHS_);  
00262       for (int i=0; i<numRHS_; ++i) {
00263   state.V[i] = V_[i];
00264   state.H[i] = H_[i];
00265   state.Z[i] = Z_[i];
00266         state.sn[i] = sn_[i];
00267         state.cs[i] = cs_[i];
00268       }
00269       return state;
00270     }
00271     
00273     
00274     
00276 
00277     
00279     int getNumIters() const { return iter_; }
00280     
00282     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00283     
00286     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00287     
00289 
00294     Teuchos::RCP<MV> getCurrentUpdate() const;
00295     
00297 
00300     void updateLSQR( int dim = -1 );
00301     
00303     int getCurSubspaceDim() const { 
00304       if (!initialized_) return 0;
00305       return curDim_;
00306     };
00307     
00309     int getMaxSubspaceDim() const { return numBlocks_; }
00310     
00312     
00313     
00315 
00316     
00318     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00319     
00321     int getBlockSize() const { return 1; }
00322     
00324     void setBlockSize(int blockSize) { 
00325       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00326        "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
00327     }
00328     
00330     int getNumBlocks() const { return numBlocks_; }
00331     
00333     void setNumBlocks(int numBlocks);
00334     
00336     bool isInitialized() { return initialized_; }
00337     
00339     
00340   private:
00341     
00342     //
00343     // Classes inputed through constructor that define the linear problem to be solved.
00344     //
00345     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00346     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00347     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00348     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00349     
00350     //
00351     // Algorithmic parameters
00352     //  
00353     // numRHS_ is the current number of linear systems being solved.
00354     int numRHS_;
00355     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00356     int numBlocks_; 
00357     
00358     // Storage for QR factorization of the least squares system.
00359     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
00360     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
00361     
00362     // Pointers to a work std::vector used to improve aggregate performance.
00363     RCP<MV> U_vec_, AU_vec_;    
00364 
00365     // Pointers to the current right-hand side and solution multivecs being solved for.
00366     RCP<MV> cur_block_rhs_, cur_block_sol_;
00367 
00368     // 
00369     // Current solver state
00370     //
00371     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00372     // is capable of running; _initialize is controlled  by the initialize() member method
00373     // For the implications of the state of initialized_, please see documentation for initialize()
00374     bool initialized_;
00375     
00376     // Current subspace dimension, and number of iterations performed.
00377     int curDim_, iter_;
00378 
00379     // 
00380     // State Storage
00381     //
00382     std::vector<Teuchos::RCP<MV> > V_;
00383     //
00384     // Projected matrices
00385     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00386     //
00387     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
00388     // 
00389     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00390     // R_: Upper triangular reduction of H
00391     // Z_: Q applied to right-hand side of the least squares system
00392     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
00393     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;  
00394   };
00395   
00397   // Constructor.
00398   template<class ScalarType, class MV, class OP>
00399   PseudoBlockGmresIter<ScalarType,MV,OP>::PseudoBlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00400                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00401                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00402                      const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00403                      Teuchos::ParameterList &params ):
00404     lp_(problem),
00405     om_(printer),
00406     stest_(tester),
00407     ortho_(ortho),
00408     numRHS_(0),
00409     numBlocks_(0),
00410     initialized_(false),
00411     curDim_(0),
00412     iter_(0)
00413   {
00414     // Get the maximum number of blocks allowed for each Krylov subspace
00415     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00416                        "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00417     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00418 
00419     setNumBlocks( nb );
00420   }
00421   
00423   // Set the block size and make necessary adjustments.
00424   template <class ScalarType, class MV, class OP>
00425   void PseudoBlockGmresIter<ScalarType,MV,OP>::setNumBlocks (int numBlocks)
00426   {
00427     // This routine only allocates space; it doesn't not perform any computation
00428     // any change in size will invalidate the state of the solver.
00429     
00430     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
00431 
00432     numBlocks_ = numBlocks;
00433     curDim_ = 0;
00434 
00435     initialized_ = false;
00436   }
00437 
00439   // Get the current update from this subspace.
00440   template <class ScalarType, class MV, class OP>
00441   Teuchos::RCP<MV> PseudoBlockGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00442   {
00443     //
00444     // If this is the first iteration of the Arnoldi factorization, 
00445     // there is no update, so return Teuchos::null. 
00446     //
00447     RCP<MV> currentUpdate = Teuchos::null;
00448     if (curDim_==0) { 
00449       return currentUpdate; 
00450     } else {
00451       currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
00452       std::vector<int> index(1), index2(curDim_);
00453       for (int i=0; i<curDim_; ++i) {
00454         index2[i] = i;
00455       }
00456       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00457       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00458       Teuchos::BLAS<int,ScalarType> blas;
00459       
00460       for (int i=0; i<numRHS_; ++i) {
00461         index[0] = i;
00462         RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
00463         //
00464         //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00465         //
00466         Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
00467         //
00468         //  Solve the least squares problem and compute current solutions.
00469         //
00470         blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00471              Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00472        H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
00473   
00474   RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
00475   MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
00476       }
00477     }
00478     return currentUpdate;
00479   }
00480   
00481 
00483   // Get the native residuals stored in this iteration.  
00484   // Note:  No residual std::vector will be returned by Gmres.
00485   template <class ScalarType, class MV, class OP>
00486   Teuchos::RCP<const MV> PseudoBlockGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00487   {
00488     //
00489     // NOTE: Make sure the incoming std::vector is the correct size!
00490     //
00491     if ( norms && (int)norms->size() < numRHS_ )                         
00492       norms->resize( numRHS_ );                                          
00493 
00494     if (norms) {
00495       Teuchos::BLAS<int,ScalarType> blas;
00496       for (int j=0; j<numRHS_; j++) {
00497         (*norms)[j] = Teuchos::ScalarTraits<ScalarType>::magnitude( (*Z_[j])(curDim_) );
00498       }
00499     }
00500     return Teuchos::null;
00501   }
00502 
00503   
00504 
00506   // Initialize this iteration object
00507   template <class ScalarType, class MV, class OP>
00508   void PseudoBlockGmresIter<ScalarType,MV,OP>::initialize(PseudoBlockGmresIterState<ScalarType,MV> newstate)
00509   {
00510     // Get the number of right-hand sides we're solving for now.
00511     int numRHS = MVT::GetNumberVecs(*(lp_->getCurrLHSVec()));
00512     numRHS_ = numRHS;
00513 
00514     // NOTE:  In PseudoBlockGmresIter, V and Z are required!!!  
00515     // inconsitent multivectors widths and lengths will not be tolerated, and
00516     // will be treated with exceptions.
00517     //
00518     std::string errstr("Belos::PseudoBlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
00519 
00520     // Check that we have V and Z.
00521     TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0, std::invalid_argument,
00522                        "Belos::PseudoBlockGmresIter::initialize(): V and/or Z is not specified.");
00523 
00524     // Get the multivector that is not null.
00525     Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00526     Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00527     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00528     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00529                        "Belos::PseudoBlockGmresIter::initialize(): linear problem does not specify multivectors to clone from.");
00530 
00531     // Check the new dimension is not more that the maximum number of allowable blocks.  
00532     TEST_FOR_EXCEPTION( newstate.curDim > numBlocks_+1,
00533                         std::invalid_argument, errstr );
00534     curDim_ = newstate.curDim;
00535 
00536     // Initialize the state storage
00537     // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00538     V_.resize(numRHS_);
00539     for (int i=0; i<numRHS_; ++i) {
00540       // Create a new std::vector if we need to.
00541       if (V_[i] == Teuchos::null || MVT::GetNumberVecs(*V_[i]) < numBlocks_+1 ) {
00542         V_[i] = MVT::Clone(*tmp,numBlocks_+1);
00543       }
00544       // Check that the newstate std::vector is consistent.
00545       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V[i]) != MVT::GetVecLength(*V_[i]),
00546                           std::invalid_argument, errstr );
00547       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
00548                           std::invalid_argument, errstr );
00549       
00550       int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
00551       if (newstate.V[i] != V_[i]) {
00552         // Cnly copy over the first block and print a warning.
00553         if (curDim_ == 0 && lclDim > 1) {
00554           om_->stream(Warnings) << "Belos::PseudoBlockGmresIter::initialize(): the solver was initialized with a kernel of " 
00555                                 << lclDim << std::endl << "The block size however is only " << 1 << std::endl
00556                                 << "The last " << lclDim - 1 << " vectors will be discarded." << std::endl;
00557         }
00558         std::vector<int> nevind(curDim_+1);
00559         for (int j=0; j<curDim_+1; j++) nevind[j] = j;
00560         Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V[i], nevind );
00561         Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
00562         const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00563         const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00564         MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
00565 
00566         // Done with local pointers
00567         lclV = Teuchos::null;
00568       }
00569     }
00570 
00571 
00572     // Check size of Z
00573     Z_.resize(numRHS_);
00574     for (int i=0; i<numRHS_; ++i) {
00575       // Create a std::vector if we need to.
00576       if (Z_[i] == Teuchos::null) {
00577   Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
00578       }
00579       if (Z_[i]->length() < numBlocks_+1) {
00580   Z_[i]->shapeUninitialized(numBlocks_+1, 1); 
00581       }
00582       
00583       // Check that the newstate std::vector is consistent.
00584       TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
00585       
00586       // Put data into Z_, make sure old information is not still hanging around.
00587       if (newstate.Z[i] != Z_[i]) {
00588   if (curDim_==0)
00589     Z_[i]->putScalar();
00590   
00591         Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
00592         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
00593         lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
00594         lclZ->assign(newZ);
00595   
00596         // Done with local pointers
00597   lclZ = Teuchos::null;
00598       }
00599     }
00600 
00601 
00602     // Check size of H
00603     H_.resize(numRHS_);
00604     for (int i=0; i<numRHS_; ++i) {
00605       // Create a matrix if we need to.
00606       if (H_[i] == Teuchos::null) {
00607   H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00608       }
00609       if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
00610   H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
00611       }
00612       
00613       // Put data into H_ if it exists, make sure old information is not still hanging around.
00614       if ((int)newstate.H.size() == numRHS_) {
00615   
00616   // Check that the newstate matrix is consistent.
00617   TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument, 
00618          "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
00619   
00620   if (newstate.H[i] != H_[i]) {
00621     // H_[i]->putScalar();
00622     
00623     Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
00624     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00625     lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
00626     lclH->assign(newH);
00627     
00628     // Done with local pointers
00629     lclH = Teuchos::null;
00630   }
00631       }
00632     }      
00633     
00635     // Reinitialize storage for least squares solve
00636     //
00637     cs_.resize(numRHS_);
00638     sn_.resize(numRHS_);
00639       
00640     // Copy over rotation angles if they exist
00641     if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
00642       for (int i=0; i<numRHS_; ++i) {
00643         if (cs_[i] != newstate.cs[i])
00644           cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
00645         if (sn_[i] != newstate.sn[i])
00646           sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
00647       }
00648     } 
00649       
00650     // Resize or create the vectors as necessary
00651     for (int i=0; i<numRHS_; ++i) {
00652       if (cs_[i] == Teuchos::null) 
00653         cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
00654       else
00655         cs_[i]->resize(numBlocks_+1);
00656       if (sn_[i] == Teuchos::null)
00657         sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
00658       else
00659         sn_[i]->resize(numBlocks_+1);
00660     }
00661 
00662     // the solver is initialized
00663     initialized_ = true;
00664       
00665       /*
00666   if (om_->isVerbosity( Debug ) ) {
00667       // Check almost everything here
00668       CheckList chk;
00669       chk.checkV = true;
00670       chk.checkArn = true;
00671       chk.checkAux = true;
00672       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00673     }
00674     */
00675 
00676   }
00677 
00678 
00680   // Iterate until the status test informs us we should stop.
00681   template <class ScalarType, class MV, class OP>
00682   void PseudoBlockGmresIter<ScalarType,MV,OP>::iterate()
00683   {
00684     //
00685     // Allocate/initialize data structures
00686     //
00687     if (initialized_ == false) {
00688       initialize();
00689     }
00690 
00691     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00692     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00693     
00694     // Compute the current search dimension. 
00695     int searchDim = numBlocks_;
00696     //
00697     // Associate each initial block of V_[i] with U_vec[i]
00698     // Reset the index std::vector (this might have been changed if there was a restart)
00699     //
00700     std::vector<int> index(1);
00701     std::vector<int> index2(1);
00702     index[0] = curDim_;
00703     Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
00704 
00705     // Create AU_vec to hold A*U_vec.
00706     Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
00707 
00708     for (int i=0; i<numRHS_; ++i) {
00709       index2[0] = i;
00710       RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
00711       RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
00712       MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
00713     }
00714     
00716     // iterate until the status test tells us to stop.
00717     //
00718     // also break if our basis is full
00719     //
00720     while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
00721 
00722       iter_++;
00723       //
00724       // Apply the operator to _work_vector
00725       //
00726       lp_->apply( *U_vec, *AU_vec );
00727       //
00728       //
00729       // Resize index.
00730       //
00731       int num_prev = curDim_+1;
00732       index.resize( num_prev );
00733       for (int i=0; i<num_prev; ++i) { 
00734   index[i] = i; 
00735       }
00736       //
00737       // Orthogonalize next Krylov std::vector for each right-hand side.
00738       //
00739       for (int i=0; i<numRHS_; ++i) {
00740   //
00741   // Get previous Krylov vectors.
00742   //
00743   RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
00744   Teuchos::Array< RCP<const MV> > V_array( 1, V_prev );
00745   //
00746   // Get a view of the new candidate std::vector.
00747   //
00748   index2[0] = i;
00749   RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
00750   //
00751   // Get a view of the current part of the upper-hessenberg matrix.
00752   //
00753   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new 
00754     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00755         ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
00756   Teuchos::Array< RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
00757   
00758   RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
00759     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00760         ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
00761   //
00762   // Orthonormalize the new block of the Krylov expansion
00763   // NOTE:  Rank deficiencies are not checked because this is a single-std::vector Krylov method.
00764   //
00765   ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
00766   //
00767   // NOTE:  V_new is a copy of the iter+1 std::vector in V_[i], so the normalized std::vector has to be
00768   // be copied back in when V_new is changed.  
00769   //
00770   index2[0] = curDim_+1;
00771   RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
00772   MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
00773       }
00774       // 
00775       // Now _AU_vec is the new _U_vec, so swap these two vectors.
00776       // NOTE: This alleviates the need for allocating a std::vector for AU_vec each iteration.
00777       // 
00778       RCP<MV> tmp_AU_vec = U_vec;
00779       U_vec = AU_vec;
00780       AU_vec = tmp_AU_vec;
00781       //
00782       // V has been extended, and H has been extended. 
00783       //
00784       // Update the QR factorization of the upper Hessenberg matrix
00785       //
00786       updateLSQR();
00787       //
00788       // Update basis dim and release all pointers.
00789       //
00790       curDim_ += 1;
00791       //        
00792       /*      
00793       // When required, monitor some orthogonalities
00794       if (om_->isVerbosity( Debug ) ) {
00795       // Check almost everything here
00796       CheckList chk;
00797       chk.checkV = true;
00798       chk.checkArn = true;
00799       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00800       }
00801       else if (om_->isVerbosity( OrthoDetails ) ) {
00802         CheckList chk;
00803         chk.checkV = true;
00804         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00805       }
00806       */ 
00807       
00808     } // end while (statusTest == false)
00809    
00810   }
00811 
00813   // Update the least squares solution for each right-hand side.
00814   template<class ScalarType, class MV, class OP>
00815   void PseudoBlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00816   {
00817     // Get correct dimension based on input "dim"
00818     // Remember that ortho failures result in an exit before updateLSQR() is called.
00819     // Therefore, it is possible that dim == curDim_.
00820     int curDim = curDim_;
00821     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00822       curDim = dim;
00823     }
00824     
00825     int i, j;
00826     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00827 
00828     Teuchos::BLAS<int, ScalarType> blas;
00829     
00830     for (i=0; i<numRHS_; ++i) {
00831       //
00832       // Update the least-squares QR for each linear system.
00833       //
00834       // QR factorization of Least-Squares system with Givens rotations
00835       //
00836       for (j=0; j<curDim; j++) {
00837   //
00838   // Apply previous Givens rotations to new column of Hessenberg matrix
00839   //
00840   blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
00841       }
00842       //
00843       // Calculate new Givens rotation
00844       //
00845       blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00846       (*H_[i])(curDim+1,curDim) = zero;
00847       //
00848       // Update RHS w/ new transformation
00849       //
00850       blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00851     }
00852 
00853   } // end updateLSQR()
00854 
00855 } // end Belos namespace
00856 
00857 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines