BelosTFQMR.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 algorithm
00030 // for solving non-Hermitian linear systems of equations Ax = b, 
00031 // where b is a single-vector and x is the corresponding solution.
00032 //
00033 // The implementation is a slight modification on the TFQMR algorithm
00034 // found in Saad's "Iterative Methods for Sparse Linear Systems".
00035 //
00036 
00037 #ifndef BELOS_TFQMR_HPP
00038 #define BELOS_TFQMR_HPP
00039 
00047 #include "BelosConfigDefs.hpp"
00048 #include "BelosIterativeSolver.hpp"
00049 #include "BelosLinearProblem.hpp"
00050 #include "BelosOutputManager.hpp"
00051 #include "BelosOperator.hpp"
00052 #include "BelosStatusTest.hpp"
00053 #include "BelosMultiVecTraits.hpp"
00054 
00055 #include "Teuchos_LAPACK.hpp"
00056 #include "Teuchos_SerialDenseMatrix.hpp"
00057 #include "Teuchos_ScalarTraits.hpp"
00058 #include "Teuchos_TimeMonitor.hpp"
00059 
00069 namespace Belos {
00070   
00071   template <class ScalarType, class MV, class OP>
00072   class TFQMR : public IterativeSolver<ScalarType,MV,OP> { 
00073   public:
00074     //
00075     // Convenience typedefs
00076     //
00077     typedef MultiVecTraits<ScalarType,MV> MVT;
00078     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00079     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00080     typedef typename SCT::magnitudeType MagnitudeType;
00081     
00083 
00084 
00086     TFQMR(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp, 
00087     const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest, 
00088     const RefCountPtr<OutputManager<ScalarType> > &om,
00089     const RefCountPtr<ParameterList> &pl = Teuchos::null
00090     );
00091     
00093     virtual ~TFQMR() {};
00095     
00097 
00098     
00100     int GetNumIters() const { return( _iter ); }
00101     
00103     int GetNumRestarts() const { return(0); }
00104     
00106 
00109     RefCountPtr<const MV> GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const;
00110     
00112 
00116     RefCountPtr<MV> GetCurrentSoln() { return MVT::CloneCopy( *_cur_block_sol ); };
00117   
00119 
00121     RefCountPtr<LinearProblem<ScalarType,MV,OP> > GetLinearProblem() const { return( _lp ); }
00122     
00123     RefCountPtr<StatusTest<ScalarType,MV,OP> > GetStatusTest() const { return( _stest ); }
00124 
00126     
00128 
00129     
00134     void Solve();
00136 
00139     
00141     std::string description() const;
00142     
00144     
00145   private:
00146     
00148     RefCountPtr<LinearProblem<ScalarType,MV,OP> > _lp; 
00149     
00151     RefCountPtr<StatusTest<ScalarType,MV,OP> > _stest; 
00152     
00154     RefCountPtr<OutputManager<ScalarType> > _om;
00155     
00157     RefCountPtr<ParameterList> _pl;     
00158 
00160     RefCountPtr<MV> _cur_block_sol;
00161     
00163     RefCountPtr<MV> _cur_block_rhs; 
00164     
00166     RefCountPtr<MV> _residvec;
00167     
00169     RefCountPtr<ostream> _os;
00170     
00172     int _iter;
00173 
00175     std::vector<MagnitudeType> _tau;
00176 
00178     bool _restartTimers;
00179 
00181     Teuchos::RefCountPtr<Teuchos::Time> _timerOp, _timerTotal;
00182   };
00183   
00184   //
00185   // Implementation
00186   //
00187   
00188   template <class ScalarType, class MV, class OP>
00189   TFQMR<ScalarType,MV,OP>::TFQMR(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp,
00190          const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest,
00191          const RefCountPtr<OutputManager<ScalarType> > &om,
00192          const RefCountPtr<ParameterList> &pl 
00193          ) : 
00194     _lp(lp), 
00195     _stest(stest),
00196     _om(om),
00197     _pl(pl),
00198     _os(om->GetOStream()),
00199     _iter(0),
00200     _tau(1),
00201     _restartTimers(true),
00202     _timerOp(Teuchos::TimeMonitor::getNewTimer("Operation Op*x")),
00203     _timerTotal(Teuchos::TimeMonitor::getNewTimer("Total time"))
00204   { 
00205   }
00206   
00207   template <class ScalarType, class MV, class OP>
00208   RefCountPtr<const MV> 
00209   TFQMR<ScalarType,MV,OP>::GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00210   {
00211     MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
00212     if (normvec)
00213       (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( _iter + one )*_tau[0];
00214 
00215     return Teuchos::null;
00216   }
00217   
00218   template <class ScalarType, class MV, class OP>
00219   void 
00220   TFQMR<ScalarType,MV,OP>::Solve () 
00221   {
00222     Teuchos::TimeMonitor LocalTimer(*_timerTotal,_restartTimers);
00223     
00224     if ( _restartTimers ) {
00225       _timerOp->reset();
00226     }
00227     
00228     // Get the output stream from the OutputManager
00229     _os = _om->GetOStream();
00230     //
00231     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00232     const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
00233     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00234     bool exit_flg = false;
00235     Teuchos::SerialDenseMatrix<int,ScalarType> _alpha( 1, 1 );
00236     Teuchos::SerialDenseMatrix<int,ScalarType> _rho( 1, 1 ), _rho_old( 1, 1 );
00237     RefCountPtr<MV> _v, _w, _u, _Au, _d, _rtilde;
00238     ScalarType _eta = STzero, _beta = STzero;
00239     std::vector<MagnitudeType> _cs(1,MTzero), _theta(1,MTzero);
00240     //
00241     // Retrieve the first linear system to be solved.
00242     //
00243     _cur_block_sol = _lp->GetCurrLHSVec();
00244     _cur_block_rhs = _lp->GetCurrRHSVec();
00245     //
00246     //  Start executable statements. 
00247     //
00248     while (_cur_block_sol.get() && _cur_block_rhs.get() ) {
00249       //
00250       // Only continue if the linear system is single-vector.
00251       //
00252       if ( _lp->GetBlockSize() > 1 ) return;
00253       //
00254       if (_om->isVerbosityAndPrint( IterationDetails )) {
00255   *_os << endl;
00256   *_os << "===================================================" << endl;
00257   *_os << "Solving linear system(s):  " << _lp->GetRHSIndex() << " through " << _lp->GetRHSIndex()+_lp->GetNumToSolve() << endl;
00258   *_os << endl;
00259       } 
00260       //
00261       _d = MVT::Clone( *_cur_block_sol, 1 );
00262       _v = MVT::Clone( *_cur_block_sol, 1 );
00263       MVT::MvInit( *_d );
00264       //
00265       // ************ Compute the initial residual ********************************
00266       //
00267       // w0 = y0 = r0 = cur_block_rhs - A * cur_block_sol 
00268       //
00269       _residvec = MVT::Clone( *_cur_block_sol, 1 );
00270       _lp->ComputeResVec( &*_residvec, &*_cur_block_sol, &*_cur_block_rhs );
00271       //
00272       _w = MVT::CloneCopy( *_residvec ); 
00273       _u = MVT::CloneCopy( *_residvec ); 
00274       _rtilde = MVT::CloneCopy( *_residvec ); 
00275       //
00276       // Multiply the current residual by Op and store in _v
00277       //       _v = _Op*_residvec 
00278       //
00279       {
00280   Teuchos::TimeMonitor OpTimer(*_timerOp);
00281   _lp->Apply( *_residvec, *_v );
00282       }
00283       _Au = MVT::CloneCopy( *_v ); 
00284       //
00285       // Compute initial scalars: theta, eta, tau, rho_old
00286       //
00287       MVT::MvNorm( *_residvec, &_tau );                         // tau = ||r_0||
00288       MVT::MvTransMv( one, *_residvec, *_rtilde, _rho_old );    // rho = (r_0, r_tilde)
00289       //
00290       // ***************************************************************************
00291       // ************************Main TFQMR Loop***************************************
00292       // ***************************************************************************
00293       // 
00294       for (_iter=0; _stest->CheckStatus(this) == Unconverged && !exit_flg; _iter++) 
00295   {
00296     //
00297     //--------------------------------------------------------
00298     // Compute the new alpha if we need to
00299     //--------------------------------------------------------
00300     //
00301     if (_iter%2 == 0) {
00302       MVT::MvTransMv( one, *_v, *_rtilde, _alpha );      //   alpha = rho / (v, r_tilde) 
00303       _alpha(0,0) = _rho_old(0,0)/_alpha(0,0);
00304     }
00305     //
00306     //--------------------------------------------------------
00307     // Update d.
00308     //   d = u + (theta^2/alpha)eta*d
00309     //--------------------------------------------------------
00310     //
00311     MVT::MvAddMv( one, *_u, (_theta[0]*_theta[0]/_alpha(0,0))*_eta, *_d, *_d );
00312     //
00313     //--------------------------------------------------------
00314     // Update w.
00315     //   w = w - alpha*Au
00316     //--------------------------------------------------------
00317     //
00318     MVT::MvAddMv( one, *_w, -_alpha(0,0), *_Au, *_w );
00319     //
00320     //--------------------------------------------------------
00321     // Update u if we need to.
00322     //   u = u - alpha*v
00323     //   
00324     // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
00325     //--------------------------------------------------------
00326     //
00327     if (_iter%2 == 0) {
00328       MVT::MvAddMv( one, *_u, -_alpha(0,0), *_v, *_u );
00329     }
00330     //
00331     //--------------------------------------------------------
00332     // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
00333     //--------------------------------------------------------
00334     //
00335     MVT::MvNorm( *_w, &_theta );     // theta = ||w|| / tau
00336     _theta[0] /= _tau[0];
00337 
00338     // cs = sqrt( 1 + theta^2 )
00339     _cs[0] = Teuchos::ScalarTraits<ScalarType>::squareroot(one + _theta[0]*_theta[0]);
00340 
00341     _tau[0] *= _theta[0]*_cs[0];     // tau = tau * theta * cs
00342     _eta = _cs[0]*_cs[0]*_alpha(0,0);     // eta = cs^2 * alpha
00343 
00344     //
00345     //--------------------------------------------------------
00346     // Update the solution.
00347     //--------------------------------------------------------
00348     //
00349     _lp->SolutionUpdated( &*_d, _eta );
00350     //
00351     if (_iter%2) {
00352       //
00353       //--------------------------------------------------------
00354       // Compute the new rho, beta if we need to.
00355       //--------------------------------------------------------
00356       //
00357       MVT::MvTransMv( one, *_w, *_rtilde, _rho );       // rho = ( w, r_tilde )
00358       _beta = _rho(0,0)/_rho_old(0,0);                  // beta = rho / rho_old
00359       _rho_old(0,0) = _rho(0,0);                        // rho_old = rho
00360       //
00361       //--------------------------------------------------------
00362       // Update u, v, and Au if we need to.
00363       // Note: We are updating v in two stages to be memory efficient
00364       //--------------------------------------------------------
00365       //
00366       MVT::MvAddMv( one, *_w, _beta, *_u, *_u );       // u = w + beta*u
00367 
00368       // First stage of v update.
00369       MVT::MvAddMv( one, *_Au, _beta, *_v, *_v );      // v = Au + beta*v 
00370 
00371       // Update Au.
00372       {
00373         Teuchos::TimeMonitor OpTimer(*_timerOp);
00374         _lp->Apply( *_u, *_Au );                       // Au = A*u
00375       }
00376 
00377       // Second stage of v update.
00378       MVT::MvAddMv( one, *_Au, _beta, *_v, *_v );      // v = Au + beta*v
00379     }
00380 
00381   } // end of the main TFQMR loop -- for(_iter = 0;...)
00382       // *******************************************************************************
00383       //
00384       // Inform the linear problem manager that we are done with the current block of linear systems.
00385       //
00386       _lp->SetCurrLSVec();
00387       //
00388       // Get the next block of linear systems, if it returns the null pointer we are done.
00389       //
00390       _cur_block_sol = _lp->GetCurrLHSVec();
00391       _cur_block_rhs = _lp->GetCurrRHSVec();
00392       //
00393       // Print out solver status.
00394       //
00395       if (_om->isVerbosityAndPrint( FinalSummary )) {
00396   _stest->Print(*_os);
00397       }
00398       //
00399     } // end while ( _cur_block_sol && _cur_block_rhs )
00400     // **********************************************************************************
00401     //
00402     // Print timing details 
00403 
00404     // Stop timer.
00405     _timerTotal->stop();
00406 
00407     // Reset format that will be used to print the summary
00408     Teuchos::TimeMonitor::format().setPageWidth(54);
00409 
00410     if (_om->isVerbosity( Belos::TimingDetails )) {
00411       if (_om->doPrint())
00412         *_os <<"********************TIMING DETAILS********************"<<endl;
00413       Teuchos::TimeMonitor::summarize( *_os );
00414       if (_om->doPrint())
00415         *_os <<"******************************************************"<<endl;
00416     }
00417   } // end TFQMRSolve()
00418 
00419 
00420   // Overridden from Teuchos::Describable
00421   template<class ScalarType, class MV, class OP>
00422   std::string 
00423   TFQMR<ScalarType,MV,OP>::description() const
00424   {
00425     std::ostringstream oss;
00426     oss << "Belos::TFQMR<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00427     oss << "{";
00428     oss << "}";
00429     return oss.str();
00430   }
00431 
00432   //
00433 } // namespace Belos
00434 //
00435 #endif // BELOS_TFQMR_HPP
00436 //
00437 // End of file BelosTFQMR.hpp
00438 
00439 

Generated on Thu Sep 18 12:30:12 2008 for Belos by doxygen 1.3.9.1