Belos Package Browser (Single Doxygen Collection) Development
BelosTFQMRIter.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 // This file contains an implementation of the TFQMR iteration
00043 // for solving non-Hermitian linear systems of equations Ax = b, 
00044 // where b is a single-std::vector and x is the corresponding solution.
00045 //
00046 // The implementation is a slight modification on the TFQMR iteration
00047 // found in Saad's "Iterative Methods for Sparse Linear Systems".
00048 //
00049 
00050 #ifndef BELOS_TFQMR_ITER_HPP
00051 #define BELOS_TFQMR_ITER_HPP
00052 
00060 #include "BelosConfigDefs.hpp"
00061 #include "BelosIteration.hpp"
00062 #include "BelosTypes.hpp" 
00063 
00064 #include "BelosLinearProblem.hpp"
00065 #include "BelosMatOrthoManager.hpp"
00066 #include "BelosOutputManager.hpp"
00067 #include "BelosStatusTest.hpp"
00068 #include "BelosOperatorTraits.hpp"
00069 #include "BelosMultiVecTraits.hpp"
00070 
00071 #include "Teuchos_BLAS.hpp"
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 
00089 namespace Belos {
00090 
00095   template <class ScalarType, class MV>
00096   struct TFQMRIterState {
00097 
00099     Teuchos::RCP<const MV> R;
00100     Teuchos::RCP<const MV> W;
00101     Teuchos::RCP<const MV> U;
00102     Teuchos::RCP<const MV> Rtilde;
00103     Teuchos::RCP<const MV> D;
00104     Teuchos::RCP<const MV> V;
00105 
00106     TFQMRIterState() : R(Teuchos::null), W(Teuchos::null), U(Teuchos::null),
00107                        Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
00108     {}
00109   };
00110 
00111   
00113 
00114   
00126   class TFQMRIterInitFailure : public BelosError {public:
00127     TFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00128     {}};
00129   
00136   class TFQMRIterateFailure : public BelosError {public:
00137     TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
00138     {}};
00139   
00141 
00142 
00143   template <class ScalarType, class MV, class OP>
00144   class TFQMRIter : public Iteration<ScalarType,MV,OP> { 
00145   public:
00146     //
00147     // Convenience typedefs
00148     //
00149     typedef MultiVecTraits<ScalarType,MV> MVT;
00150     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00151     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00152     typedef typename SCT::magnitudeType MagnitudeType;
00153     
00155 
00156 
00158     TFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00159          const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00160          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00161          Teuchos::ParameterList &params );
00162     
00164     virtual ~TFQMRIter() {};
00166   
00167 
00169 
00170     
00192     void iterate();
00193 
00215     void initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate);
00216     
00220     void initialize()
00221     {
00222       TFQMRIterState<ScalarType,MV> empty;
00223       initializeTFQMR(empty);
00224     }
00225     
00233     TFQMRIterState<ScalarType,MV> getState() const {
00234       TFQMRIterState<ScalarType,MV> state;
00235       state.R = R_;
00236       state.W = W_;
00237       state.U = U_;
00238       state.Rtilde = Rtilde_;
00239       state.D = D_;
00240       state.V = V_;
00241       return state;
00242     }
00243     
00245     
00246     
00248 
00249     
00251     int getNumIters() const { return iter_; }
00252     
00254     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00255     
00258     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00259     
00261 
00263     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00264 
00266 
00267 
00269 
00270     
00272     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00273     
00275     int getBlockSize() const { return 1; }
00276     
00278     void setBlockSize(int blockSize) {
00279       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00280                        "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
00281     }
00282     
00284     bool isInitialized() { return initialized_; }
00285     
00287     
00288   
00289   private:
00290 
00291     //
00292     // Internal methods
00293     //
00295     void setStateSize();
00296     
00297     //
00298     // Classes inputed through constructor that define the linear problem to be solved.
00299     //
00300     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00301     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00302     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00303     
00304     //
00305     // Algorithmic parameters
00306     //      
00307 
00308     // Storage for QR factorization of the least squares system.
00309     Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
00310     std::vector<MagnitudeType> tau_, cs_, theta_;
00311     
00312     // 
00313     // Current solver state
00314     //
00315     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00316     // is capable of running; _initialize is controlled  by the initialize() member method
00317     // For the implications of the state of initialized_, please see documentation for initialize()
00318     bool initialized_;
00319     
00320     // stateStorageInitialized_ specifies that the state storage has be initialized to the current
00321     // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00322     // generated without the right-hand side or solution vectors.
00323     bool stateStorageInitialized_;
00324     
00325      // Current subspace dimension, and number of iterations performed.
00326     int iter_;
00327     
00328     // 
00329     // State Storage
00330     //
00331     Teuchos::RCP<MV> R_;
00332     Teuchos::RCP<MV> W_;
00333     Teuchos::RCP<MV> U_, AU_;
00334     Teuchos::RCP<MV> Rtilde_;
00335     Teuchos::RCP<MV> D_;
00336     Teuchos::RCP<MV> V_;
00337 
00338   };
00339   
00340   
00341   //
00342   // Implementation
00343   //
00344   
00346   // Constructor.
00347   template <class ScalarType, class MV, class OP>
00348   TFQMRIter<ScalarType,MV,OP>::TFQMRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00349            const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00350            const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00351            Teuchos::ParameterList &params 
00352            ) : 
00353     lp_(problem), 
00354     om_(printer),
00355     stest_(tester),
00356     alpha_(1,1),
00357     rho_(1,1),
00358     rho_old_(1,1),
00359     tau_(1),
00360     cs_(1),
00361     theta_(1),
00362     initialized_(false),
00363     stateStorageInitialized_(false),
00364     iter_(0)
00365   { 
00366   }
00367 
00369   // Compute native residual from TFQMR recurrence.
00370   template <class ScalarType, class MV, class OP>
00371   Teuchos::RCP<const MV> 
00372   TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00373   {
00374     MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00375     if (normvec)
00376       (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( iter_ + one )*tau_[0];
00377 
00378     return Teuchos::null;
00379   }
00380   
00381 
00383   // Setup the state storage.
00384   template <class ScalarType, class MV, class OP>
00385   void TFQMRIter<ScalarType,MV,OP>::setStateSize ()
00386   {
00387     if (!stateStorageInitialized_) {
00388 
00389       // Check if there is any multivector to clone from.
00390       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00391       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00392       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00393         stateStorageInitialized_ = false;
00394         return;
00395       }
00396       else {
00397 
00398         // Initialize the state storage
00399         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00400         if (R_ == Teuchos::null) {
00401           // Get the multivector that is not null.
00402           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00403           TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00404                              "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00405           R_ = MVT::Clone( *tmp, 1 );
00406     AU_ = MVT::Clone( *tmp, 1 );
00407     D_ = MVT::Clone( *tmp, 1 );
00408           V_ = MVT::Clone( *tmp, 1 );
00409         }
00410 
00411         // State storage has now been initialized.
00412         stateStorageInitialized_ = true;
00413       }
00414     }
00415   }
00416 
00418   // Initialize this iteration object
00419   template <class ScalarType, class MV, class OP>
00420   void TFQMRIter<ScalarType,MV,OP>::initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate)
00421   {
00422     // Initialize the state storage if it isn't already.
00423     if (!stateStorageInitialized_)
00424       setStateSize();
00425 
00426     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00427                        "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
00428 
00429     // NOTE:  In TFQMRIter R_, the initial residual, is required!!!
00430     //
00431     std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
00432 
00433     // Create convenience variables for zero and one.
00434     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00435     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00436     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00437 
00438     if (newstate.R != Teuchos::null) {
00439 
00440       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00441                           std::invalid_argument, errstr );
00442       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00443                           std::invalid_argument, errstr );
00444 
00445       // Copy basis vectors from newstate into V
00446       if (newstate.R != R_) {
00447         // copy over the initial residual (unpreconditioned).
00448         MVT::MvAddMv( one, *newstate.R, STzero, *newstate.R, *R_ );
00449       }
00450 
00451       // Compute initial vectors
00452       // Initially, they are set to the preconditioned residuals
00453       //
00454       W_ = MVT::CloneCopy( *R_ );
00455       U_ = MVT::CloneCopy( *R_ );
00456       Rtilde_ = MVT::CloneCopy( *R_ );
00457       MVT::MvInit( *D_ );
00458       // Multiply the current residual by Op and store in V_
00459       //       V_ = Op * R_ 
00460       //
00461       lp_->apply( *U_, *V_ );
00462       AU_ = MVT::CloneCopy( *V_ ); 
00463       //
00464       // Compute initial scalars: theta, eta, tau, rho_old
00465       //
00466       theta_[0] = MTzero;
00467       MVT::MvNorm( *R_, tau_ );                         // tau = ||r_0||
00468       MVT::MvTransMv( one, *Rtilde_, *R_, rho_old_ );   // rho = (r_0, r_tilde)
00469     }
00470     else {
00471 
00472       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00473                          "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
00474     }
00475 
00476     // The solver is initialized
00477     initialized_ = true;
00478   }
00479 
00480 
00482   // Iterate until the status test informs us we should stop.
00483   template <class ScalarType, class MV, class OP>
00484   void TFQMRIter<ScalarType,MV,OP>::iterate() 
00485   {
00486     //
00487     // Allocate/initialize data structures
00488     //
00489     if (initialized_ == false) {
00490       initialize();
00491     }
00492  
00493     // Create convenience variables for zero and one.
00494     const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
00495     const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00496     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00497     ScalarType eta = STzero, beta = STzero;
00498     //
00499     //  Start executable statements. 
00500     //
00501     // Get the current solution vector.
00502     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00503 
00504     // Check that the current solution vector only has one column.
00505     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
00506                         "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
00507 
00508 
00510     // Iterate until the status test tells us to stop.
00511     //
00512     while (stest_->checkStatus(this) != Passed) {
00513 
00514       //
00515       //--------------------------------------------------------
00516       // Compute the new alpha if we need to
00517       //--------------------------------------------------------
00518       //
00519       if (iter_%2 == 0) {
00520   MVT::MvTransMv( STone, *V_, *Rtilde_, alpha_ );      //   alpha = rho / (v, r_tilde) 
00521   alpha_(0,0) = rho_old_(0,0)/alpha_(0,0);
00522       }
00523       //
00524       //--------------------------------------------------------
00525       // Update w.
00526       //   w = w - alpha*Au
00527       //--------------------------------------------------------
00528       //
00529       MVT::MvAddMv( STone, *W_, -alpha_(0,0), *AU_, *W_ );
00530       //
00531       //--------------------------------------------------------
00532       // Update d.
00533       //   d = u + (theta^2/alpha)eta*d
00534       //--------------------------------------------------------
00535       //
00536       MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_(0,0))*eta, *D_, *D_ );
00537       //
00538       //--------------------------------------------------------
00539       // Update u if we need to.
00540       //   u = u - alpha*v
00541       //   
00542       // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
00543       //--------------------------------------------------------
00544       //
00545       if (iter_%2 == 0) {
00546         // Compute new U.
00547   MVT::MvAddMv( STone, *U_, -alpha_(0,0), *V_, *U_ );
00548 
00549   // Update Au for the next iteration.
00550   lp_->apply( *U_, *AU_ );                       
00551       }
00552       //
00553       //--------------------------------------------------------
00554       // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
00555       //--------------------------------------------------------
00556       //
00557       MVT::MvNorm( *W_, theta_ );     // theta = ||w|| / tau
00558       theta_[0] /= tau_[0];
00559       // cs = 1.0 / sqrt(1.0 + theta^2)
00560       cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);  
00561       tau_[0] *= theta_[0]*cs_[0];     // tau = tau * theta * cs
00562       eta = cs_[0]*cs_[0]*alpha_(0,0);     // eta = cs^2 * alpha
00563       
00564       //
00565       //--------------------------------------------------------
00566       // Update the solution.
00567       //--------------------------------------------------------
00568       //
00569       lp_->updateSolution( D_, true, eta );
00570       //
00571       if (iter_%2) {
00572   //
00573   //--------------------------------------------------------
00574   // Compute the new rho, beta if we need to.
00575   //--------------------------------------------------------
00576   //
00577   MVT::MvTransMv( STone, *W_, *Rtilde_, rho_ );     // rho = ( w, r_tilde )
00578   beta = rho_(0,0)/rho_old_(0,0);                   // beta = rho / rho_old
00579   rho_old_(0,0) = rho_(0,0);                        // rho_old = rho
00580   //
00581   //--------------------------------------------------------
00582   // Update u, v, and Au if we need to.
00583   // Note: We are updating v in two stages to be memory efficient
00584   //--------------------------------------------------------
00585   //
00586   MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );       // u = w + beta*u
00587   
00588   // First stage of v update.
00589   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v 
00590   
00591   // Update Au.
00592   lp_->apply( *U_, *AU_ );                          // Au = A*u
00593   
00594   // Second stage of v update.
00595   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v
00596       }
00597 
00598       // Increment the iteration
00599       iter_++;
00600       
00601     } // end while (sTest_->checkStatus(this) != Passed)
00602   } 
00603 
00604 } // namespace Belos
00605 //
00606 #endif // BELOS_TFQMR_ITER_HPP
00607 //
00608 // End of file BelosTFQMRIter.hpp
00609 
00610 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines