Belos Version of the Day
BelosPseudoBlockStochasticCGIter.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_STOCHASTIC_CG_ITER_HPP
00043 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosStochasticCGIteration.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 
00082 namespace Belos {
00083   
00084   template<class ScalarType, class MV, class OP>
00085   class PseudoBlockStochasticCGIter : virtual public StochasticCGIteration<ScalarType,MV,OP> {
00086     
00087   public:
00088     
00089     //
00090     // Convenience typedefs
00091     //
00092     typedef MultiVecTraits<ScalarType,MV> MVT;
00093     typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00094     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00095     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00096     typedef typename SCT::magnitudeType MagnitudeType;
00097     
00099 
00100     
00106     PseudoBlockStochasticCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00107          const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00108          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00109          Teuchos::ParameterList &params );
00110     
00112     virtual ~PseudoBlockStochasticCGIter() {};
00114     
00115     
00117 
00118     
00132     void iterate();
00133     
00154     void initializeCG(StochasticCGIterationState<ScalarType,MV> newstate);
00155     
00159     void initialize()
00160     {
00161       StochasticCGIterationState<ScalarType,MV> empty;
00162       initializeCG(empty);
00163     }
00164     
00172     StochasticCGIterationState<ScalarType,MV> getState() const {
00173       StochasticCGIterationState<ScalarType,MV> state;
00174       state.R = R_;
00175       state.P = P_;
00176       state.AP = AP_;
00177       state.Z = Z_;
00178       state.Y = Y_;
00179       return state;
00180     }
00181     
00183     
00184     
00186 
00187     
00189     int getNumIters() const { return iter_; }
00190     
00192     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00193     
00196     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00197     
00199 
00201     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00202     
00204     Teuchos::RCP<MV> getStochasticVector() const { return Y_; }
00205 
00207     
00209 
00210     
00212     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00213     
00215     int getBlockSize() const { return 1; }
00216     
00218     void setBlockSize(int blockSize) { 
00219       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00220        "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
00221     }
00222     
00224     bool isInitialized() { return initialized_; }
00225     
00227     
00228   private:
00229 
00231     inline ScalarType normal() {
00232       // Do all of the calculations with doubles, because that is what the Odeh and Evans 1974 constants are for.
00233       // Then cast to ScalarType.
00234       
00235       const double p0 = -0.322232431088;
00236       const double p1 = -1.0;
00237       const double p2 = -0.342242088547; 
00238       const double p3 = -0.204231210245e-1;
00239       const double p4 = -0.453642210148e-4;
00240       const double q0 =  0.993484626060e-1;
00241       const double q1 =  0.588581570495;
00242       const double q2 =  0.531103462366;
00243       const double q3 =  0.103537752850;
00244       const double q4 =  0.38560700634e-2;
00245       double r,p,q,y,z;
00246       
00247       // Get a random number (-1,1) and rescale to (0,1). 
00248       r=0.5*(Teuchos::ScalarTraits<double>::random() + 1.0);
00249       
00250       // Odeh and Evans algorithm (as modified by Park & Geyer)
00251       if(r < 0.5) y=std::sqrt(-2.0 * log(r));
00252       else y=std::sqrt(-2.0 * log(1.0 - r));
00253       
00254       p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
00255       q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
00256       
00257       if(r < 0.5) z = (p / q) - y;
00258       else z = y - (p / q);
00259       
00260       return Teuchos::as<ScalarType,double>(z);
00261     }
00262        
00263     //
00264     // Classes inputed through constructor that define the linear problem to be solved.
00265     //
00266     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00267     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00268     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00269     
00270     //
00271     // Algorithmic parameters
00272     //  
00273     // numRHS_ is the current number of linear systems being solved.
00274     int numRHS_;
00275     
00276     // 
00277     // Current solver state
00278     //
00279     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00280     // is capable of running; _initialize is controlled  by the initialize() member method
00281     // For the implications of the state of initialized_, please see documentation for initialize()
00282     bool initialized_;
00283     
00284     // Current number of iterations performed.
00285     int iter_;
00286 
00287     // Current number of iterations performed.
00288     bool assertPositiveDefiniteness_;
00289 
00290     // 
00291     // State Storage
00292     //
00293     // Residual
00294     Teuchos::RCP<MV> R_;
00295     //
00296     // Preconditioned residual
00297     Teuchos::RCP<MV> Z_;
00298     //
00299     // Direction vector
00300     Teuchos::RCP<MV> P_;
00301     //
00302     // Operator applied to direction vector
00303     Teuchos::RCP<MV> AP_;
00304     //
00305     // Stochastic recurrence vector
00306     Teuchos::RCP<MV> Y_;
00307 
00308 
00309   };
00310   
00312   // Constructor.
00313   template<class ScalarType, class MV, class OP>
00314   PseudoBlockStochasticCGIter<ScalarType,MV,OP>::PseudoBlockStochasticCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00315                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00316                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00317                      Teuchos::ParameterList &params ):
00318     lp_(problem),
00319     om_(printer),
00320     stest_(tester),
00321     numRHS_(0),
00322     initialized_(false),
00323     iter_(0),
00324     assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
00325   {
00326   }
00327   
00328 
00330   // Initialize this iteration object
00331   template <class ScalarType, class MV, class OP>
00332   void PseudoBlockStochasticCGIter<ScalarType,MV,OP>::initializeCG(StochasticCGIterationState<ScalarType,MV> newstate)
00333   {
00334     // Check if there is any multivector to clone from.
00335     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
00336     Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
00337     TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
00338            "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
00339 
00340     // Get the multivector that is not null.
00341     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00342 
00343     // Get the number of right-hand sides we're solving for now.
00344     int numRHS = MVT::GetNumberVecs(*tmp);
00345     numRHS_ = numRHS;
00346   
00347     // Initialize the state storage
00348     // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
00349     if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
00350       R_  = MVT::Clone( *tmp, numRHS_ );
00351       Z_  = MVT::Clone( *tmp, numRHS_ );
00352       P_  = MVT::Clone( *tmp, numRHS_ );
00353       AP_ = MVT::Clone( *tmp, numRHS_ );
00354       Y_  = MVT::Clone( *tmp, numRHS_ );
00355     }
00356   
00357     // NOTE:  In StochasticCGIter R_, the initial residual, is required!!!  
00358     //
00359     std::string errstr("Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00360 
00361     // Create convenience variables for zero and one.
00362     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00363     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00364 
00365     if (!Teuchos::is_null(newstate.R)) {
00366 
00367       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.R) != MVText::GetGlobalLength(*R_),
00368                           std::invalid_argument, errstr );
00369       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
00370                           std::invalid_argument, errstr );
00371 
00372       // Copy basis vectors from newstate into V
00373       if (newstate.R != R_) {
00374         // copy over the initial residual (unpreconditioned).
00375   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00376       }
00377 
00378       // Compute initial direction vectors
00379       // Initially, they are set to the preconditioned residuals
00380       
00381       if ( lp_->getLeftPrec() != Teuchos::null ) {
00382         lp_->applyLeftPrec( *R_, *Z_ );
00383         if ( lp_->getRightPrec() != Teuchos::null ) {
00384           Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
00385           lp_->applyRightPrec( *Z_, *tmp2 );
00386           Z_ = tmp2;
00387         }
00388       }
00389       else if ( lp_->getRightPrec() != Teuchos::null ) {
00390   lp_->applyRightPrec( *R_, *Z_ ); 
00391       } 
00392       else {
00393   Z_ = R_;
00394       }
00395       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00396     }
00397     else {
00398 
00399       TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
00400                          "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
00401     }
00402 
00403     // The solver is initialized
00404     initialized_ = true;
00405   }
00406 
00407 
00409   // Iterate until the status test informs us we should stop.
00410   template <class ScalarType, class MV, class OP>
00411   void PseudoBlockStochasticCGIter<ScalarType,MV,OP>::iterate()
00412   {
00413     //
00414     // Allocate/initialize data structures
00415     //
00416     if (initialized_ == false) {
00417       initialize();
00418     }
00419 
00420     // Allocate memory for scalars.
00421     int i=0;
00422     std::vector<int> index(1);
00423     std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
00424     Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(numRHS_,numRHS_);
00425 
00426     // Create convenience variables for zero and one.
00427     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00428     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00429     
00430     // Get the current solution std::vector.
00431     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00432 
00433     // Compute first <r,z> a.k.a. rHz
00434     MVT::MvDot( *R_, *Z_, rHz );
00435 
00436     if ( assertPositiveDefiniteness_ )
00437         for (i=0; i<numRHS_; ++i)
00438             TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00439                                 CGIterateFailure,
00440                                 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
00441 
00443     // Iterate until the status test tells us to stop.
00444     //
00445     while (stest_->checkStatus(this) != Passed) {
00446       
00447       // Increment the iteration
00448       iter_++;
00449 
00450       // Multiply the current direction std::vector by A and store in AP_
00451       lp_->applyOp( *P_, *AP_ );
00452       
00453       // Compute alpha := <R_,Z_> / <P_,AP_>
00454       MVT::MvDot( *P_, *AP_, pAp );
00455 
00456       for (i=0; i<numRHS_; ++i) {
00457         if ( assertPositiveDefiniteness_ )
00458             // Check that pAp[i] is a positive number!
00459             TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
00460                                 CGIterateFailure,
00461                                 "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00462 
00463         alpha(i,i) = rHz[i] / pAp[i];
00464 
00465 
00466   // Compute the scaling parameter for the stochastic vector
00467   ScalarType z = normal();
00468   zeta(i,i) = z / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
00469       }
00470 
00471       //
00472       // Update the solution std::vector x := x + alpha * P_
00473       //
00474       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00475       lp_->updateSolution();
00476 
00477       // Updates the stochastic vector y := y + zeta * P_
00478       MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
00479 
00480       //
00481       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00482       //
00483       for (i=0; i<numRHS_; ++i) {
00484         rHz_old[i] = rHz[i];
00485       }
00486       //
00487       // Compute the new residual R_ := R_ - alpha * AP_
00488       //
00489       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00490       //
00491       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00492       // and the new direction std::vector p.
00493       //
00494       if ( lp_->getLeftPrec() != Teuchos::null ) {
00495         lp_->applyLeftPrec( *R_, *Z_ );
00496         if ( lp_->getRightPrec() != Teuchos::null ) {
00497           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
00498           lp_->applyRightPrec( *Z_, *tmp );
00499           Z_ = tmp;
00500         }
00501       }
00502       else if ( lp_->getRightPrec() != Teuchos::null ) {
00503         lp_->applyRightPrec( *R_, *Z_ );
00504       } 
00505       else {
00506         Z_ = R_;
00507       }
00508       //
00509       MVT::MvDot( *R_, *Z_, rHz );
00510       if ( assertPositiveDefiniteness_ )
00511           for (i=0; i<numRHS_; ++i)
00512               TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00513                                   CGIterateFailure,
00514                                   "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
00515       //
00516       // Update the search directions.
00517       for (i=0; i<numRHS_; ++i) {
00518         beta(i,i) = rHz[i] / rHz_old[i];
00519   index[0] = i;
00520   Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
00521   Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
00522         MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
00523       }
00524       //      
00525     } // end while (sTest_->checkStatus(this) != Passed)
00526   }
00527 
00528 } // end Belos namespace
00529 
00530 #endif /* BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines