Belos  Version of the Day
BelosTFQMRIter.hpp
Go to the documentation of this file.
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
7 // *****************************************************************************
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  1.8.5