Belos Version of the Day
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-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 MultiVecTraitsExt<ScalarType,MV> MVText;
00151     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00152     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00153     typedef typename SCT::magnitudeType MagnitudeType;
00154     
00156 
00157 
00159     TFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00160          const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00161          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00162          Teuchos::ParameterList &params );
00163     
00165     virtual ~TFQMRIter() {};
00167   
00168 
00170 
00171     
00193     void iterate();
00194 
00216     void initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate);
00217     
00221     void initialize()
00222     {
00223       TFQMRIterState<ScalarType,MV> empty;
00224       initializeTFQMR(empty);
00225     }
00226     
00234     TFQMRIterState<ScalarType,MV> getState() const {
00235       TFQMRIterState<ScalarType,MV> state;
00236       state.R = R_;
00237       state.W = W_;
00238       state.U = U_;
00239       state.Rtilde = Rtilde_;
00240       state.D = D_;
00241       state.V = V_;
00242       return state;
00243     }
00244     
00246     
00247     
00249 
00250     
00252     int getNumIters() const { return iter_; }
00253     
00255     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00256     
00259     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00260     
00262 
00264     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00265 
00267 
00268 
00270 
00271     
00273     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00274     
00276     int getBlockSize() const { return 1; }
00277     
00279     void setBlockSize(int blockSize) {
00280       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00281                        "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
00282     }
00283     
00285     bool isInitialized() { return initialized_; }
00286     
00288     
00289   
00290   private:
00291 
00292     //
00293     // Internal methods
00294     //
00296     void setStateSize();
00297     
00298     //
00299     // Classes inputed through constructor that define the linear problem to be solved.
00300     //
00301     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00302     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00303     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00304     
00305     //
00306     // Algorithmic parameters
00307     //      
00308 
00309     // Storage for QR factorization of the least squares system.
00310     Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
00311     std::vector<MagnitudeType> tau_, cs_, theta_;
00312     
00313     // 
00314     // Current solver state
00315     //
00316     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00317     // is capable of running; _initialize is controlled  by the initialize() member method
00318     // For the implications of the state of initialized_, please see documentation for initialize()
00319     bool initialized_;
00320     
00321     // stateStorageInitialized_ specifies that the state storage has be initialized to the current
00322     // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
00323     // generated without the right-hand side or solution vectors.
00324     bool stateStorageInitialized_;
00325     
00326      // Current subspace dimension, and number of iterations performed.
00327     int iter_;
00328     
00329     // 
00330     // State Storage
00331     //
00332     Teuchos::RCP<MV> R_;
00333     Teuchos::RCP<MV> W_;
00334     Teuchos::RCP<MV> U_, AU_;
00335     Teuchos::RCP<MV> Rtilde_;
00336     Teuchos::RCP<MV> D_;
00337     Teuchos::RCP<MV> V_;
00338 
00339   };
00340   
00341   
00342   //
00343   // Implementation
00344   //
00345   
00347   // Constructor.
00348   template <class ScalarType, class MV, class OP>
00349   TFQMRIter<ScalarType,MV,OP>::TFQMRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00350            const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00351            const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00352            Teuchos::ParameterList &params 
00353            ) : 
00354     lp_(problem), 
00355     om_(printer),
00356     stest_(tester),
00357     alpha_(1,1),
00358     rho_(1,1),
00359     rho_old_(1,1),
00360     tau_(1),
00361     cs_(1),
00362     theta_(1),
00363     initialized_(false),
00364     stateStorageInitialized_(false),
00365     iter_(0)
00366   { 
00367   }
00368 
00370   // Compute native residual from TFQMR recurrence.
00371   template <class ScalarType, class MV, class OP>
00372   Teuchos::RCP<const MV> 
00373   TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00374   {
00375     MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00376     if (normvec)
00377       (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( iter_ + one )*tau_[0];
00378 
00379     return Teuchos::null;
00380   }
00381   
00382 
00384   // Setup the state storage.
00385   template <class ScalarType, class MV, class OP>
00386   void TFQMRIter<ScalarType,MV,OP>::setStateSize ()
00387   {
00388     if (!stateStorageInitialized_) {
00389 
00390       // Check if there is any multivector to clone from.
00391       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00392       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00393       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00394         stateStorageInitialized_ = false;
00395         return;
00396       }
00397       else {
00398 
00399         // Initialize the state storage
00400         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00401         if (R_ == Teuchos::null) {
00402           // Get the multivector that is not null.
00403           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00404           TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00405                              "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00406           R_ = MVT::Clone( *tmp, 1 );
00407     AU_ = MVT::Clone( *tmp, 1 );
00408     D_ = MVT::Clone( *tmp, 1 );
00409           V_ = MVT::Clone( *tmp, 1 );
00410         }
00411 
00412         // State storage has now been initialized.
00413         stateStorageInitialized_ = true;
00414       }
00415     }
00416   }
00417 
00419   // Initialize this iteration object
00420   template <class ScalarType, class MV, class OP>
00421   void TFQMRIter<ScalarType,MV,OP>::initializeTFQMR(TFQMRIterState<ScalarType,MV> newstate)
00422   {
00423     // Initialize the state storage if it isn't already.
00424     if (!stateStorageInitialized_)
00425       setStateSize();
00426 
00427     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00428                        "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
00429 
00430     // NOTE:  In TFQMRIter R_, the initial residual, is required!!!
00431     //
00432     std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
00433 
00434     // Create convenience variables for zero and one.
00435     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00436     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00437     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00438 
00439     if (newstate.R != Teuchos::null) {
00440 
00441       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.R) != MVText::GetGlobalLength(*R_),
00442                           std::invalid_argument, errstr );
00443       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00444                           std::invalid_argument, errstr );
00445 
00446       // Copy basis vectors from newstate into V
00447       if (newstate.R != R_) {
00448         // copy over the initial residual (unpreconditioned).
00449         MVT::MvAddMv( one, *newstate.R, STzero, *newstate.R, *R_ );
00450       }
00451 
00452       // Compute initial vectors
00453       // Initially, they are set to the preconditioned residuals
00454       //
00455       W_ = MVT::CloneCopy( *R_ );
00456       U_ = MVT::CloneCopy( *R_ );
00457       Rtilde_ = MVT::CloneCopy( *R_ );
00458       MVT::MvInit( *D_ );
00459       // Multiply the current residual by Op and store in V_
00460       //       V_ = Op * R_ 
00461       //
00462       lp_->apply( *U_, *V_ );
00463       AU_ = MVT::CloneCopy( *V_ ); 
00464       //
00465       // Compute initial scalars: theta, eta, tau, rho_old
00466       //
00467       theta_[0] = MTzero;
00468       MVT::MvNorm( *R_, tau_ );                         // tau = ||r_0||
00469       MVT::MvTransMv( one, *Rtilde_, *R_, rho_old_ );   // rho = (r_tilde, r0)
00470     }
00471     else {
00472 
00473       TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00474                          "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
00475     }
00476 
00477     // The solver is initialized
00478     initialized_ = true;
00479   }
00480 
00481 
00483   // Iterate until the status test informs us we should stop.
00484   template <class ScalarType, class MV, class OP>
00485   void TFQMRIter<ScalarType,MV,OP>::iterate() 
00486   {
00487     //
00488     // Allocate/initialize data structures
00489     //
00490     if (initialized_ == false) {
00491       initialize();
00492     }
00493  
00494     // Create convenience variables for zero and one.
00495     const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
00496     const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
00497     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00498     ScalarType eta = STzero, beta = STzero;
00499     //
00500     //  Start executable statements. 
00501     //
00502     // Get the current solution vector.
00503     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00504 
00505     // Check that the current solution vector only has one column.
00506     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
00507                         "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
00508 
00509 
00511     // Iterate until the status test tells us to stop.
00512     //
00513     while (stest_->checkStatus(this) != Passed) {
00514 
00515       //
00516       //--------------------------------------------------------
00517       // Compute the new alpha if we need to
00518       //--------------------------------------------------------
00519       //
00520       if (iter_%2 == 0) {
00521   MVT::MvTransMv( STone, *Rtilde_, *V_, alpha_ );      //   alpha = rho / (r_tilde, v) 
00522   alpha_(0,0) = rho_old_(0,0)/alpha_(0,0);
00523       }
00524       //
00525       //--------------------------------------------------------
00526       // Update w.
00527       //   w = w - alpha*Au
00528       //--------------------------------------------------------
00529       //
00530       MVT::MvAddMv( STone, *W_, -alpha_(0,0), *AU_, *W_ );
00531       //
00532       //--------------------------------------------------------
00533       // Update d.
00534       //   d = u + (theta^2/alpha)eta*d
00535       //--------------------------------------------------------
00536       //
00537       MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_(0,0))*eta, *D_, *D_ );
00538       //
00539       //--------------------------------------------------------
00540       // Update u if we need to.
00541       //   u = u - alpha*v
00542       //   
00543       // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
00544       //--------------------------------------------------------
00545       //
00546       if (iter_%2 == 0) {
00547         // Compute new U.
00548   MVT::MvAddMv( STone, *U_, -alpha_(0,0), *V_, *U_ );
00549 
00550   // Update Au for the next iteration.
00551   lp_->apply( *U_, *AU_ );                       
00552       }
00553       //
00554       //--------------------------------------------------------
00555       // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
00556       //--------------------------------------------------------
00557       //
00558       MVT::MvNorm( *W_, theta_ );     // theta = ||w|| / tau
00559       theta_[0] /= tau_[0];
00560       // cs = 1.0 / sqrt(1.0 + theta^2)
00561       cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);  
00562       tau_[0] *= theta_[0]*cs_[0];     // tau = tau * theta * cs
00563       eta = cs_[0]*cs_[0]*alpha_(0,0);     // eta = cs^2 * alpha
00564       
00565       //
00566       //--------------------------------------------------------
00567       // Update the solution.
00568       //--------------------------------------------------------
00569       //
00570       lp_->updateSolution( D_, true, eta );
00571       //
00572       if (iter_%2) {
00573   //
00574   //--------------------------------------------------------
00575   // Compute the new rho, beta if we need to.
00576   //--------------------------------------------------------
00577   //
00578   MVT::MvTransMv( STone, *Rtilde_, *W_, rho_ );     // rho = (r_tilde, w)
00579   beta = rho_(0,0)/rho_old_(0,0);                   // beta = rho / rho_old
00580   rho_old_(0,0) = rho_(0,0);                        // rho_old = rho
00581   //
00582   //--------------------------------------------------------
00583   // Update u, v, and Au if we need to.
00584   // Note: We are updating v in two stages to be memory efficient
00585   //--------------------------------------------------------
00586   //
00587   MVT::MvAddMv( STone, *W_, beta, *U_, *U_ );       // u = w + beta*u
00588   
00589   // First stage of v update.
00590   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v 
00591   
00592   // Update Au.
00593   lp_->apply( *U_, *AU_ );                          // Au = A*u
00594   
00595   // Second stage of v update.
00596   MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ );      // v = Au + beta*v
00597       }
00598 
00599       // Increment the iteration
00600       iter_++;
00601       
00602     } // end while (sTest_->checkStatus(this) != Passed)
00603   } 
00604 
00605 } // namespace Belos
00606 //
00607 #endif // BELOS_TFQMR_ITER_HPP
00608 //
00609 // End of file BelosTFQMRIter.hpp
00610 
00611 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines