Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockTFQMRIter.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_PSEUDO_BLOCK_TFQMR_ITER_HPP
19 #define BELOS_PSEUDO_BLOCK_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 
38 #include "Teuchos_ScalarTraits.hpp"
40 #include "Teuchos_TimeMonitor.hpp"
41 
53 namespace Belos {
54 
59  template <class ScalarType, class MV>
61 
64 
72  std::vector<ScalarType> alpha, eta, rho;
73  std::vector<MagnitudeType> tau, theta;
74 
75 
76  PseudoBlockTFQMRIterState() : W(Teuchos::null), U(Teuchos::null), AU(Teuchos::null),
77  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
78  {}
79  };
80 
81  template <class ScalarType, class MV, class OP>
82  class PseudoBlockTFQMRIter : public Iteration<ScalarType,MV,OP> {
83  public:
84  //
85  // Convenience typedefs
86  //
91 
93 
94 
97  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
99  Teuchos::ParameterList &params );
100 
102  virtual ~PseudoBlockTFQMRIter() {};
104 
105 
107 
108 
119  void iterate();
120 
143 
147  void initialize()
148  {
150  initializeTFQMR(empty);
151  }
152 
162 
163  // Copy over the vectors.
164  state.W = W_;
165  state.U = U_;
166  state.AU = AU_;
167  state.Rtilde = Rtilde_;
168  state.D = D_;
169  state.V = V_;
170 
171  // Copy over the scalars.
172  state.alpha = alpha_;
173  state.eta = eta_;
174  state.rho = rho_;
175  state.tau = tau_;
176  state.theta = theta_;
177 
178  return state;
179  }
180 
182 
183 
185 
186 
188  int getNumIters() const { return iter_; }
189 
191  void resetNumIters( int iter = 0 ) { iter_ = iter; }
192 
195  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
196 
198 
201  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
202 
204 
205 
207 
208 
210  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
211 
213  int getBlockSize() const { return 1; }
214 
216  void setBlockSize(int blockSize) {
217  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
218  "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
219  }
220 
222  bool isInitialized() { return initialized_; }
223 
225 
226 
227  private:
228 
229  //
230  // Classes inputed through constructor that define the linear problem to be solved.
231  //
235 
236  //
237  // Algorithmic parameters
238  //
239 
240  // numRHS_ is the current number of linear systems being solved.
241  int numRHS_;
242 
243  // Storage for QR factorization of the least squares system.
244  std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
245  std::vector<MagnitudeType> tau_, theta_;
246 
247  //
248  // Current solver state
249  //
250  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
251  // is capable of running; _initialize is controlled by the initialize() member method
252  // For the implications of the state of initialized_, please see documentation for initialize()
253  bool initialized_;
254 
255  // Current subspace dimension, and number of iterations performed.
256  int iter_;
257 
258  //
259  // State Storage
260  //
261  Teuchos::RCP<MV> W_;
262  Teuchos::RCP<MV> U_, AU_;
263  Teuchos::RCP<MV> Rtilde_;
264  Teuchos::RCP<MV> D_;
265  Teuchos::RCP<MV> V_;
266  Teuchos::RCP<MV> solnUpdate_;
267 
268  };
269 
270 
271  //
272  // Implementation
273  //
274 
276  // Constructor.
277  template <class ScalarType, class MV, class OP>
279  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
281  Teuchos::ParameterList &/* params */
282  ) :
283  lp_(problem),
284  om_(printer),
285  stest_(tester),
286  numRHS_(0),
287  initialized_(false),
288  iter_(0)
289  {
290  }
291 
293  // Compute native residual from TFQMR recurrence.
294  template <class ScalarType, class MV, class OP>
296  PseudoBlockTFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
297  {
299  if (normvec) {
300  // Resize the vector passed in, if it is too small.
301  if ((int) normvec->size() < numRHS_)
302  normvec->resize( numRHS_ );
303 
304  // Compute the native residuals.
305  for (int i=0; i<numRHS_; i++) {
306  (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
307  }
308  }
309 
310  return Teuchos::null;
311  }
312 
314  // Initialize this iteration object
315  template <class ScalarType, class MV, class OP>
317  {
318  // Create convenience variables for zero and one.
319  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
320  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
322 
323  // NOTE: In PseudoBlockTFQMRIter Rtilde_, the initial residual, is required!!!
324  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Rtilde == Teuchos::null,std::invalid_argument,
325  "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
326 
327  // Get the number of right-hand sides we're solving for now.
328  int numRHS = MVT::GetNumberVecs(*newstate.Rtilde);
329  numRHS_ = numRHS;
330 
331  // Initialize the state storage
332  // If the subspace has not be initialized before or we are reusing this solver object, generate it using Rtilde.
333  if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
334  {
335  // Create and/or initialize D_.
336  if ( Teuchos::is_null(D_) )
337  D_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
338  MVT::MvInit( *D_, STzero );
339 
340  // Create and/or initialize solnUpdate_;
341  if ( Teuchos::is_null(solnUpdate_) )
342  solnUpdate_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
343  MVT::MvInit( *solnUpdate_, STzero );
344 
345  // Create Rtilde_.
346  if (newstate.Rtilde != Rtilde_)
347  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
348  W_ = MVT::CloneCopy( *Rtilde_ );
349  U_ = MVT::CloneCopy( *Rtilde_ );
350  V_ = MVT::Clone( *Rtilde_, numRHS_ );
351 
352  // Multiply the current residual by Op and store in V_
353  // V_ = Op * R_
354  lp_->apply( *U_, *V_ );
355  AU_ = MVT::CloneCopy( *V_ );
356 
357  // Resize work vectors.
358  alpha_.resize( numRHS_, STone );
359  eta_.resize( numRHS_, STzero );
360  rho_.resize( numRHS_ );
361  rho_old_.resize( numRHS_ );
362  tau_.resize( numRHS_ );
363  theta_.resize( numRHS_, MTzero );
364 
365  MVT::MvNorm( *Rtilde_, tau_ ); // tau = ||r_0||
366  MVT::MvDot( *Rtilde_, *Rtilde_, rho_ ); // rho = (r_tilde, r0)
367  }
368  else
369  {
370  // If the subspace has changed sizes, clone it from the incoming state.
371  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
372  W_ = MVT::CloneCopy( *newstate.W );
373  U_ = MVT::CloneCopy( *newstate.U );
374  AU_ = MVT::CloneCopy( *newstate.AU );
375  D_ = MVT::CloneCopy( *newstate.D );
376  V_ = MVT::CloneCopy( *newstate.V );
377 
378  // The solution update is just set to zero, since the current update has already
379  // been added to the solution during deflation.
380  solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
381  MVT::MvInit( *solnUpdate_, STzero );
382 
383  // Copy work vectors.
384  alpha_ = newstate.alpha;
385  eta_ = newstate.eta;
386  rho_ = newstate.rho;
387  tau_ = newstate.tau;
388  theta_ = newstate.theta;
389  }
390 
391  // The solver is initialized
392  initialized_ = true;
393  }
394 
395 
397  // Iterate until the status test informs us we should stop.
398  template <class ScalarType, class MV, class OP>
400  {
401  //
402  // Allocate/initialize data structures
403  //
404  if (initialized_ == false) {
405  initialize();
406  }
407 
408  // Create convenience variables for zero and one.
409  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
410  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
413  std::vector< ScalarType > beta(numRHS_,STzero);
414  std::vector<int> index(1);
415  //
416  // Start executable statements.
417  //
418 
420  // Iterate until the status test tells us to stop.
421  //
422  while (stest_->checkStatus(this) != Passed) {
423 
424  for (int iIter=0; iIter<2; iIter++)
425  {
426  //
427  //--------------------------------------------------------
428  // Compute the new alpha if we need to
429  //--------------------------------------------------------
430  //
431  if (iIter == 0) {
432  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
433  for (int i=0; i<numRHS_; i++) {
434  rho_old_[i] = rho_[i]; // rho_old = rho
435  alpha_[i] = rho_old_[i]/alpha_[i];
436  }
437  }
438  //
439  //--------------------------------------------------------
440  // Loop over all RHS and compute updates.
441  //--------------------------------------------------------
442  //
443  for (int i=0; i<numRHS_; ++i) {
444  index[0] = i;
445 
446  //
447  //--------------------------------------------------------
448  // Update w.
449  // w = w - alpha*Au
450  //--------------------------------------------------------
451  //
452  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
453  Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
454  MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
455  //
456  //--------------------------------------------------------
457  // Update d.
458  // d = u + (theta^2/alpha)eta*d
459  //--------------------------------------------------------
460  //
461  Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
462  Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
463  MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
464  //
465  //--------------------------------------------------------
466  // Update u if we need to.
467  // u = u - alpha*v
468  //
469  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
470  //--------------------------------------------------------
471  //
472  if (iIter == 0) {
473  // Compute new U.
474  Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
475  Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
476  MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
477  }
478  }
479  //
480  //--------------------------------------------------------
481  // Update Au for the next iteration.
482  //--------------------------------------------------------
483  //
484  if (iIter == 0) {
485  lp_->apply( *U_, *AU_ );
486  }
487  //
488  //--------------------------------------------------------
489  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
490  //--------------------------------------------------------
491  //
492  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
493 
494  bool breakdown=false;
495  for (int i=0; i<numRHS_; ++i) {
496  theta_[i] /= tau_[i];
497  // cs = 1.0 / sqrt(1.0 + theta^2)
498  MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
499  tau_[i] *= theta_[i]*cs; // tau = tau * theta * cs
500  if ( tau_[i] == MTzero ) {
501  breakdown = true;
502  }
503  eta_[i] = cs*cs*alpha_[i]; // eta = cs^2 * alpha
504  }
505  //
506  //--------------------------------------------------------
507  // Accumulate the update for the solution x := x + eta*D_
508  //--------------------------------------------------------
509  //
510  for (int i=0; i<numRHS_; ++i) {
511  index[0]=i;
512  Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
513  Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
514  MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
515  }
516  //
517  //--------------------------------------------------------
518  // Breakdown was detected above, return to status test to
519  // remove converged solutions.
520  //--------------------------------------------------------
521  if ( breakdown ) {
522  break;
523  }
524  //
525  if (iIter == 1) {
526  //
527  //--------------------------------------------------------
528  // Compute the new rho, beta if we need to.
529  //--------------------------------------------------------
530  //
531  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
532 
533  for (int i=0; i<numRHS_; ++i) {
534  beta[i] = rho_[i]/rho_old_[i]; // beta = rho / rho_old
535 
536  //
537  //--------------------------------------------------------
538  // Update u, v, and Au if we need to.
539  // Note: We are updating v in two stages to be memory efficient
540  //--------------------------------------------------------
541  //
542  index[0]=i;
543  Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
544  Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
545  MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i ); // u = w + beta*u
546 
547  // First stage of v update.
548  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
549  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
550  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
551  }
552 
553  // Update Au.
554  lp_->apply( *U_, *AU_ ); // Au = A*u
555 
556  // Second stage of v update.
557  for (int i=0; i<numRHS_; ++i) {
558  index[0]=i;
559  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
560  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
561  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
562  }
563  }
564  }
565 
566  // Increment the iteration
567  iter_++;
568 
569  } // end while (sTest_->checkStatus(this) != Passed)
570  }
571 
572 } // namespace Belos
573 //
574 #endif // BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
575 //
576 // End of file BelosPseudoBlockTFQMRIter.hpp
577 
578 
Collection of types and exceptions used within the Belos solvers.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
OperatorTraits< ScalarType, MV, OP > OPT
Class which manages the output and verbosity of the Belos solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
bool is_null(const std::shared_ptr< T > &p)
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getNumIters() const
Get the current iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Traits class which defines basic operations on multivectors.
PseudoBlockTFQMRIter(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::PseudoBlockTFQMRIter constructor.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setBlockSize(int blockSize)
Set the blocksize.
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
Teuchos::RCP< const MV > W
The current residual basis.
Class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType > SCT
void iterate()
This method performs pseudo-block TFQMR iterations until the status test indicates the need to stop o...
void resetNumIters(int iter=0)
Reset the iteration count.
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool isInitialized()
States whether the solver has been initialized or not.
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.

Generated on Fri Dec 20 2024 09:24:49 for Belos by doxygen 1.8.5