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