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     
00301     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00302     
00304 
00309     Teuchos::RCP<MV> getCurrentUpdate() const;
00310     
00312 
00315     void updateLSQR( int dim = -1 );
00316     
00318     int getCurSubspaceDim() const { 
00319       if (!initialized_) return 0;
00320       return curDim_;
00321     };
00322     
00324     int getMaxSubspaceDim() const { return numBlocks_; }
00325     
00327     
00328     
00330 
00331     
00333     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00334     
00336     int getBlockSize() const { return 1; }
00337     
00339     void setBlockSize(int blockSize) { 
00340       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00341        "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
00342     }
00343     
00345     int getNumBlocks() const { return numBlocks_; }
00346     
00348     void setNumBlocks(int numBlocks);
00349     
00351     bool isInitialized() { return initialized_; }
00352     
00354     
00355   private:
00356     
00357     //
00358     // Classes inputed through constructor that define the linear problem to be solved.
00359     //
00360     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00361     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00362     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00363     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00364     
00365     //
00366     // Algorithmic parameters
00367     //  
00368     // numRHS_ is the current number of linear systems being solved.
00369     int numRHS_;
00370     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00371     int numBlocks_; 
00372     
00373     // Storage for QR factorization of the least squares system.
00374     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
00375     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
00376     
00377     // Pointers to a work vector used to improve aggregate performance.
00378     Teuchos::RCP<MV> U_vec_, AU_vec_;    
00379 
00380     // Pointers to the current right-hand side and solution multivecs being solved for.
00381     Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
00382 
00383     // 
00384     // Current solver state
00385     //
00386     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00387     // is capable of running; _initialize is controlled  by the initialize() member method
00388     // For the implications of the state of initialized_, please see documentation for initialize()
00389     bool initialized_;
00390     
00391     // Current subspace dimension, and number of iterations performed.
00392     int curDim_, iter_;
00393 
00394     // 
00395     // State Storage
00396     //
00397     std::vector<Teuchos::RCP<MV> > V_;
00398     //
00399     // Projected matrices
00400     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00401     //
00402     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
00403     // 
00404     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00405     // R_: Upper triangular reduction of H
00406     // Z_: Q applied to right-hand side of the least squares system
00407     std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
00408     std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;  
00409   };
00410   
00412   // Constructor.
00413   template<class ScalarType, class MV, class OP>
00414   PseudoBlockGmresIter<ScalarType,MV,OP>::PseudoBlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00415                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00416                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00417                      const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00418                      Teuchos::ParameterList &params ):
00419     lp_(problem),
00420     om_(printer),
00421     stest_(tester),
00422     ortho_(ortho),
00423     numRHS_(0),
00424     numBlocks_(0),
00425     initialized_(false),
00426     curDim_(0),
00427     iter_(0)
00428   {
00429     // Get the maximum number of blocks allowed for each Krylov subspace
00430     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00431                        "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
00432     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00433 
00434     setNumBlocks( nb );
00435   }
00436   
00438   // Set the block size and make necessary adjustments.
00439   template <class ScalarType, class MV, class OP>
00440   void PseudoBlockGmresIter<ScalarType,MV,OP>::setNumBlocks (int numBlocks)
00441   {
00442     // This routine only allocates space; it doesn't not perform any computation
00443     // any change in size will invalidate the state of the solver.
00444     
00445     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
00446 
00447     numBlocks_ = numBlocks;
00448     curDim_ = 0;
00449 
00450     initialized_ = false;
00451   }
00452 
00454   // Get the current update from this subspace.
00455   template <class ScalarType, class MV, class OP>
00456   Teuchos::RCP<MV> PseudoBlockGmresIter<ScalarType,MV,OP>::getCurrentUpdate() const
00457   {
00458     //
00459     // If this is the first iteration of the Arnoldi factorization, 
00460     // there is no update, so return Teuchos::null. 
00461     //
00462     Teuchos::RCP<MV> currentUpdate = Teuchos::null;
00463     if (curDim_==0) { 
00464       return currentUpdate; 
00465     } else {
00466       currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
00467       std::vector<int> index(1), index2(curDim_);
00468       for (int i=0; i<curDim_; ++i) {
00469         index2[i] = i;
00470       }
00471       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00472       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00473       Teuchos::BLAS<int,ScalarType> blas;
00474       
00475       for (int i=0; i<numRHS_; ++i) {
00476         index[0] = i;
00477         Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
00478         //
00479         //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00480         //
00481         Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
00482         //
00483         //  Solve the least squares problem and compute current solutions.
00484         //
00485         blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00486              Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00487        H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
00488   
00489   Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
00490   MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
00491       }
00492     }
00493     return currentUpdate;
00494   }
00495   
00496 
00498   // Get the native residuals stored in this iteration.  
00499   // Note:  No residual vector will be returned by Gmres.
00500   template <class ScalarType, class MV, class OP>
00501   Teuchos::RCP<const MV> 
00502   PseudoBlockGmresIter<ScalarType,MV,OP>::
00503   getNativeResiduals (std::vector<MagnitudeType> *norms) const 
00504   {
00505     typedef typename Teuchos::ScalarTraits<ScalarType> STS;
00506 
00507     if (norms)
00508       { // Resize the incoming std::vector if necessary.  The type
00509   // cast avoids the compiler warning resulting from a signed /
00510   // unsigned integer comparison.
00511   if (static_cast<int> (norms->size()) < numRHS_)
00512     norms->resize (numRHS_); 
00513 
00514   Teuchos::BLAS<int, ScalarType> blas;
00515   for (int j = 0; j < numRHS_; ++j) 
00516     {
00517       const ScalarType curNativeResid = (*Z_[j])(curDim_);
00518       (*norms)[j] = STS::magnitude (curNativeResid);
00519     }
00520     }
00521     return Teuchos::null;
00522   }
00523 
00524   
00525   template <class ScalarType, class MV, class OP>
00526   void 
00527   PseudoBlockGmresIter<ScalarType,MV,OP>::
00528   initialize (PseudoBlockGmresIterState<ScalarType,MV> newstate)
00529   {
00530     using Teuchos::RCP;
00531 
00532     // (Re)set the number of right-hand sides, by interrogating the
00533     // current LinearProblem to solve.
00534     this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
00535 
00536     // NOTE:  In PseudoBlockGmresIter, V and Z are required!!!  
00537     // Inconsistent multivectors widths and lengths will not be tolerated, and
00538     // will be treated with exceptions.
00539     //
00540     std::string errstr ("Belos::PseudoBlockGmresIter::initialize(): "
00541       "Specified multivectors must have a consistent "
00542       "length and width.");
00543 
00544     // Check that newstate has V and Z arrays with nonzero length.
00545     TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0, 
00546            std::invalid_argument,
00547                        "Belos::PseudoBlockGmresIter::initialize(): "
00548            "V and/or Z was not specified in the input state; "
00549            "the V and/or Z arrays have length zero.");
00550 
00551     // In order to create basis multivectors, we have to clone them
00552     // from some already existing multivector.  We require that at
00553     // least one of the right-hand side B and left-hand side X in the
00554     // LinearProblem be non-null.  Thus, we can clone from either B or
00555     // X.  We prefer to close from B, since B is in the range of the
00556     // operator A and the basis vectors should also be in the range of
00557     // A (the first basis vector is a scaled residual vector).
00558     // However, if B is not given, we will try our best by cloning
00559     // from X.
00560     RCP<const MV> lhsMV = lp_->getLHS();
00561     RCP<const MV> rhsMV = lp_->getRHS();
00562 
00563     // If the right-hand side is null, we make do with the left-hand
00564     // side, otherwise we use the right-hand side.
00565     RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
00566     //RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00567 
00568     TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(), 
00569            std::invalid_argument,
00570                        "Belos::PseudoBlockGmresIter::initialize(): "
00571            "The linear problem to solve does not specify multi"
00572            "vectors from which we can clone basis vectors.  The "
00573            "right-hand side(s), left-hand side(s), or both should "
00574            "be nonnull.");
00575 
00576     // Check the new dimension is not more that the maximum number of
00577     // allowable blocks.
00578     TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1, 
00579            std::invalid_argument, 
00580            errstr);
00581     curDim_ = newstate.curDim;
00582 
00583     // Initialize the state storage.  If the subspace has not be
00584     // initialized before, generate it using the right-hand side or
00585     // left-hand side from the LinearProblem lp_ to solve.
00586     V_.resize(numRHS_);
00587     for (int i=0; i<numRHS_; ++i) {
00588       // Create a new vector if we need to.  We "need to" if the
00589       // current vector V_[i] is null, or if it doesn't have enough
00590       // columns.
00591       if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
00592         V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
00593       }
00594       // Check that the newstate vector newstate.V[i] has dimensions
00595       // consistent with those of V_[i].
00596       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V[i]) != MVT::GetVecLength(*V_[i]),
00597                           std::invalid_argument, errstr );
00598       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
00599                           std::invalid_argument, errstr );
00600       //
00601       // If newstate.V[i] and V_[i] are not identically the same
00602       // vector, then copy newstate.V[i] into V_[i].  
00603       //
00604       int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
00605       if (newstate.V[i] != V_[i]) {
00606         // Only copy over the first block and print a warning.
00607         if (curDim_ == 0 && lclDim > 1) {
00608           om_->stream(Warnings) 
00609       << "Belos::PseudoBlockGmresIter::initialize(): the solver was "
00610       << "initialized with a kernel of " << lclDim 
00611       << std::endl 
00612       << "The block size however is only " << 1 
00613       << std::endl
00614       << "The last " << lclDim - 1 << " vectors will be discarded." 
00615       << std::endl;
00616         }
00617         std::vector<int> nevind (curDim_ + 1);
00618         for (int j = 0; j < curDim_ + 1; ++j) 
00619     nevind[j] = j;
00620 
00621         RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
00622         RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
00623         const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00624         const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00625         MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
00626 
00627         // Done with local pointers
00628         lclV = Teuchos::null;
00629       }
00630     }
00631 
00632 
00633     // Check size of Z
00634     Z_.resize(numRHS_);
00635     for (int i=0; i<numRHS_; ++i) {
00636       // Create a vector if we need to.
00637       if (Z_[i] == Teuchos::null) {
00638   Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
00639       }
00640       if (Z_[i]->length() < numBlocks_+1) {
00641   Z_[i]->shapeUninitialized(numBlocks_+1, 1); 
00642       }
00643       
00644       // Check that the newstate vector is consistent.
00645       TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
00646       
00647       // Put data into Z_, make sure old information is not still hanging around.
00648       if (newstate.Z[i] != Z_[i]) {
00649   if (curDim_==0)
00650     Z_[i]->putScalar();
00651   
00652         Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
00653         Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
00654         lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
00655         lclZ->assign(newZ);
00656   
00657         // Done with local pointers
00658   lclZ = Teuchos::null;
00659       }
00660     }
00661 
00662 
00663     // Check size of H
00664     H_.resize(numRHS_);
00665     for (int i=0; i<numRHS_; ++i) {
00666       // Create a matrix if we need to.
00667       if (H_[i] == Teuchos::null) {
00668   H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
00669       }
00670       if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
00671   H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
00672       }
00673       
00674       // Put data into H_ if it exists, make sure old information is not still hanging around.
00675       if ((int)newstate.H.size() == numRHS_) {
00676   
00677   // Check that the newstate matrix is consistent.
00678   TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument, 
00679          "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
00680   
00681   if (newstate.H[i] != H_[i]) {
00682     // H_[i]->putScalar();
00683     
00684     Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
00685     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00686     lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
00687     lclH->assign(newH);
00688     
00689     // Done with local pointers
00690     lclH = Teuchos::null;
00691   }
00692       }
00693     }      
00694     
00696     // Reinitialize storage for least squares solve
00697     //
00698     cs_.resize(numRHS_);
00699     sn_.resize(numRHS_);
00700       
00701     // Copy over rotation angles if they exist
00702     if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
00703       for (int i=0; i<numRHS_; ++i) {
00704         if (cs_[i] != newstate.cs[i])
00705           cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
00706         if (sn_[i] != newstate.sn[i])
00707           sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
00708       }
00709     } 
00710       
00711     // Resize or create the vectors as necessary
00712     for (int i=0; i<numRHS_; ++i) {
00713       if (cs_[i] == Teuchos::null) 
00714         cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
00715       else
00716         cs_[i]->resize(numBlocks_+1);
00717       if (sn_[i] == Teuchos::null)
00718         sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
00719       else
00720         sn_[i]->resize(numBlocks_+1);
00721     }
00722 
00723     // the solver is initialized
00724     initialized_ = true;
00725       
00726       /*
00727   if (om_->isVerbosity( Debug ) ) {
00728       // Check almost everything here
00729       CheckList chk;
00730       chk.checkV = true;
00731       chk.checkArn = true;
00732       chk.checkAux = true;
00733       om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00734     }
00735     */
00736 
00737   }
00738 
00739 
00741   // Iterate until the status test informs us we should stop.
00742   template <class ScalarType, class MV, class OP>
00743   void PseudoBlockGmresIter<ScalarType,MV,OP>::iterate()
00744   {
00745     //
00746     // Allocate/initialize data structures
00747     //
00748     if (initialized_ == false) {
00749       initialize();
00750     }
00751 
00752     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00753     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00754     
00755     // Compute the current search dimension. 
00756     int searchDim = numBlocks_;
00757     //
00758     // Associate each initial block of V_[i] with U_vec[i]
00759     // Reset the index vector (this might have been changed if there was a restart)
00760     //
00761     std::vector<int> index(1);
00762     std::vector<int> index2(1);
00763     index[0] = curDim_;
00764     Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
00765 
00766     // Create AU_vec to hold A*U_vec.
00767     Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
00768 
00769     for (int i=0; i<numRHS_; ++i) {
00770       index2[0] = i;
00771       Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
00772       Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
00773       MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
00774     }
00775     
00777     // iterate until the status test tells us to stop.
00778     //
00779     // also break if our basis is full
00780     //
00781     while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
00782 
00783       iter_++;
00784       //
00785       // Apply the operator to _work_vector
00786       //
00787       lp_->apply( *U_vec, *AU_vec );
00788       //
00789       //
00790       // Resize index.
00791       //
00792       int num_prev = curDim_+1;
00793       index.resize( num_prev );
00794       for (int i=0; i<num_prev; ++i) { 
00795   index[i] = i; 
00796       }
00797       //
00798       // Orthogonalize next Krylov vector for each right-hand side.
00799       //
00800       for (int i=0; i<numRHS_; ++i) {
00801   //
00802   // Get previous Krylov vectors.
00803   //
00804   Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
00805   Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
00806   //
00807   // Get a view of the new candidate std::vector.
00808   //
00809   index2[0] = i;
00810   Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
00811   //
00812   // Get a view of the current part of the upper-hessenberg matrix.
00813   //
00814   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new 
00815     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00816         ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
00817   Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
00818   
00819   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
00820     = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
00821         ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
00822   //
00823   // Orthonormalize the new block of the Krylov expansion
00824   // NOTE:  Rank deficiencies are not checked because this is a single-std::vector Krylov method.
00825   //
00826   ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
00827   //
00828   // NOTE:  V_new is a copy of the iter+1 vector in V_[i], so the normalized vector has to be
00829   // be copied back in when V_new is changed.  
00830   //
00831   index2[0] = curDim_+1;
00832   Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
00833   MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
00834       }
00835       // 
00836       // Now _AU_vec is the new _U_vec, so swap these two vectors.
00837       // NOTE: This alleviates the need for allocating a vector for AU_vec each iteration.
00838       // 
00839       Teuchos::RCP<MV> tmp_AU_vec = U_vec;
00840       U_vec = AU_vec;
00841       AU_vec = tmp_AU_vec;
00842       //
00843       // V has been extended, and H has been extended. 
00844       //
00845       // Update the QR factorization of the upper Hessenberg matrix
00846       //
00847       updateLSQR();
00848       //
00849       // Update basis dim and release all pointers.
00850       //
00851       curDim_ += 1;
00852       //        
00853       /*      
00854       // When required, monitor some orthogonalities
00855       if (om_->isVerbosity( Debug ) ) {
00856       // Check almost everything here
00857       CheckList chk;
00858       chk.checkV = true;
00859       chk.checkArn = true;
00860       om_->print( Debug, accuracyCheck(chk, ": after local update") );
00861       }
00862       else if (om_->isVerbosity( OrthoDetails ) ) {
00863         CheckList chk;
00864         chk.checkV = true;
00865         om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
00866       }
00867       */ 
00868       
00869     } // end while (statusTest == false)
00870    
00871   }
00872 
00874   // Update the least squares solution for each right-hand side.
00875   template<class ScalarType, class MV, class OP>
00876   void PseudoBlockGmresIter<ScalarType,MV,OP>::updateLSQR( int dim )
00877   {
00878     // Get correct dimension based on input "dim"
00879     // Remember that ortho failures result in an exit before updateLSQR() is called.
00880     // Therefore, it is possible that dim == curDim_.
00881     int curDim = curDim_;
00882     if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
00883       curDim = dim;
00884     }
00885     
00886     int i, j;
00887     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00888 
00889     Teuchos::BLAS<int, ScalarType> blas;
00890     
00891     for (i=0; i<numRHS_; ++i) {
00892       //
00893       // Update the least-squares QR for each linear system.
00894       //
00895       // QR factorization of Least-Squares system with Givens rotations
00896       //
00897       for (j=0; j<curDim; j++) {
00898   //
00899   // Apply previous Givens rotations to new column of Hessenberg matrix
00900   //
00901   blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
00902       }
00903       //
00904       // Calculate new Givens rotation
00905       //
00906       blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00907       (*H_[i])(curDim+1,curDim) = zero;
00908       //
00909       // Update RHS w/ new transformation
00910       //
00911       blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
00912     }
00913 
00914   } // end updateLSQR()
00915 
00916 } // end Belos namespace
00917 
00918 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines