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 
00067 namespace Belos {
00068   
00070 
00071   
00076   template <class ScalarType, class MV>
00077   struct RCGIterState {
00082     int curDim;
00083     
00085     Teuchos::RCP<MV> P;
00086 
00088     Teuchos::RCP<MV> Ap;
00089 
00091     Teuchos::RCP<MV> r;
00092  
00094     Teuchos::RCP<MV> z;
00095  
00097     bool existU;
00098    
00100     Teuchos::RCP<MV> U, AU;
00101 
00104     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha;
00105     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta;
00106     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D;
00107     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old;
00108 
00110     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta;
00111 
00113     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU;
00115     Teuchos::RCP<std::vector<int> > ipiv;
00116 
00117 
00118     RCGIterState() : curDim(0), P(Teuchos::null), Ap(Teuchos::null), r(Teuchos::null), 
00119                      z(Teuchos::null),
00120                      existU(false),
00121          U(Teuchos::null), AU(Teuchos::null),
00122          Alpha(Teuchos::null), Beta(Teuchos::null), D(Teuchos::null), rTz_old(Teuchos::null),
00123          Delta(Teuchos::null), LUUTAU(Teuchos::null), ipiv(Teuchos::null)
00124     {}
00125   };
00126   
00128   
00130 
00131   
00143   class RCGIterInitFailure : public BelosError {public:
00144     RCGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00145     {}};
00146 
00153   class RCGIterFailure : public BelosError {public:
00154     RCGIterFailure(const std::string& what_arg) : BelosError(what_arg)
00155     {}};
00156   
00163   class RCGIterOrthoFailure : public BelosError {public:
00164     RCGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00165     {}};
00166   
00173   class RCGIterLAPACKFailure : public BelosError {public:
00174     RCGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00175     {}};
00176   
00178   
00179   
00180   template<class ScalarType, class MV, class OP>
00181   class RCGIter : virtual public Iteration<ScalarType,MV,OP> {
00182     
00183   public:
00184     
00185     //
00186     // Convenience typedefs
00187     //
00188     typedef MultiVecTraits<ScalarType,MV> MVT;
00189     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00190     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00191     typedef typename SCT::magnitudeType MagnitudeType;
00192     
00194 
00195     
00203     RCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00204     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00205     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00206     Teuchos::ParameterList &params );
00207     
00209     virtual ~RCGIter() {};
00211     
00212     
00214 
00215     
00228     void iterate();
00229     
00244     void initialize(RCGIterState<ScalarType,MV> &newstate);
00245  
00249     void initialize()
00250     {
00251       RCGIterState<ScalarType,MV> empty;
00252       initialize(empty);
00253     }
00254 
00256     
00258 
00259     
00261     int getNumIters() const { return iter_; }
00262     
00264     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00265     
00268     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return r_; }
00269     
00271 
00276     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00277     
00279     int getCurSubspaceDim() const { 
00280       if (!initialized_) return 0;
00281       return curDim_;
00282     };
00283     
00285     int getMaxSubspaceDim() const { return numBlocks_+1; }
00286     
00288     
00289     
00291 
00292     
00294     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00295     
00297     int getNumBlocks() const { return numBlocks_; }
00298     
00300     void setNumBlocks(int numBlocks) { setSize( recycleBlocks_, numBlocks ); };
00301     
00303     int getRecycledBlocks() const { return recycleBlocks_; }
00304 
00306     void setRecycledBlocks(int recycleBlocks) { setSize( recycleBlocks, numBlocks_ ); };
00307     
00309     int getBlockSize() const { return 1; }
00310     
00312     void setBlockSize(int blockSize) {
00313       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00314        "Belos::RCGIter::setBlockSize(): Cannot use a block size that is not one.");
00315     }
00316 
00318     void setSize( int recycleBlocks, int numBlocks );
00319 
00321     bool isInitialized() { return initialized_; }
00322     
00324     
00325   private:
00326     
00327     //
00328     // Internal methods
00329     //
00330     
00331     //
00332     // Classes input through constructor that define the linear problem to be solved.
00333     //
00334     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00335     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00336     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00337     
00338     //
00339     // Algorithmic parameters
00340     //  
00341     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00342     int numBlocks_; 
00343 
00344     // recycleBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
00345     int recycleBlocks_; 
00346     
00347     // 
00348     // Current solver state
00349     //
00350     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00351     // is capable of running; _initialize is controlled  by the initialize() member method
00352     // For the implications of the state of initialized_, please see documentation for initialize()
00353     bool initialized_;
00354     
00355     // Current subspace dimension, and number of iterations performed.
00356     int curDim_, iter_;
00357     
00358     // 
00359     // State Storage
00360     //
00361     // Search vectors
00362     Teuchos::RCP<MV> P_;
00363     //
00364     // A times current search vector
00365     Teuchos::RCP<MV> Ap_;
00366     //
00367     // Residual vector
00368     Teuchos::RCP<MV> r_;
00369     //
00370     // Preconditioned residual
00371     Teuchos::RCP<MV> z_;
00372     //
00373     // Flag to indicate that the recycle space should be used
00374     bool existU_;
00375     // Recycled subspace and its image
00376     Teuchos::RCP<MV> U_, AU_;
00377     //
00378     // Coefficients arising in RCG iteration
00379     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_,Beta_,D_;
00380     //
00381     // Solutions to local least-squares problems
00382     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
00383     //
00384     // The LU factorization of the matrix U^T A U
00385     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
00386     //
00387     // Data from LU factorization of UTAU
00388     Teuchos::RCP<std::vector<int> > ipiv_;
00389     //
00390     // The scalar r'*z
00391     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
00392   };
00393   
00395   // Constructor.
00396   template<class ScalarType, class MV, class OP>
00397   RCGIter<ScalarType,MV,OP>::RCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00398              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00399              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00400              Teuchos::ParameterList &params ):
00401     lp_(problem),
00402     om_(printer),
00403     stest_(tester),
00404     numBlocks_(0),
00405     recycleBlocks_(0),
00406     initialized_(false),
00407     curDim_(0),
00408     iter_(0)
00409   {
00410     // Get the maximum number of blocks allowed for this Krylov subspace
00411     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00412                        "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00413     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00414 
00415     TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,
00416                        "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00417     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00418 
00419     // Set the number of blocks and allocate data
00420     setSize( rb, nb );
00421   }
00422   
00424   // Set the block size and make necessary adjustments.
00425   template <class ScalarType, class MV, class OP>
00426   void RCGIter<ScalarType,MV,OP>::setSize( int recycleBlocks, int numBlocks )
00427   {
00428 
00429     TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
00430     TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
00431     TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument, "Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
00432 
00433     numBlocks_ = numBlocks;
00434     recycleBlocks_ = recycleBlocks;
00435 
00436   }
00437 
00439   // Initialize this iteration object
00440   template <class ScalarType, class MV, class OP>
00441   void RCGIter<ScalarType,MV,OP>::initialize(RCGIterState<ScalarType,MV> &newstate)
00442   {
00443 
00444     if (newstate.P != Teuchos::null && 
00445         newstate.Ap != Teuchos::null && 
00446         newstate.r != Teuchos::null && 
00447         newstate.z != Teuchos::null && 
00448         newstate.U != Teuchos::null && 
00449         newstate.AU != Teuchos::null && 
00450         newstate.Alpha != Teuchos::null && 
00451         newstate.Beta != Teuchos::null && 
00452         newstate.D != Teuchos::null && 
00453         newstate.Delta != Teuchos::null && 
00454         newstate.LUUTAU != Teuchos::null &&
00455         newstate.ipiv != Teuchos::null &&
00456         newstate.rTz_old != Teuchos::null) {
00457 
00458       curDim_ = newstate.curDim;
00459       P_ = newstate.P;
00460       Ap_ = newstate.Ap;
00461       r_ = newstate.r;
00462       z_ = newstate.z;
00463       existU_ = newstate.existU;
00464       U_ = newstate.U;
00465       AU_ = newstate.AU;
00466       Alpha_ = newstate.Alpha;
00467       Beta_ = newstate.Beta;
00468       D_ = newstate.D;
00469       Delta_ = newstate.Delta;
00470       LUUTAU_ = newstate.LUUTAU;
00471       ipiv_ = newstate.ipiv;
00472       rTz_old_ = newstate.rTz_old;
00473     }
00474     else {
00475 
00476       TEST_FOR_EXCEPTION(newstate.P == Teuchos::null,std::invalid_argument,
00477                          "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
00478 
00479       TEST_FOR_EXCEPTION(newstate.Ap == Teuchos::null,std::invalid_argument,
00480                          "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
00481 
00482       TEST_FOR_EXCEPTION(newstate.r == Teuchos::null,std::invalid_argument,
00483                          "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
00484 
00485       TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00486                          "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
00487 
00488       TEST_FOR_EXCEPTION(newstate.U == Teuchos::null,std::invalid_argument,
00489                          "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
00490 
00491       TEST_FOR_EXCEPTION(newstate.AU == Teuchos::null,std::invalid_argument,
00492                          "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
00493 
00494       TEST_FOR_EXCEPTION(newstate.Alpha == Teuchos::null,std::invalid_argument,
00495                          "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
00496 
00497       TEST_FOR_EXCEPTION(newstate.Beta == Teuchos::null,std::invalid_argument,
00498                          "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
00499 
00500       TEST_FOR_EXCEPTION(newstate.D == Teuchos::null,std::invalid_argument,
00501                          "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
00502 
00503       TEST_FOR_EXCEPTION(newstate.Delta == Teuchos::null,std::invalid_argument,
00504                          "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
00505 
00506       TEST_FOR_EXCEPTION(newstate.LUUTAU == Teuchos::null,std::invalid_argument,
00507                          "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
00508 
00509       TEST_FOR_EXCEPTION(newstate.ipiv == Teuchos::null,std::invalid_argument,
00510                          "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
00511 
00512       TEST_FOR_EXCEPTION(newstate.rTz_old == Teuchos::null,std::invalid_argument,
00513                          "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
00514 
00515     }
00516 
00517     // the solver is initialized
00518     initialized_ = true;
00519 
00520   }
00521 
00523   // Iterate until the status test informs us we should stop.
00524   template <class ScalarType, class MV, class OP>
00525   void RCGIter<ScalarType,MV,OP>::iterate()
00526   {
00527     TEST_FOR_EXCEPTION( initialized_ == false, RCGIterFailure,
00528                         "Belos::RCGIter::iterate(): RCGIter class not initialized." );
00529     
00530     // We'll need LAPACK
00531     Teuchos::LAPACK<int,ScalarType> lapack;
00532 
00533     // Create convenience variables for zero and one.
00534     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00535     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00536 
00537     // Allocate memory for scalars
00538     std::vector<int> index(1);
00539     Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
00540 
00541     // Get the current solution std::vector.
00542     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00543  
00544     // Check that the current solution std::vector only has one column.
00545     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, RCGIterFailure,
00546                         "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
00547     
00548     // Compute the current search dimension. 
00549     int searchDim = numBlocks_+1;
00550 
00551     // index of iteration within current cycle
00552     int i_ = 0;
00553 
00555     // iterate until the status test tells us to stop.
00556     //
00557     // also break if our basis is full
00558     //
00559     Teuchos::RCP<MV> p_ = Teuchos::null;
00560     Teuchos::RCP<MV> pnext_ = Teuchos::null;
00561     while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {
00562 
00563       // Ap = A*p;
00564       index.resize( 1 );
00565       index[0] = i_;
00566       p_  = MVT::CloneView( *P_,  index );
00567       lp_->applyOp( *p_, *Ap_ );
00568 
00569       // d = p'*Ap;
00570       MVT::MvTransMv( one, *p_, *Ap_, pAp );
00571       (*D_)(i_,0) = pAp(0,0);
00572 
00573       // alpha = rTz_old / pAp
00574       (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
00575 
00576       // Check that alpha is a positive number
00577       TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero, RCGIterFailure, "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00578 
00579       // x = x + (alpha * p);
00580       MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
00581       lp_->updateSolution();
00582 
00583       // r = r - (alpha * Ap);
00584       MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
00585 
00586       std::vector<MagnitudeType> norm(1);
00587       MVT::MvNorm( *r_, norm );
00588 //printf("i = %i\tnorm(r) = %e\n",i_,norm[0]);
00589 
00590       // z = M\r
00591       if ( lp_->getLeftPrec() != Teuchos::null ) {
00592         lp_->applyLeftPrec( *r_, *z_ );
00593       } else {
00594         z_ = r_;
00595       }
00596 
00597       // rTz_new = r'*z;
00598       MVT::MvTransMv( one, *r_, *z_, rTz );
00599 
00600       // beta = rTz_new/rTz_old;
00601       (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
00602 
00603       // rTz_old = rTz_new;
00604       (*rTz_old_)(0,0) = rTz(0,0);
00605 
00606       // get pointer for next p
00607       index.resize( 1 );
00608       index[0] = i_+1;
00609       pnext_ = MVT::CloneView( *P_,  index );
00610 
00611       if (existU_) {
00612         // mu = UTAU \ (AU'*z);
00613         Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
00614         MVT::MvTransMv( one, *AU_, *z_, mu );
00615         char TRANS = 'N';
00616         int info;
00617         lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
00618         TEST_FOR_EXCEPTION(info != 0, RCGIterLAPACKFailure,
00619                            "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
00620         // p = -(U*mu) + (beta*p) + z (in two steps)
00621         // p = (beta*p) + z;
00622         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00623         // pnext = -(U*mu) + (one)*pnext;
00624         MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
00625       }
00626       else {
00627         // p = (beta*p) + z;
00628         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00629       }
00630 
00631       // Done with this view; release pointer
00632       p_ = Teuchos::null;
00633       pnext_ = Teuchos::null;
00634 
00635       // increment iteration count and dimension index  
00636       i_++;
00637       iter_++;
00638       curDim_++;
00639 
00640     } // end while (statusTest == false)
00641     
00642    }
00643 
00644 } // end Belos namespace
00645  
00646 #endif /* BELOS_RCG_ITER_HPP */

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7