BelosRCGIter.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_RCG_ITER_HPP
00030 #define BELOS_RCG_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosMatOrthoManager.hpp"
00041 #include "BelosOutputManager.hpp"
00042 #include "BelosStatusTest.hpp"
00043 #include "BelosOperatorTraits.hpp"
00044 #include "BelosMultiVecTraits.hpp"
00045 
00046 #include "Teuchos_BLAS.hpp"
00047 #include "Teuchos_LAPACK.hpp"
00048 #include "Teuchos_SerialDenseMatrix.hpp"
00049 #include "Teuchos_SerialDenseVector.hpp"
00050 #include "Teuchos_ScalarTraits.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00054 // MLP Remove after debugging
00055 #include <fstream>
00056 #include <iomanip>
00057 
00069 namespace Belos {
00070   
00072 
00073   
00078   template <class ScalarType, class MV>
00079   struct RCGIterState {
00084     int curDim;
00085     
00087     Teuchos::RCP<MV> P;
00088 
00090     Teuchos::RCP<MV> Ap;
00091 
00093     Teuchos::RCP<MV> r;
00094  
00096     Teuchos::RCP<MV> z;
00097  
00099     bool existU;
00100    
00102     Teuchos::RCP<MV> U, AU;
00103 
00106     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha;
00107     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta;
00108     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D;
00109     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old;
00110 
00112     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta;
00113 
00115     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU;
00117     Teuchos::RCP<std::vector<int> > ipiv;
00118 
00119 
00120     RCGIterState() : curDim(0), P(Teuchos::null), Ap(Teuchos::null), r(Teuchos::null), 
00121                      z(Teuchos::null),
00122                      existU(false),
00123          U(Teuchos::null), AU(Teuchos::null),
00124          Alpha(Teuchos::null), Beta(Teuchos::null), D(Teuchos::null), rTz_old(Teuchos::null),
00125          Delta(Teuchos::null), LUUTAU(Teuchos::null), ipiv(Teuchos::null)
00126     {}
00127   };
00128   
00130   
00132 
00133   
00145   class RCGIterInitFailure : public BelosError {public:
00146     RCGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00147     {}};
00148 
00155   class RCGIterFailure : public BelosError {public:
00156     RCGIterFailure(const std::string& what_arg) : BelosError(what_arg)
00157     {}};
00158   
00165   class RCGIterOrthoFailure : public BelosError {public:
00166     RCGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00167     {}};
00168   
00175   class RCGIterLAPACKFailure : public BelosError {public:
00176     RCGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00177     {}};
00178   
00180   
00181   
00182   template<class ScalarType, class MV, class OP>
00183   class RCGIter : virtual public Iteration<ScalarType,MV,OP> {
00184     
00185   public:
00186     
00187     //
00188     // Convenience typedefs
00189     //
00190     typedef MultiVecTraits<ScalarType,MV> MVT;
00191     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00192     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00193     typedef typename SCT::magnitudeType MagnitudeType;
00194     
00196 
00197     
00205     RCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00206     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00207     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00208     Teuchos::ParameterList &params );
00209     
00211     virtual ~RCGIter() {};
00213     
00214     
00216 
00217     
00230     void iterate();
00231     
00246     void initialize(RCGIterState<ScalarType,MV> &newstate);
00247  
00251     void initialize()
00252     {
00253       RCGIterState<ScalarType,MV> empty;
00254       initialize(empty);
00255     }
00256 
00258     
00260 
00261     
00263     int getNumIters() const { return iter_; }
00264     
00266     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00267     
00270     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return r_; }
00271     
00273 
00278     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00279     
00281     int getCurSubspaceDim() const { 
00282       if (!initialized_) return 0;
00283       return curDim_;
00284     };
00285     
00287     int getMaxSubspaceDim() const { return numBlocks_+1; }
00288     
00290     
00291     
00293 
00294     
00296     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00297     
00299     int getNumBlocks() const { return numBlocks_; }
00300     
00302     void setNumBlocks(int numBlocks) { setSize( recycleBlocks_, numBlocks ); };
00303     
00305     int getRecycledBlocks() const { return recycleBlocks_; }
00306 
00308     void setRecycledBlocks(int recycleBlocks) { setSize( recycleBlocks, numBlocks_ ); };
00309     
00311     int getBlockSize() const { return 1; }
00312     
00314     void setBlockSize(int blockSize) {
00315       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00316        "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
00317     }
00318 
00320     void setSize( int recycleBlocks, int numBlocks );
00321 
00323     bool isInitialized() { return initialized_; }
00324     
00326     
00327   private:
00328     
00329     //
00330     // Internal methods
00331     //
00332     
00333     //
00334     // Classes input through constructor that define the linear problem to be solved.
00335     //
00336     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00337     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00338     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00339     
00340     //
00341     // Algorithmic parameters
00342     //  
00343     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00344     int numBlocks_; 
00345 
00346     // recycleBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
00347     int recycleBlocks_; 
00348     
00349     // 
00350     // Current solver state
00351     //
00352     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00353     // is capable of running; _initialize is controlled  by the initialize() member method
00354     // For the implications of the state of initialized_, please see documentation for initialize()
00355     bool initialized_;
00356     
00357     // Current subspace dimension, and number of iterations performed.
00358     int curDim_, iter_;
00359     
00360     // 
00361     // State Storage
00362     //
00363     // Search vectors
00364     Teuchos::RCP<MV> P_;
00365     //
00366     // A times current search vector
00367     Teuchos::RCP<MV> Ap_;
00368     //
00369     // Residual vector
00370     Teuchos::RCP<MV> r_;
00371     //
00372     // Preconditioned residual
00373     Teuchos::RCP<MV> z_;
00374     //
00375     // Flag to indicate that the recycle space should be used
00376     bool existU_;
00377     // Recycled subspace and its image
00378     Teuchos::RCP<MV> U_, AU_;
00379     //
00380     // Coefficients arising in RCG iteration
00381     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_,Beta_,D_;
00382     //
00383     // Solutions to local least-squares problems
00384     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
00385     //
00386     // The LU factorization of the matrix U^T A U
00387     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
00388     //
00389     // Data from LU factorization of UTAU
00390     Teuchos::RCP<std::vector<int> > ipiv_;
00391     //
00392     // The scalar r'*z
00393     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
00394   };
00395   
00397   // Constructor.
00398   template<class ScalarType, class MV, class OP>
00399   RCGIter<ScalarType,MV,OP>::RCGIter(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              Teuchos::ParameterList &params ):
00403     lp_(problem),
00404     om_(printer),
00405     stest_(tester),
00406     numBlocks_(0),
00407     recycleBlocks_(0),
00408     initialized_(false),
00409     curDim_(0),
00410     iter_(0)
00411   {
00412     // Get the maximum number of blocks allowed for this Krylov subspace
00413     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00414                        "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00415     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00416 
00417     TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,
00418                        "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00419     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00420 
00421     // Set the number of blocks and allocate data
00422     setSize( rb, nb );
00423   }
00424   
00426   // Set the block size and make necessary adjustments.
00427   template <class ScalarType, class MV, class OP>
00428   void RCGIter<ScalarType,MV,OP>::setSize( int recycleBlocks, int numBlocks )
00429   {
00430 
00431     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
00432     TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
00433     TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument, "Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
00434 
00435     numBlocks_ = numBlocks;
00436     recycleBlocks_ = recycleBlocks;
00437 
00438   }
00439 
00441   // Initialize this iteration object
00442   template <class ScalarType, class MV, class OP>
00443   void RCGIter<ScalarType,MV,OP>::initialize(RCGIterState<ScalarType,MV> &newstate)
00444   {
00445 
00446     if (newstate.P != Teuchos::null && 
00447         newstate.Ap != Teuchos::null && 
00448         newstate.r != Teuchos::null && 
00449         newstate.z != Teuchos::null && 
00450         newstate.U != Teuchos::null && 
00451         newstate.AU != Teuchos::null && 
00452         newstate.Alpha != Teuchos::null && 
00453         newstate.Beta != Teuchos::null && 
00454         newstate.D != Teuchos::null && 
00455         newstate.Delta != Teuchos::null && 
00456         newstate.LUUTAU != Teuchos::null &&
00457         newstate.ipiv != Teuchos::null &&
00458         newstate.rTz_old != Teuchos::null) {
00459 
00460       curDim_ = newstate.curDim;
00461       P_ = newstate.P;
00462       Ap_ = newstate.Ap;
00463       r_ = newstate.r;
00464       z_ = newstate.z;
00465       existU_ = newstate.existU;
00466       U_ = newstate.U;
00467       AU_ = newstate.AU;
00468       Alpha_ = newstate.Alpha;
00469       Beta_ = newstate.Beta;
00470       D_ = newstate.D;
00471       Delta_ = newstate.Delta;
00472       LUUTAU_ = newstate.LUUTAU;
00473       ipiv_ = newstate.ipiv;
00474       rTz_old_ = newstate.rTz_old;
00475     }
00476     else {
00477 
00478       TEST_FOR_EXCEPTION(newstate.P == Teuchos::null,std::invalid_argument,
00479                          "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
00480 
00481       TEST_FOR_EXCEPTION(newstate.Ap == Teuchos::null,std::invalid_argument,
00482                          "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
00483 
00484       TEST_FOR_EXCEPTION(newstate.r == Teuchos::null,std::invalid_argument,
00485                          "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
00486 
00487       TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00488                          "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
00489 
00490       TEST_FOR_EXCEPTION(newstate.U == Teuchos::null,std::invalid_argument,
00491                          "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
00492 
00493       TEST_FOR_EXCEPTION(newstate.AU == Teuchos::null,std::invalid_argument,
00494                          "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
00495 
00496       TEST_FOR_EXCEPTION(newstate.Alpha == Teuchos::null,std::invalid_argument,
00497                          "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
00498 
00499       TEST_FOR_EXCEPTION(newstate.Beta == Teuchos::null,std::invalid_argument,
00500                          "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
00501 
00502       TEST_FOR_EXCEPTION(newstate.D == Teuchos::null,std::invalid_argument,
00503                          "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
00504 
00505       TEST_FOR_EXCEPTION(newstate.Delta == Teuchos::null,std::invalid_argument,
00506                          "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
00507 
00508       TEST_FOR_EXCEPTION(newstate.LUUTAU == Teuchos::null,std::invalid_argument,
00509                          "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
00510 
00511       TEST_FOR_EXCEPTION(newstate.ipiv == Teuchos::null,std::invalid_argument,
00512                          "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
00513 
00514       TEST_FOR_EXCEPTION(newstate.rTz_old == Teuchos::null,std::invalid_argument,
00515                          "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
00516 
00517     }
00518 
00519     // the solver is initialized
00520     initialized_ = true;
00521 
00522   }
00523 
00525   // Iterate until the status test informs us we should stop.
00526   template <class ScalarType, class MV, class OP>
00527   void RCGIter<ScalarType,MV,OP>::iterate()
00528   {
00529     TEST_FOR_EXCEPTION( initialized_ == false, RCGIterFailure,
00530                         "Belos::RCGIter::iterate(): RCGIter class not initialized." );
00531     
00532     // We'll need LAPACK
00533     Teuchos::LAPACK<int,ScalarType> lapack;
00534 
00535     // Create convenience variables for zero and one.
00536     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00537     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00538 
00539     // Allocate memory for scalars
00540     std::vector<int> index(1);
00541     Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
00542 
00543     // Get the current solution std::vector.
00544     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00545  
00546     // Check that the current solution std::vector only has one column.
00547     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, RCGIterFailure,
00548                         "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
00549     
00550     // Compute the current search dimension. 
00551     int searchDim = numBlocks_+1;
00552 
00553     // index of iteration within current cycle
00554     int i_ = 0;
00555 
00557     // iterate until the status test tells us to stop.
00558     //
00559     // also break if our basis is full
00560     //
00561     Teuchos::RCP<MV> p_ = Teuchos::null;
00562     Teuchos::RCP<MV> pnext_ = Teuchos::null;
00563     while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {
00564 
00565       // Ap = A*p;
00566       index.resize( 1 );
00567       index[0] = i_;
00568       p_  = MVT::CloneView( *P_,  index );
00569       lp_->applyOp( *p_, *Ap_ );
00570 
00571       // d = p'*Ap;
00572       MVT::MvTransMv( one, *p_, *Ap_, pAp );
00573       (*D_)(i_,0) = pAp(0,0);
00574 
00575       // alpha = rTz_old / pAp
00576       (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
00577 
00578       // Check that alpha is a positive number
00579       TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero, RCGIterFailure, "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00580 
00581       // x = x + (alpha * p);
00582       MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
00583       lp_->updateSolution();
00584 
00585       // r = r - (alpha * Ap);
00586       MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
00587 
00588       std::vector<MagnitudeType> norm(1);
00589       MVT::MvNorm( *r_, norm );
00590 //printf("i = %i\tnorm(r) = %e\n",i_,norm[0]);
00591 
00592       // z = M\r
00593       if ( lp_->getLeftPrec() != Teuchos::null ) {
00594         lp_->applyLeftPrec( *r_, *z_ );
00595       } else {
00596         z_ = r_;
00597       }
00598 
00599       // rTz_new = r'*z;
00600       MVT::MvTransMv( one, *r_, *z_, rTz );
00601 
00602       // beta = rTz_new/rTz_old;
00603       (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
00604 
00605       // rTz_old = rTz_new;
00606       (*rTz_old_)(0,0) = rTz(0,0);
00607 
00608       // get pointer for next p
00609       index.resize( 1 );
00610       index[0] = i_+1;
00611       pnext_ = MVT::CloneView( *P_,  index );
00612 
00613       if (existU_) {
00614         // mu = UTAU \ (AU'*z);
00615         Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
00616         MVT::MvTransMv( one, *AU_, *z_, mu );
00617         char TRANS = 'N';
00618         int info;
00619         lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
00620         TEST_FOR_EXCEPTION(info != 0, RCGIterLAPACKFailure,
00621                            "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
00622         // p = -(U*mu) + (beta*p) + z (in two steps)
00623         // p = (beta*p) + z;
00624         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00625         // pnext = -(U*mu) + (one)*pnext;
00626         MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
00627       }
00628       else {
00629         // p = (beta*p) + z;
00630         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00631       }
00632 
00633       // Done with this view; release pointer
00634       p_ = Teuchos::null;
00635       pnext_ = Teuchos::null;
00636 
00637       // increment iteration count and dimension index  
00638       i_++;
00639       iter_++;
00640       curDim_++;
00641 
00642     } // end while (statusTest == false)
00643     
00644    }
00645 
00646 } // end Belos namespace
00647  
00648 #endif /* BELOS_RCG_ITER_HPP */

Generated on Tue Jul 13 09:27:03 2010 for Belos by  doxygen 1.4.7