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