Belos Package Browser (Single Doxygen Collection) Development
BelosMinresIter.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_MINRES_ITER_HPP
00043 #define BELOS_MINRES_ITER_HPP
00044 
00061 
00062 #include "BelosConfigDefs.hpp"
00063 #include "BelosTypes.hpp"
00064 #include "BelosMinresIteration.hpp"
00065 
00066 #include "BelosLinearProblem.hpp"
00067 #include "BelosOutputManager.hpp"
00068 #include "BelosStatusTest.hpp"
00069 #include "BelosOperatorTraits.hpp"
00070 #include "BelosMultiVecTraits.hpp"
00071 
00072 #include "Teuchos_SerialDenseMatrix.hpp"
00073 #include "Teuchos_SerialDenseVector.hpp"
00074 #include "Teuchos_ScalarTraits.hpp"
00075 #include "Teuchos_ParameterList.hpp"
00076 #include "Teuchos_TimeMonitor.hpp"
00077 #include "Teuchos_BLAS.hpp"
00078 
00079 namespace Belos {
00080 
00094 template<class ScalarType, class MV, class OP>
00095 class MinresIter : virtual public MinresIteration<ScalarType,MV,OP> {
00096 
00097   public:
00098 
00099   //
00100   // Convenience typedefs
00101   //
00102   typedef MultiVecTraits< ScalarType, MV > MVT;
00103   typedef OperatorTraits< ScalarType, MV, OP > OPT;
00104   typedef Teuchos::ScalarTraits< ScalarType > SCT;
00105   typedef typename SCT::magnitudeType MagnitudeType;
00106   typedef Teuchos::ScalarTraits< MagnitudeType > SMT;
00107 
00109 
00110 
00119   MinresIter (const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > >& problem,
00120         const Teuchos::RCP< OutputManager< ScalarType > > &        printer,
00121         const Teuchos::RCP< StatusTest< ScalarType, MV, OP > >&    tester,
00122         const Teuchos::ParameterList& params);
00123 
00125   virtual ~MinresIter() {};
00127 
00128 
00130 
00131 
00146   void iterate();
00147 
00162   void initializeMinres (MinresIterationState<ScalarType,MV> newstate);
00163 
00169   void initialize()
00170   {
00171     MinresIterationState<ScalarType,MV> empty;
00172     initializeMinres(empty);
00173   }
00174 
00181   MinresIterationState<ScalarType,MV> getState() const {
00182     if (! isInitialized())
00183       throw std::logic_error("getState() cannot be called unless "
00184            "the state has been initialized");
00185     MinresIterationState<ScalarType,MV> state;
00186     state.Y = Y_;
00187     state.R1 = R1_;
00188     state.R2 = R2_;
00189     state.W = W_;
00190     state.W1 = W1_;
00191     state.W2 = W2_;
00192     return state;
00193   }
00194 
00196 
00197   
00199 
00200 
00202   int getNumIters() const { return iter_; }
00203   
00205   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00206 
00209   Teuchos::RCP<const MV> 
00210   getNativeResiduals( std::vector<MagnitudeType> *norms ) const 
00211   { 
00212     if (norms != NULL)
00213       {
00214   std::vector<MagnitudeType>& theNorms = *norms;
00215   if (theNorms.size() < 1)
00216     theNorms.resize(1);
00217   theNorms[0] = phibar_;
00218       }
00219     return Teuchos::null; 
00220   }
00221 
00223 
00225   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00226 
00228   void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
00229 
00231 
00233 
00234 
00236   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00237 
00239   int getBlockSize() const { return 1; }
00240 
00242   void setBlockSize(int blockSize) {
00243     TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00244                        "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
00245   }
00246 
00248   bool isInitialized() const { return initialized_; }
00249   bool isInitialized() { return initialized_; }
00250 
00252 
00253   private:
00254 
00255   //
00256   // Internal methods
00257   //
00259   void setStateSize();
00260 
00261   //
00262   // Classes inputed through constructor that define the linear problem to be solved.
00263   //
00264   const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_;
00265   const Teuchos::RCP< OutputManager< ScalarType > >         om_;
00266   const Teuchos::RCP< StatusTest< ScalarType, MV, OP > >    stest_;
00267 
00268 
00276   bool initialized_;
00277 
00284   bool stateStorageInitialized_;
00285 
00287   int iter_;
00288 
00293   MagnitudeType phibar_;
00294 
00295   // 
00296   // State Storage
00297   //
00298 
00300   Teuchos::RCP< MV > Y_;
00302   Teuchos::RCP< MV > R1_;
00304   Teuchos::RCP< MV > R2_;
00306   Teuchos::RCP< MV > W_;
00308   Teuchos::RCP< MV > W1_;
00310   Teuchos::RCP< MV > W2_;
00311 
00319   Teuchos::SerialDenseMatrix<int,ScalarType> beta1_;
00320 
00321 };
00322 
00324   // Constructor.
00325   template<class ScalarType, class MV, class OP>
00326   MinresIter<ScalarType,MV,OP>::MinresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00327                                                    const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00328                                                    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00329                                                    const Teuchos::ParameterList &params ):
00330     lp_(problem),
00331     om_(printer),
00332     stest_(tester),
00333     initialized_(false),
00334     stateStorageInitialized_(false),
00335     iter_(0)
00336   {
00337   }
00338 
00340   // Setup the state storage.
00341   template <class ScalarType, class MV, class OP>
00342   void MinresIter<ScalarType,MV,OP>::setStateSize ()
00343   {
00344     if (!stateStorageInitialized_) {
00345 
00346       // Check if there is any multivector to clone from.
00347       Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
00348       Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
00349       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00350         stateStorageInitialized_ = false;
00351         return;
00352       }
00353       else {
00354 
00355         // Initialize the state storage
00356         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00357         if (Y_ == Teuchos::null) {
00358           // Get the multivector that is not null.
00359           Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00360           TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
00361                               std::invalid_argument,
00362                               "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00363           Y_  = MVT::Clone( *tmp, 1 );
00364           R1_ = MVT::Clone( *tmp, 1 );
00365           R2_ = MVT::Clone( *tmp, 1 );
00366           W_  = MVT::Clone( *tmp, 1 );
00367           W1_ = MVT::Clone( *tmp, 1 );
00368           W2_ = MVT::Clone( *tmp, 1 );
00369         }
00370         // State storage has now been initialized.
00371         stateStorageInitialized_ = true;
00372       }
00373     }
00374   }
00375 
00376 
00378   // Initialize this iteration object
00379   template <class ScalarType, class MV, class OP>
00380   void MinresIter<ScalarType,MV,OP>::initializeMinres(MinresIterationState<ScalarType,MV> newstate)
00381   {
00382     // Initialize the state storage if it isn't already.
00383     if (!stateStorageInitialized_) 
00384       setStateSize();
00385 
00386     TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
00387                         std::invalid_argument,
00388                         "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
00389 
00390     TEUCHOS_TEST_FOR_EXCEPTION( newstate.Y == Teuchos::null,
00391                         std::invalid_argument,
00392                         "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
00393 
00394     std::string errstr("Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
00395     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.Y) != MVT::GetVecLength(*Y_),
00396                         std::invalid_argument,
00397                         errstr );
00398     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.Y) != 1,
00399                         std::invalid_argument,
00400                         errstr );
00401 
00402     // Create convenience variables for zero, one.
00403     const ScalarType one = SCT::one();
00404     const MagnitudeType zero = SMT::zero();
00405 
00406     // Set up y and v for the first Lanczos vector v_1.
00407     // y  =  beta1_ P' v1,  where  P = C**(-1).
00408     // v is really P' v1.
00409     MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R2_ );
00410     MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R1_ );
00411 
00412     // Initialize the W's to 0.
00413     MVT::MvInit ( *W_ );
00414     MVT::MvInit ( *W2_ );
00415 
00416     if ( lp_->getLeftPrec() != Teuchos::null ) {
00417       lp_->applyLeftPrec( *newstate.Y, *Y_ );
00418     } 
00419     else {
00420       if (newstate.Y != Y_) {
00421         // copy over the initial residual (unpreconditioned).
00422         MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *Y_ );
00423       }
00424     }
00425 
00426     // beta1_ = b'*y;
00427     beta1_ = Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 );
00428     MVT::MvTransMv( one, *newstate.Y, *Y_, beta1_ );
00429 
00430     TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < zero,
00431                         std::invalid_argument,
00432                         "The preconditioner is not positive definite." );
00433 
00434     if( SCT::magnitude(beta1_(0,0)) == zero )
00435     {
00436         // X = 0
00437         Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00438         MVT::MvInit( *cur_soln_vec );
00439     }
00440 
00441     beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
00442 
00443     // The solver is initialized
00444     initialized_ = true;
00445   }
00446 
00447 
00449   // Iterate until the status test informs us we should stop.
00450   template <class ScalarType, class MV, class OP>
00451   void MinresIter<ScalarType,MV,OP>::iterate()
00452   {
00453     //
00454     // Allocate/initialize data structures
00455     //
00456     if (initialized_ == false) {
00457       initialize();
00458     }
00459 
00460     Teuchos::BLAS<int,ScalarType> blas;
00461 
00462     // Create convenience variables for zero and one.
00463     const ScalarType one = SCT::one();
00464     const MagnitudeType zero = SMT::zero();
00465 
00466     // Allocate memory for scalars.
00467     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00468     Teuchos::SerialDenseMatrix<int,ScalarType> beta( beta1_ );
00469     phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
00470     ScalarType shift = zero; // TODO Allow for proper shift.
00471 
00472     // Initialize a few variables.
00473     ScalarType oldBeta = zero;
00474     ScalarType epsln = zero;
00475     ScalarType cs = -one;
00476     ScalarType sn = zero;
00477     ScalarType dbar = zero;
00478 
00479     // Declare a few others that will be initialized in the loop.
00480     ScalarType oldeps;
00481     ScalarType delta;
00482     ScalarType gbar;
00483     ScalarType phi;
00484     ScalarType gamma;
00485 
00486     // Allocate workspace.
00487     Teuchos::RCP<MV> V    = MVT::Clone( *Y_, 1 );
00488     Teuchos::RCP<MV> tmpY, tmpW;  // Not allocated, just used to transfer ownership.
00489 
00490     // Get the current solution vector.
00491     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00492 
00493     // Check that the current solution vector only has one column. 
00494     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
00495                         MinresIterateFailure,
00496                         "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
00497 
00499     // Iterate until the status test tells us to stop.
00500     //
00501     while (stest_->checkStatus(this) != Passed) {
00502 
00503       // Increment the iteration
00504       iter_++;
00505 
00506       // Normalize previous vector.
00507       //   v = y / beta(0,0);
00508       MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
00509 
00510       // Apply operator.
00511       lp_->applyOp (*V, *Y_);
00512 
00513       // Apply shift
00514       if (shift != zero)
00515   MVT::MvAddMv (one, *Y_, -shift, *V, *Y_);
00516 
00517       if (iter_ > 1)
00518   MVT::MvAddMv (one, *Y_, -beta(0,0)/oldBeta, *R1_, *Y_);
00519 
00520       // alpha := dot(V, Y_)
00521       MVT::MvTransMv (one, *V, *Y_, alpha);
00522 
00523       // y := y - alpha/beta r2
00524       MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
00525 
00526       // r1 = r2;
00527       // r2 = y;
00528       tmpY = R1_;
00529       R1_ = R2_;
00530       R2_ = Y_;
00531       Y_ = tmpY;
00532 
00533       // apply left preconditioner
00534       if ( lp_->getLeftPrec() != Teuchos::null ) {
00535         lp_->applyLeftPrec( *R2_, *Y_ );
00536       } // else "y = r2"
00537       else {
00538         MVT::MvAddMv( one, *R2_, zero, *R2_, *Y_ );
00539       }
00540 
00541       // Get new beta.
00542       oldBeta = beta(0,0);
00543       MVT::MvTransMv( one, *R2_, *Y_, beta );
00544 
00545       // Intercept beta <= 0.
00546       //
00547       // Note: we don't try to test for nonzero imaginary component of
00548       // beta, because (a) it could be small and nonzero due to
00549       // rounding error in computing the inner product, and (b) it's
00550       // hard to tell how big "not small" should be, without computing
00551       // some error bounds (for example, by modifying the linear
00552       // algebra library to compute a posteriori rounding error bounds
00553       // for the inner product, and then changing
00554       // Belos::MultiVecTraits to make this information available).
00555       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) <= zero,
00556                           MinresIterateFailure,
00557                           "Belos::MinresIter::iterate(): Encountered nonpositi"
00558         "ve value " << beta(0,0) << " for r2^H*M*r2 at itera"
00559         "tion " << iter_ << ": MINRES cannot continue." );
00560       beta(0,0) = SCT::squareroot( beta(0,0) );
00561 
00562       // Apply previous rotation Q_{k-1} to get
00563       //
00564       //    [delta_k epsln_{k+1}] = [cs  sn][dbar_k  0         ]
00565       //    [gbar_k  dbar_{k+1} ]   [-sn cs][alpha_k beta_{k+1}].
00566       //
00567       oldeps = epsln;
00568       delta  = cs*dbar + sn*alpha(0,0);
00569       gbar   = sn*dbar - cs*alpha(0,0);
00570       epsln  =           sn*beta(0,0);
00571       dbar   =         - cs*beta(0,0);
00572 
00573       // Compute the next plane rotation Q_k.
00574       this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
00575 
00576       phi    = cs * phibar_; // phi_k
00577       phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ ); // phibar_{k+1}
00578 
00579       //  w1 = w2;
00580       //  w2 = w;
00581       MVT::MvAddMv( one, *W_, zero, *W_, *W1_ );
00582       tmpW = W1_;
00583       W1_ = W2_;
00584       W2_ = W_;
00585       W_ = tmpW;
00586 
00587       //  w = (v - oldeps*w1 - delta*w2) / gamma;
00588       MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
00589       MVT::MvAddMv( one, *W_, -delta,  *W2_, *W_ );
00590       MVT::MvScale( *W_, one / gamma );
00591 
00592       // Update x:
00593       // x = x + phi*w;
00594       MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
00595       lp_->updateSolution();
00596     } // end while (sTest_->checkStatus(this) != Passed)
00597   }
00598 
00599 
00601   // Compute the next plane rotation Qk.
00602   //   r = norm([a b]);
00603   //   c = a / r;
00604   //   s = b / r;
00605   template <class ScalarType, class MV, class OP>
00606   void MinresIter<ScalarType,MV,OP>::symOrtho( ScalarType a, ScalarType b,
00607                                                ScalarType *c, ScalarType *s, ScalarType *r
00608                                              )
00609   {
00610     const ScalarType one = SCT::one();
00611     const ScalarType zero = SCT::zero();
00612     const MagnitudeType m_zero = SMT::zero();
00613     const MagnitudeType absA = SCT::magnitude( a );
00614     const MagnitudeType absB = SCT::magnitude( b );
00615     if ( absB == m_zero ) {
00616         *s = zero;
00617         *r = absA;
00618         if ( absA == m_zero )
00619             *c = one;
00620         else
00621             *c = a / absA;
00622     } else if ( absA == m_zero ) {
00623         *c = zero;
00624         *s = b / absB;
00625         *r = absB;
00626     } else if ( absB >= absA ) { // && a!=0 && b!=0
00627         ScalarType tau = a / b;
00628         if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
00629             *s = -one / SCT::squareroot( one+tau*tau );
00630         else
00631             *s =  one / SCT::squareroot( one+tau*tau );
00632         *c = *s * tau;
00633         *r = b / *s;
00634     } else { // (absA > absB) && a!=0 && b!=0
00635         ScalarType tau = b / a;
00636         if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
00637             *c = -one / SCT::squareroot( one+tau*tau );
00638         else
00639             *c =  one / SCT::squareroot( one+tau*tau );
00640         *s = *c * tau;
00641         *r = a / *c;
00642     }
00643   }
00644 
00645 } // end Belos namespace
00646 
00647 #endif /* BELOS_MINRES_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines