Belos Package Browser (Single Doxygen Collection) Development
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     existU_(false)
00425   {
00426     // Get the maximum number of blocks allowed for this Krylov subspace
00427     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
00428                        "Belos::RCGIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00429     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00430 
00431     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,
00432                        "Belos::RCGIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00433     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00434 
00435     // Set the number of blocks and allocate data
00436     setSize( rb, nb );
00437   }
00438 
00440   // Set the block size and make necessary adjustments.
00441   template <class ScalarType, class MV, class OP>
00442   void RCGIter<ScalarType,MV,OP>::setSize( int recycleBlocks, int numBlocks )
00443   {
00444 
00445     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Num Blocks\".");
00446     TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks <= 0, std::invalid_argument, "Belos::RCGIter::setSize() was passed a non-positive argument for \"Recycled Blocks\".");
00447     TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks >= numBlocks, std::invalid_argument, "Belos::RCGIter::setSize() the number of recycled blocks is larger than the allowable subspace.");
00448 
00449     numBlocks_ = numBlocks;
00450     recycleBlocks_ = recycleBlocks;
00451 
00452   }
00453 
00455   // Initialize this iteration object
00456   template <class ScalarType, class MV, class OP>
00457   void RCGIter<ScalarType,MV,OP>::initialize(RCGIterState<ScalarType,MV> &newstate)
00458   {
00459 
00460     if (newstate.P != Teuchos::null &&
00461         newstate.Ap != Teuchos::null &&
00462         newstate.r != Teuchos::null &&
00463         newstate.z != Teuchos::null &&
00464         newstate.U != Teuchos::null &&
00465         newstate.AU != Teuchos::null &&
00466         newstate.Alpha != Teuchos::null &&
00467         newstate.Beta != Teuchos::null &&
00468         newstate.D != Teuchos::null &&
00469         newstate.Delta != Teuchos::null &&
00470         newstate.LUUTAU != Teuchos::null &&
00471         newstate.ipiv != Teuchos::null &&
00472         newstate.rTz_old != Teuchos::null) {
00473 
00474       curDim_ = newstate.curDim;
00475       P_ = newstate.P;
00476       Ap_ = newstate.Ap;
00477       r_ = newstate.r;
00478       z_ = newstate.z;
00479       existU_ = newstate.existU;
00480       U_ = newstate.U;
00481       AU_ = newstate.AU;
00482       Alpha_ = newstate.Alpha;
00483       Beta_ = newstate.Beta;
00484       D_ = newstate.D;
00485       Delta_ = newstate.Delta;
00486       LUUTAU_ = newstate.LUUTAU;
00487       ipiv_ = newstate.ipiv;
00488       rTz_old_ = newstate.rTz_old;
00489     }
00490     else {
00491 
00492       TEUCHOS_TEST_FOR_EXCEPTION(newstate.P == Teuchos::null,std::invalid_argument,
00493                          "Belos::RCGIter::initialize(): RCGIterState does not have P initialized.");
00494 
00495       TEUCHOS_TEST_FOR_EXCEPTION(newstate.Ap == Teuchos::null,std::invalid_argument,
00496                          "Belos::RCGIter::initialize(): RCGIterState does not have Ap initialized.");
00497 
00498       TEUCHOS_TEST_FOR_EXCEPTION(newstate.r == Teuchos::null,std::invalid_argument,
00499                          "Belos::RCGIter::initialize(): RCGIterState does not have r initialized.");
00500 
00501       TEUCHOS_TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
00502                          "Belos::RCGIter::initialize(): RCGIterState does not have z initialized.");
00503 
00504       TEUCHOS_TEST_FOR_EXCEPTION(newstate.U == Teuchos::null,std::invalid_argument,
00505                          "Belos::RCGIter::initialize(): RCGIterState does not have U initialized.");
00506 
00507       TEUCHOS_TEST_FOR_EXCEPTION(newstate.AU == Teuchos::null,std::invalid_argument,
00508                          "Belos::RCGIter::initialize(): RCGIterState does not have AU initialized.");
00509 
00510       TEUCHOS_TEST_FOR_EXCEPTION(newstate.Alpha == Teuchos::null,std::invalid_argument,
00511                          "Belos::RCGIter::initialize(): RCGIterState does not have Alpha initialized.");
00512 
00513       TEUCHOS_TEST_FOR_EXCEPTION(newstate.Beta == Teuchos::null,std::invalid_argument,
00514                          "Belos::RCGIter::initialize(): RCGIterState does not have Beta initialized.");
00515 
00516       TEUCHOS_TEST_FOR_EXCEPTION(newstate.D == Teuchos::null,std::invalid_argument,
00517                          "Belos::RCGIter::initialize(): RCGIterState does not have D initialized.");
00518 
00519       TEUCHOS_TEST_FOR_EXCEPTION(newstate.Delta == Teuchos::null,std::invalid_argument,
00520                          "Belos::RCGIter::initialize(): RCGIterState does not have Delta initialized.");
00521 
00522       TEUCHOS_TEST_FOR_EXCEPTION(newstate.LUUTAU == Teuchos::null,std::invalid_argument,
00523                          "Belos::RCGIter::initialize(): RCGIterState does not have LUUTAU initialized.");
00524 
00525       TEUCHOS_TEST_FOR_EXCEPTION(newstate.ipiv == Teuchos::null,std::invalid_argument,
00526                          "Belos::RCGIter::initialize(): RCGIterState does not have ipiv initialized.");
00527 
00528       TEUCHOS_TEST_FOR_EXCEPTION(newstate.rTz_old == Teuchos::null,std::invalid_argument,
00529                          "Belos::RCGIter::initialize(): RCGIterState does not have rTz_old initialized.");
00530 
00531     }
00532 
00533     // the solver is initialized
00534     initialized_ = true;
00535 
00536   }
00537 
00539   // Iterate until the status test informs us we should stop.
00540   template <class ScalarType, class MV, class OP>
00541   void RCGIter<ScalarType,MV,OP>::iterate()
00542   {
00543     TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, RCGIterFailure,
00544                         "Belos::RCGIter::iterate(): RCGIter class not initialized." );
00545 
00546     // We'll need LAPACK
00547     Teuchos::LAPACK<int,ScalarType> lapack;
00548 
00549     // Create convenience variables for zero and one.
00550     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00551     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00552 
00553     // Allocate memory for scalars
00554     std::vector<int> index(1);
00555     Teuchos::SerialDenseMatrix<int,ScalarType> pAp(1,1), rTz(1,1);
00556 
00557     // Get the current solution std::vector.
00558     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00559 
00560     // Check that the current solution std::vector only has one column.
00561     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, RCGIterFailure,
00562                         "Belos::RCGIter::iterate(): current linear system has more than one std::vector!" );
00563 
00564     // Compute the current search dimension.
00565     int searchDim = numBlocks_+1;
00566 
00567     // index of iteration within current cycle
00568     int i_ = 0;
00569 
00571     // iterate until the status test tells us to stop.
00572     //
00573     // also break if our basis is full
00574     //
00575     Teuchos::RCP<const MV> p_ = Teuchos::null;
00576     Teuchos::RCP<MV> pnext_ = Teuchos::null;
00577     while (stest_->checkStatus(this) != Passed && curDim_+1 <= searchDim) {
00578 
00579       // Ap = A*p;
00580       index.resize( 1 );
00581       index[0] = i_;
00582       p_  = MVT::CloneView( *P_,  index );
00583       lp_->applyOp( *p_, *Ap_ );
00584 
00585       // d = p'*Ap;
00586       MVT::MvTransMv( one, *p_, *Ap_, pAp );
00587       (*D_)(i_,0) = pAp(0,0);
00588 
00589       // alpha = rTz_old / pAp
00590       (*Alpha_)(i_,0) = (*rTz_old_)(0,0) / pAp(0,0);
00591 
00592       // Check that alpha is a positive number
00593       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp(0,0)) <= zero, RCGIterFailure, "Belos::RCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00594 
00595       // x = x + (alpha * p);
00596       MVT::MvAddMv( one, *cur_soln_vec, (*Alpha_)(i_,0), *p_, *cur_soln_vec );
00597       lp_->updateSolution();
00598 
00599       // r = r - (alpha * Ap);
00600       MVT::MvAddMv( one, *r_, -(*Alpha_)(i_,0), *Ap_, *r_ );
00601 
00602       std::vector<MagnitudeType> norm(1);
00603       MVT::MvNorm( *r_, norm );
00604 //printf("i = %i\tnorm(r) = %e\n",i_,norm[0]);
00605 
00606       // z = M\r
00607       if ( lp_->getLeftPrec() != Teuchos::null ) {
00608         lp_->applyLeftPrec( *r_, *z_ );
00609       }
00610       else if ( lp_->getRightPrec() != Teuchos::null ) {
00611         lp_->applyRightPrec( *r_, *z_ );
00612       }
00613       else {
00614         z_ = r_;
00615       }
00616 
00617       // rTz_new = r'*z;
00618       MVT::MvTransMv( one, *r_, *z_, rTz );
00619 
00620       // beta = rTz_new/rTz_old;
00621       (*Beta_)(i_,0) = rTz(0,0) / (*rTz_old_)(0,0);
00622 
00623       // rTz_old = rTz_new;
00624       (*rTz_old_)(0,0) = rTz(0,0);
00625 
00626       // get pointer for next p
00627       index.resize( 1 );
00628       index[0] = i_+1;
00629       pnext_ = MVT::CloneViewNonConst( *P_,  index );
00630 
00631       if (existU_) {
00632         // mu = UTAU \ (AU'*z);
00633         Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, i_ );
00634         MVT::MvTransMv( one, *AU_, *z_, mu );
00635         char TRANS = 'N';
00636         int info;
00637         lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAU_->values(), LUUTAU_->stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
00638         TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGIterLAPACKFailure,
00639                            "Belos::RCGIter::solve(): LAPACK GETRS failed to compute a solution.");
00640         // p = -(U*mu) + (beta*p) + z (in two steps)
00641         // p = (beta*p) + z;
00642         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00643         // pnext = -(U*mu) + (one)*pnext;
00644         MVT::MvTimesMatAddMv( -one, *U_, mu, one, *pnext_ );
00645       }
00646       else {
00647         // p = (beta*p) + z;
00648         MVT::MvAddMv( (*Beta_)(i_,0), *p_, one, *z_, *pnext_ );
00649       }
00650 
00651       // Done with this view; release pointer
00652       p_ = Teuchos::null;
00653       pnext_ = Teuchos::null;
00654 
00655       // increment iteration count and dimension index
00656       i_++;
00657       iter_++;
00658       curDim_++;
00659 
00660     } // end while (statusTest == false)
00661 
00662    }
00663 
00664 } // end Belos namespace
00665 
00666 #endif /* BELOS_RCG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines