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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 // This file contains an implementation of the TFQMR iteration
00030 // for solving non-Hermitian linear systems of equations Ax = b, 
00031 // where b is a single-std::vector and x is the corresponding solution.
00032 //
00033 // The implementation is a slight modification on the TFQMR iteration
00034 // found in Saad's "Iterative Methods for Sparse Linear Systems".
00035 //
00036 
00037 #ifndef BELOS_TFQMR_ITER_HPP
00038 #define BELOS_TFQMR_ITER_HPP
00039 
00047 #include "BelosConfigDefs.hpp"
00048 #include "BelosIteration.hpp"
00049 #include "BelosTypes.hpp" 
00050 
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosMatOrthoManager.hpp"
00053 #include "BelosOutputManager.hpp"
00054 #include "BelosStatusTest.hpp"
00055 #include "BelosOperatorTraits.hpp"
00056 #include "BelosMultiVecTraits.hpp"
00057 
00058 #include "Teuchos_BLAS.hpp"
00059 #include "Teuchos_SerialDenseMatrix.hpp"
00060 #include "Teuchos_SerialDenseVector.hpp"
00061 #include "Teuchos_ScalarTraits.hpp"
00062 #include "Teuchos_ParameterList.hpp"
00063 #include "Teuchos_TimeMonitor.hpp"
00064 
00076 namespace Belos {
00077 
00082   template <class ScalarType, class MV>
00083   struct TFQMRIterState {
00084 
00086     Teuchos::RCP<const MV> R;
00087     Teuchos::RCP<const MV> W;
00088     Teuchos::RCP<const MV> U;
00089     Teuchos::RCP<const MV> Rtilde;
00090     Teuchos::RCP<const MV> D;
00091     Teuchos::RCP<const MV> V;
00092 
00093     TFQMRIterState() : R(Teuchos::null), W(Teuchos::null), U(Teuchos::null),
00094                        Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
00095     {}
00096   };
00097 
00098   
00100 
00101   
00113   class TFQMRIterInitFailure : public BelosError {public:
00114     TFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
00115     {}};
00116   
00123   class TFQMRIterateFailure : public BelosError {public:
00124     TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
00125     {}};
00126   
00128 
00129 
00130   template <class ScalarType, class MV, class OP>
00131   class TFQMRIter : public Iteration<ScalarType,MV,OP> { 
00132   public:
00133     //
00134     // Convenience typedefs
00135     //
00136     typedef MultiVecTraits<ScalarType,MV> MVT;
00137     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00138     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00139     typedef typename SCT::magnitudeType MagnitudeType;
00140     
00142 
00143 
00145     TFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00146          const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00147          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00148          Teuchos::ParameterList &params );
00149     
00151     virtual ~TFQMRIter() {};
00153   
00154 
00156 
00157     
00179     void iterate();
00180 
00202     void initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate);
00203     
00207     void initialize()
00208     {
00209       TFQMRIterState<ScalarType,MV> empty;
00210       initializeTFQMR(empty);
00211     }
00212     
00220     TFQMRIterState<ScalarType,MV> getState() const {
00221       TFQMRIterState<ScalarType,MV> state;
00222       state.R = R_;
00223       state.W = W_;
00224       state.U = U_;
00225       state.Rtilde = Rtilde_;
00226       state.D = D_;
00227       state.V = V_;
00228       return state;
00229     }
00230     
00232     
00233     
00235 
00236     
00238     int getNumIters() const { return iter_; }
00239     
00241     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00242     
00245     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00246     
00248 
00250     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00251 
00253 
00254 
00256 
00257     
00259     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00260     
00262     int getBlockSize() const { return 1; }
00263     
00265     void setBlockSize(int blockSize) {
00266       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00267                        "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
00268     }
00269     
00271     bool isInitialized() { return initialized_; }
00272     
00274     
00275   
00276   private:
00277 
00278     //
00279     // Internal methods
00280     //
00282     void setStateSize();
00283     
00284     //
00285     // Classes inputed through constructor that define the linear problem to be solved.
00286     //
00287     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00288     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00289     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00290     
00291     //
00292     // Algorithmic parameters
00293     //      
00294 
00295     // Storage for QR factorization of the least squares system.
00296     Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
00297     std::vector<MagnitudeType> tau_, cs_, theta_;
00298     
00299     // 
00300     // Current solver state
00301     //
00302     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00303     // is capable of running; _initialize is controlled  by the initialize() member method
00304     // For the implications of the state of initialized_, please see documentation for initialize()
00305     bool initialized_;
00306     
00307     // stateStorageInitialized_ specifies that the state storage has be initialized to the current
00308     // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00309     // generated without the right-hand side or solution vectors.
00310     bool stateStorageInitialized_;
00311     
00312      // Current subspace dimension, and number of iterations performed.
00313     int iter_;
00314     
00315     // 
00316     // State Storage
00317     //
00318     Teuchos::RCP<MV> R_;
00319     Teuchos::RCP<MV> W_;
00320     Teuchos::RCP<MV> U_, AU_;
00321     Teuchos::RCP<MV> Rtilde_;
00322     Teuchos::RCP<MV> D_;
00323     Teuchos::RCP<MV> V_;
00324 
00325   };
00326   
00327   
00328   //
00329   // Implementation
00330   //
00331   
00333   // Constructor.
00334   template <class ScalarType, class MV, class OP>
00335   TFQMRIter<ScalarType,MV,OP>::TFQMRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00336            const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00337            const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00338            Teuchos::ParameterList &params 
00339            ) : 
00340     lp_(problem), 
00341     om_(printer),
00342     stest_(tester),
00343     alpha_(1,1),
00344     rho_(1,1),
00345     rho_old_(1,1),
00346     tau_(1),
00347     cs_(1),
00348     theta_(1),
00349     initialized_(false),
00350     stateStorageInitialized_(false),
00351     iter_(0)
00352   { 
00353   }
00354 
00356   // Compute native residual from TFQMR recurrence.
00357   template <class ScalarType, class MV, class OP>
00358   Teuchos::RCP<const MV> 
00359   TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00360   {
00361     MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00362     if (normvec)
00363       (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( iter_ + one )*tau_[0];
00364 
00365     return Teuchos::null;
00366   }
00367   
00368 
00370   // Setup the state storage.
00371   template <class ScalarType, class MV, class OP>
00372   void TFQMRIter<ScalarType,MV,OP>::setStateSize ()
00373   {
00374     if (!stateStorageInitialized_) {
00375 
00376       // Check if there is any multivector to clone from.
00377       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00378       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00379       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00380         stateStorageInitialized_ = false;
00381         return;
00382       }
00383       else {
00384 
00385         // Initialize the state storage
00386         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00387         if (R_ == Teuchos::null) {
00388           // Get the multivector that is not null.
00389           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00390           TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00391                              "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00392           R_ = MVT::Clone( *tmp, 1 );
00393     AU_ = MVT::Clone( *tmp, 1 );
00394     D_ = MVT::Clone( *tmp, 1 );
00395           V_ = MVT::Clone( *tmp, 1 );
00396         }
00397 
00398         // State storage has now been initialized.
00399         stateStorageInitialized_ = true;
00400       }
00401     }
00402   }
00403 
00405   // Initialize this iteration object
00406   template <class ScalarType, class MV, class OP>
00407   void TFQMRIter<ScalarType,MV,OP>::initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate)
00408   {
00409     // Initialize the state storage if it isn't already.
00410     if (!stateStorageInitialized_)
00411       setStateSize();
00412 
00413     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00414                        "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
00415 
00416     // NOTE:  In TFQMRIter R_, the initial residual, is required!!!
00417     //
00418     std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
00419 
00420     // Create convenience variables for zero and one.
00421     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00422     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00423     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00424 
00425     if (newstate.R != Teuchos::null) {
00426 
00427       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00428                           std::invalid_argument, errstr );
00429       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00430                           std::invalid_argument, errstr );
00431 
00432       // Copy basis vectors from newstate into V
00433       if (newstate.R != R_) {
00434         // copy over the initial residual (unpreconditioned).
00435         MVT::MvAddMv( one, *newstate.R, STzero, *newstate.R, *R_ );
00436       }
00437 
00438       // Compute initial vectors
00439       // Initially, they are set to the preconditioned residuals
00440       //
00441       W_ = MVT::CloneCopy( *R_ );
00442       U_ = MVT::CloneCopy( *R_ );
00443       Rtilde_ = MVT::CloneCopy( *R_ );
00444       MVT::MvInit( *D_ );
00445       // Multiply the current residual by Op and store in V_
00446       //       V_ = Op * R_ 
00447       //
00448       lp_->apply( *U_, *V_ );
00449       AU_ = MVT::CloneCopy( *V_ ); 
00450       //
00451       // Compute initial scalars: theta, eta, tau, rho_old
00452       //
00453       theta_[0] = MTzero;
00454       MVT::MvNorm( *R_, tau_ );                         // tau = ||r_0||
00455       MVT::MvTransMv( one, *Rtilde_, *R_, rho_old_ );   // rho = (r_0, r_tilde)
00456     }
00457     else {
00458 
00459       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00460                          "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
00461     }
00462 
00463     // The solver is initialized
00464     initialized_ = true;
00465   }
00466 
00467 
00469   // Iterate until the status test informs us we should stop.
00470   template <class ScalarType, class MV, class OP>
00471   void TFQMRIter<ScalarType,MV,OP>::iterate() 
00472   {
00473     //
00474     // Allocate/initialize data structures
00475     //
00476     if (initialized_ == false) {
00477       initialize();
00478     }
00479  
00480     // Create convenience variables for zero and one.
00481     const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
00482     const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00483     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00484     ScalarType eta = STzero, beta = STzero;
00485     //
00486     //  Start executable statements. 
00487     //
00488     // Get the current solution vector.
00489     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00490 
00491     // Check that the current solution vector only has one column.
00492     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
00493                         "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
00494 
00495 
00497     // Iterate until the status test tells us to stop.
00498     //
00499     while (stest_->checkStatus(this) != Passed) {
00500 
00501       //
00502       //--------------------------------------------------------
00503       // Compute the new alpha if we need to
00504       //--------------------------------------------------------
00505       //
00506       if (iter_%2 == 0) {
00507   MVT::MvTransMv( STone, *V_, *Rtilde_, alpha_ );      //   alpha = rho / (v, r_tilde) 
00508   alpha_(0,0) = rho_old_(0,0)/alpha_(0,0);
00509       }
00510       //
00511       //--------------------------------------------------------
00512       // Update w.
00513       //   w = w - alpha*Au
00514       //--------------------------------------------------------
00515       //
00516       MVT::MvAddMv( STone, *W_, -alpha_(0,0), *AU_, *W_ );
00517       //
00518       //--------------------------------------------------------
00519       // Update d.
00520       //   d = u + (theta^2/alpha)eta*d
00521       //--------------------------------------------------------
00522       //
00523       MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_(0,0))*eta, *D_, *D_ );
00524       //
00525       //--------------------------------------------------------
00526       // Update u if we need to.
00527       //   u = u - alpha*v
00528       //   
00529       // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
00530       //--------------------------------------------------------
00531       //
00532       if (iter_%2 == 0) {
00533         // Compute new U.
00534   MVT::MvAddMv( STone, *U_, -alpha_(0,0), *V_, *U_ );
00535 
00536   // Update Au for the next iteration.
00537   lp_->apply( *U_, *AU_ );                       
00538       }
00539       //
00540       //--------------------------------------------------------
00541       // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
00542       //--------------------------------------------------------
00543       //
00544       MVT::MvNorm( *W_, theta_ );     // theta = ||w|| / tau
00545       theta_[0] /= tau_[0];
00546       cs_[0] = MTone / SCT::squareroot(STone + theta_[0]*theta_[0]);  // cs = 1.0 / sqrt(1.0 + theta^2)
00547       tau_[0] *= theta_[0]*cs_[0];     // tau = tau * theta * cs
00548       eta = cs_[0]*cs_[0]*alpha_(0,0);     // eta = cs^2 * alpha
00549       
00550       //
00551       //--------------------------------------------------------
00552       // Update the solution.
00553       //--------------------------------------------------------
00554       //
00555       lp_->updateSolution( D_, true, eta );
00556       //
00557       if (iter_%2) {
00558   //
00559   //--------------------------------------------------------
00560   // Compute the new rho, beta if we need to.
00561   //--------------------------------------------------------
00562   //
00563   MVT::MvTransMv( STone, *W_, *Rtilde_, rho_ );     // rho = ( w, r_tilde )
00564   beta = rho_(0,0)/rho_old_(0,0);                   // beta = rho / rho_old
00565   rho_old_(0,0) = rho_(0,0);                        // rho_old = rho
00566   //
00567   //--------------------------------------------------------
00568   // Update u, v, and Au if we need to.
00569   // Note: We are updating v in two stages to be memory efficient
00570   //--------------------------------------------------------
00571   //
00572   MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );       // u = w + beta*u
00573   
00574   // First stage of v update.
00575   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v 
00576   
00577   // Update Au.
00578   lp_->apply( *U_, *AU_ );                          // Au = A*u
00579   
00580   // Second stage of v update.
00581   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v
00582       }
00583 
00584       // Increment the iteration
00585       iter_++;
00586       
00587     } // end while (sTest_->checkStatus(this) != Passed)
00588   } 
00589 
00590 } // namespace Belos
00591 //
00592 #endif // BELOS_TFQMR_ITER_HPP
00593 //
00594 // End of file BelosTFQMRIter.hpp
00595 
00596 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:07 2011 for Belos by  doxygen 1.6.3