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