Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosTFQMRIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 //
10 // This file contains an implementation of the TFQMR iteration
11 // for solving non-Hermitian linear systems of equations Ax = b,
12 // where b is a single-vector and x is the corresponding solution.
13 //
14 // The implementation is a slight modification on the TFQMR iteration
15 // found in Saad's "Iterative Methods for Sparse Linear Systems".
16 //
17 
18 #ifndef BELOS_TFQMR_ITER_HPP
19 #define BELOS_TFQMR_ITER_HPP
20 
28 #include "BelosConfigDefs.hpp"
29 #include "BelosIteration.hpp"
30 #include "BelosTypes.hpp"
31 
32 #include "BelosLinearProblem.hpp"
33 #include "BelosOutputManager.hpp"
34 #include "BelosStatusTest.hpp"
35 #include "BelosOperatorTraits.hpp"
36 #include "BelosMultiVecTraits.hpp"
37 
40 #include "Teuchos_ScalarTraits.hpp"
42 #include "Teuchos_TimeMonitor.hpp"
43 
55 namespace Belos {
56 
61  template <class ScalarType, class MV>
62  struct TFQMRIterState {
63 
71 
72  TFQMRIterState() : R(Teuchos::null), W(Teuchos::null), U(Teuchos::null),
73  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
74  {}
75  };
76 
77 
79 
80 
87  class TFQMRIterateFailure : public BelosError {public:
88  TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
89  {}};
90 
92 
93 
94  template <class ScalarType, class MV, class OP>
95  class TFQMRIter : public Iteration<ScalarType,MV,OP> {
96  public:
97  //
98  // Convenience typedefs
99  //
104 
106 
107 
110  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
112  Teuchos::ParameterList &params );
113 
115  virtual ~TFQMRIter() {};
117 
118 
120 
121 
132  void iterate();
133 
155  void initializeTFQMR(const TFQMRIterState<ScalarType,MV> & newstate);
156 
160  void initialize()
161  {
163  initializeTFQMR(empty);
164  }
165 
175  state.R = R_;
176  state.W = W_;
177  state.U = U_;
178  state.Rtilde = Rtilde_;
179  state.D = D_;
180  state.V = V_;
181  state.solnUpdate = solnUpdate_;
182  return state;
183  }
184 
186 
187 
189 
190 
192  int getNumIters() const { return iter_; }
193 
195  void resetNumIters( int iter = 0 ) { iter_ = iter; }
196 
199  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
200 
202 
205  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
206 
208 
209 
211 
212 
214  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
215 
217  int getBlockSize() const { return 1; }
218 
220  void setBlockSize(int blockSize) {
221  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
222  "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
223  }
224 
226  bool isInitialized() { return initialized_; }
227 
229 
230 
231  private:
232 
233  //
234  // Internal methods
235  //
237  void setStateSize();
238 
239  //
240  // Classes inputed through constructor that define the linear problem to be solved.
241  //
245 
246  //
247  // Algorithmic parameters
248  //
249 
250  // Storage for QR factorization of the least squares system.
251  // Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
252  std::vector<ScalarType> alpha_, rho_, rho_old_;
253  std::vector<MagnitudeType> tau_, cs_, theta_;
254 
255  //
256  // Current solver state
257  //
258  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
259  // is capable of running; _initialize is controlled by the initialize() member method
260  // For the implications of the state of initialized_, please see documentation for initialize()
261  bool initialized_;
262 
263  // stateStorageInitialized_ specifies that the state storage has be initialized to the current
264  // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
265  // generated without the right-hand side or solution vectors.
266  bool stateStorageInitialized_;
267 
268  // Current subspace dimension, and number of iterations performed.
269  int iter_;
270 
271  //
272  // State Storage
273  //
274  Teuchos::RCP<MV> R_;
275  Teuchos::RCP<MV> W_;
276  Teuchos::RCP<MV> U_, AU_;
277  Teuchos::RCP<MV> Rtilde_;
278  Teuchos::RCP<MV> D_;
279  Teuchos::RCP<MV> V_;
280  Teuchos::RCP<MV> solnUpdate_;
281  };
282 
283 
284  //
285  // Implementation
286  //
287 
289  // Constructor.
290  template <class ScalarType, class MV, class OP>
292  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
294  Teuchos::ParameterList &/* params */
295  ) :
296  lp_(problem),
297  om_(printer),
298  stest_(tester),
299  alpha_(1),
300  rho_(1),
301  rho_old_(1),
302  tau_(1),
303  cs_(1),
304  theta_(1),
305  initialized_(false),
306  stateStorageInitialized_(false),
307  iter_(0)
308  {
309  }
310 
312  // Compute native residual from TFQMR recurrence.
313  template <class ScalarType, class MV, class OP>
315  TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
316  {
318  if (normvec)
319  (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
320 
321  return Teuchos::null;
322  }
323 
324 
326  // Setup the state storage.
327  template <class ScalarType, class MV, class OP>
329  {
330  if (!stateStorageInitialized_) {
331 
332  // Check if there is any multivector to clone from.
333  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
334  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
335  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
336  stateStorageInitialized_ = false;
337  return;
338  }
339  else {
340 
341  // Initialize the state storage
342  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
343  if (R_ == Teuchos::null) {
344  // Get the multivector that is not null.
345  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
346  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
347  "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
348  R_ = MVT::Clone( *tmp, 1 );
349  D_ = MVT::Clone( *tmp, 1 );
350  V_ = MVT::Clone( *tmp, 1 );
351  solnUpdate_ = MVT::Clone( *tmp, 1 );
352  }
353 
354  // State storage has now been initialized.
355  stateStorageInitialized_ = true;
356  }
357  }
358  }
359 
361  // Initialize this iteration object
362  template <class ScalarType, class MV, class OP>
364  {
365  // Initialize the state storage if it isn't already.
366  if (!stateStorageInitialized_)
367  setStateSize();
368 
369  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
370  "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
371 
372  // NOTE: In TFQMRIter R_, the initial residual, is required!!!
373  //
374  std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
375 
376  // Create convenience variables for zero and one.
378 
379  if (newstate.R != Teuchos::null) {
380 
381  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
382  std::invalid_argument, errstr );
383  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
384  std::invalid_argument, errstr );
385 
386  // Copy basis vectors from newstate into V
387  if (newstate.R != R_) {
388  // copy over the initial residual (unpreconditioned).
389  MVT::Assign( *newstate.R, *R_ );
390  }
391 
392  // Compute initial vectors
393  // Initially, they are set to the preconditioned residuals
394  //
395  W_ = MVT::CloneCopy( *R_ );
396  U_ = MVT::CloneCopy( *R_ );
397  Rtilde_ = MVT::CloneCopy( *R_ );
398  MVT::MvInit( *D_ );
399  MVT::MvInit( *solnUpdate_ );
400  // Multiply the current residual by Op and store in V_
401  // V_ = Op * R_
402  //
403  lp_->apply( *U_, *V_ );
404  AU_ = MVT::CloneCopy( *V_ );
405  //
406  // Compute initial scalars: theta, eta, tau, rho_old
407  //
408  theta_[0] = MTzero;
409  MVT::MvNorm( *R_, tau_ ); // tau = ||r_0||
410  MVT::MvDot( *R_, *Rtilde_, rho_old_ ); // rho = (r_tilde, r0)
411  }
412  else {
413 
414  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
415  "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
416  }
417 
418  // The solver is initialized
419  initialized_ = true;
420  }
421 
422 
424  // Iterate until the status test informs us we should stop.
425  template <class ScalarType, class MV, class OP>
427  {
428  //
429  // Allocate/initialize data structures
430  //
431  if (initialized_ == false) {
432  initialize();
433  }
434 
435  // Create convenience variables for zero and one.
436  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
439  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
440  ScalarType eta = STzero, beta = STzero;
441  //
442  // Start executable statements.
443  //
444  // Get the current solution vector.
445  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
446 
447  // Check that the current solution vector only has one column.
448  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
449  "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
450 
451 
453  // Iterate until the status test tells us to stop.
454  //
455  while (stest_->checkStatus(this) != Passed) {
456 
457  for (int iIter=0; iIter<2; iIter++)
458  {
459  //
460  //--------------------------------------------------------
461  // Compute the new alpha if we need to
462  //--------------------------------------------------------
463  //
464  if (iIter == 0) {
465  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
466  alpha_[0] = rho_old_[0]/alpha_[0];
467  }
468  //
469  //--------------------------------------------------------
470  // Update w.
471  // w = w - alpha*Au
472  //--------------------------------------------------------
473  //
474  MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
475  //
476  //--------------------------------------------------------
477  // Update d.
478  // d = u + (theta^2/alpha)eta*d
479  //--------------------------------------------------------
480  //
481  MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
482  //
483  //--------------------------------------------------------
484  // Update u if we need to.
485  // u = u - alpha*v
486  //
487  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
488  //--------------------------------------------------------
489  //
490  if (iIter == 0) {
491  // Compute new U.
492  MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
493 
494  // Update Au for the next iteration.
495  lp_->apply( *U_, *AU_ );
496  }
497  //
498  //--------------------------------------------------------
499  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
500  //--------------------------------------------------------
501  //
502  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
503  theta_[0] /= tau_[0];
504  // cs = 1.0 / sqrt(1.0 + theta^2)
505  cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
506  tau_[0] *= theta_[0]*cs_[0]; // tau = tau * theta * cs
507  eta = cs_[0]*cs_[0]*alpha_[0]; // eta = cs^2 * alpha
508  //
509  //--------------------------------------------------------
510  // Update the solution.
511  // Don't update the linear problem object, may incur additional preconditioner application.
512  //--------------------------------------------------------
513  //
514  MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
515  //
516  //--------------------------------------------------------
517  // Check for breakdown before continuing.
518  //--------------------------------------------------------
519  if ( tau_[0] == MTzero ) {
520  break;
521  }
522  //
523  if (iIter == 1) {
524  //
525  //--------------------------------------------------------
526  // Compute the new rho, beta if we need to.
527  //--------------------------------------------------------
528  //
529  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
530  beta = rho_[0]/rho_old_[0]; // beta = rho / rho_old
531  rho_old_[0] = rho_[0]; // rho_old = rho
532  //
533  //--------------------------------------------------------
534  // Update u, v, and Au if we need to.
535  // Note: We are updating v in two stages to be memory efficient
536  //--------------------------------------------------------
537  //
538  MVT::MvAddMv( STone, *W_, beta, *U_, *U_ ); // u = w + beta*u
539 
540  // First stage of v update.
541  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
542 
543  // Update Au.
544  lp_->apply( *U_, *AU_ ); // Au = A*u
545 
546  // Second stage of v update.
547  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
548  }
549 
550  }
551 
552  // Increment the iteration
553  iter_++;
554 
555  } // end while (sTest_->checkStatus(this) != Passed)
556  }
557 
558 } // namespace Belos
559 //
560 #endif // BELOS_TFQMR_ITER_HPP
561 //
562 // End of file BelosTFQMRIter.hpp
563 
564 
Teuchos::RCP< const MV > V
Collection of types and exceptions used within the Belos solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::ScalarTraits< ScalarType > SCT
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
TFQMRIterateFailure(const std::string &what_arg)
void setBlockSize(int blockSize)
Set the blocksize.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
TFQMRIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList &params)
Belos::TFQMRIter constructor.
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< const MV > U
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Rtilde
void iterate()
This method performs TFQMR iterations until the status test indicates the need to stop or an error oc...
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< const MV > W
A linear system to solve, and its associated information.
Structure to contain pointers to TFQMRIter state variables.
Class which describes the linear problem to be solved by the iterative solver.
int getNumIters() const
Get the current iteration count.
MultiVecTraits< ScalarType, MV > MVT
TFQMRIterateFailure is thrown when the TFQMRIter object is unable to compute the next iterate in the ...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Teuchos::RCP< const MV > D
Belos header file which uses auto-configuration information to include necessary C++ headers...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.

Generated on Mon Jul 15 2024 09:24:24 for Belos by doxygen 1.8.5