Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBiCGStabIter.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 #ifndef BELOS_BICGSTAB_ITER_HPP
11 #define BELOS_BICGSTAB_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 #include "BelosCGIteration.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosMatOrthoManager.hpp"
23 #include "BelosOutputManager.hpp"
24 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
30 #include "Teuchos_ScalarTraits.hpp"
32 #include "Teuchos_TimeMonitor.hpp"
33 
45 namespace Belos {
46 
48 
49 
54  template <class ScalarType, class MV>
56 
59 
62 
65 
68 
69  std::vector<ScalarType> rho_old, alpha, omega;
70 
71  BiCGStabIterationState() : R(Teuchos::null), Rhat(Teuchos::null),
72  P(Teuchos::null), V(Teuchos::null)
73  {
74  rho_old.clear();
75  alpha.clear();
76  omega.clear();
77  }
78  };
79 
80  template<class ScalarType, class MV, class OP>
81  class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP> {
82 
83  public:
84 
85  //
86  // Convenience typedefs
87  //
93 
95 
96 
103  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
105  Teuchos::ParameterList &params );
106 
108  virtual ~BiCGStabIter() {};
110 
111 
113 
114 
128  void iterate();
129 
151 
155  void initialize()
156  {
158  initializeBiCGStab(empty);
159  }
160 
170  state.R = R_;
171  state.Rhat = Rhat_;
172  state.P = P_;
173  state.V = V_;
174  state.rho_old = rho_old_;
175  state.alpha = alpha_;
176  state.omega = omega_;
177  return state;
178  }
179 
181 
182 
184 
185 
187  int getNumIters() const { return iter_; }
188 
190  void resetNumIters( int iter = 0 ) { iter_ = iter; }
191 
194  // amk TODO: are the residuals actually being set? What is a native residual?
195  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
196 
198 
200  // amk TODO: what is this supposed to be doing?
201  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
202 
204  bool breakdownDetected() { return breakdown_; }
205 
207 
209 
210 
212  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
213 
215  int getBlockSize() const { return 1; }
216 
218  void setBlockSize(int blockSize) {
219  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
220  "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
221  }
222 
224  bool isInitialized() { return initialized_; }
225 
227 
228  private:
229 
230  void axpy(const ScalarType alpha, const MV & A,
231  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus=false);
232 
233  //
234  // Classes inputed through constructor that define the linear problem to be solved.
235  //
239 
240  //
241  // Algorithmic parameters
242  //
243  // numRHS_ is the current number of linear systems being solved.
244  int numRHS_;
245 
246  //
247  // Current solver state
248  //
249  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
250  // is capable of running; _initialize is controlled by the initialize() member method
251  // For the implications of the state of initialized_, please see documentation for initialize()
252  bool initialized_;
253 
254  // Breakdown has been observed for at least one of the linear systems
255  bool breakdown_;
256 
257  // Current number of iterations performed.
258  int iter_;
259 
260  //
261  // State Storage
262  //
263  // Initial residual
264  Teuchos::RCP<MV> Rhat_;
265  //
266  // Residual
267  Teuchos::RCP<MV> R_;
268  //
269  // Direction vector 1
270  Teuchos::RCP<MV> P_;
271  //
272  // Operator applied to preconditioned direction vector 1
273  Teuchos::RCP<MV> V_;
274  //
275  std::vector<ScalarType> rho_old_, alpha_, omega_;
276  };
277 
279  // Constructor.
280  template<class ScalarType, class MV, class OP>
282  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
284  Teuchos::ParameterList &/* params */ ):
285  lp_(problem),
286  om_(printer),
287  stest_(tester),
288  numRHS_(0),
289  initialized_(false),
290  breakdown_(false),
291  iter_(0)
292  {
293  }
294 
295 
297  // Initialize this iteration object
298  template <class ScalarType, class MV, class OP>
300  {
301  // Check if there is any multivector to clone from.
302  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
303  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
304  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
305  "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
306 
307  // Get the multivector that is not null.
308  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
309 
310  // Get the number of right-hand sides we're solving for now.
311  int numRHS = MVT::GetNumberVecs(*tmp);
312  numRHS_ = numRHS;
313 
314  // Initialize the state storage
315  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
316  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
317  R_ = MVT::Clone( *tmp, numRHS_ );
318  Rhat_ = MVT::Clone( *tmp, numRHS_ );
319  P_ = MVT::Clone( *tmp, numRHS_ );
320  V_ = MVT::Clone( *tmp, numRHS_ );
321 
322  rho_old_.resize(numRHS_);
323  alpha_.resize(numRHS_);
324  omega_.resize(numRHS_);
325  }
326 
327  // Reset breakdown to false before initializing iteration
328  breakdown_ = false;
329 
330  // NOTE: In BiCGStabIter R_, the initial residual, is required!!!
331  //
332  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
333 
334  // Create convenience variable for one.
335  const ScalarType one = SCT::one();
336 
337  if (!Teuchos::is_null(newstate.R)) {
338 
339  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
340  std::invalid_argument, errstr );
341  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
342  std::invalid_argument, errstr );
343 
344  // Copy residual vectors from newstate into R
345  if (newstate.R != R_) {
346  // Assigned by the new state
347  MVT::Assign(*newstate.R, *R_);
348  }
349  else {
350  // Computed
351  lp_->computeCurrResVec(R_.get());
352  }
353 
354  // Set Rhat
355  if (!Teuchos::is_null(newstate.Rhat) && newstate.Rhat != Rhat_) {
356  // Assigned by the new state
357  MVT::Assign(*newstate.Rhat, *Rhat_);
358  }
359  else {
360  // Set to be the initial residual
361  MVT::Assign(*R_, *Rhat_);
362  }
363 
364  // Set V
365  if (!Teuchos::is_null(newstate.V) && newstate.V != V_) {
366  // Assigned by the new state
367  MVT::Assign(*newstate.V, *V_);
368  }
369  else {
370  // Initial V = 0
371  MVT::MvInit(*V_);
372  }
373 
374  // Set P
375  if (!Teuchos::is_null(newstate.P) && newstate.P != P_) {
376  // Assigned by the new state
377  MVT::Assign(*newstate.P, *P_);
378  }
379  else {
380  // Initial P = 0
381  MVT::MvInit(*P_);
382  }
383 
384  // Set rho_old
385  if (newstate.rho_old.size () == static_cast<size_t> (numRHS_)) {
386  // Assigned by the new state
387  rho_old_ = newstate.rho_old;
388  }
389  else {
390  // Initial rho = 1
391  rho_old_.assign(numRHS_,one);
392  }
393 
394  // Set alpha
395  if (newstate.alpha.size() == static_cast<size_t> (numRHS_)) {
396  // Assigned by the new state
397  alpha_ = newstate.alpha;
398  }
399  else {
400  // Initial rho = 1
401  alpha_.assign(numRHS_,one);
402  }
403 
404  // Set omega
405  if (newstate.omega.size() == static_cast<size_t> (numRHS_)) {
406  // Assigned by the new state
407  omega_ = newstate.omega;
408  }
409  else {
410  // Initial rho = 1
411  omega_.assign(numRHS_,one);
412  }
413 
414  }
415  else {
416 
417  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
418  "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
419  }
420 
421  // The solver is initialized
422  initialized_ = true;
423  }
424 
425 
427  // Iterate until the status test informs us we should stop.
428  template <class ScalarType, class MV, class OP>
430  {
431  using Teuchos::RCP;
432 
433  //
434  // Allocate/initialize data structures
435  //
436  if (initialized_ == false) {
437  initialize();
438  }
439 
440  // Allocate memory for scalars.
441  int i=0;
442  std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
443  std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
444 
445  // Create convenience variable for one.
446  const ScalarType one = SCT::one();
447 
448  // TODO: We may currently be using more space than is required
449  RCP<MV> leftPrecVec, leftPrecVec2;
450 
451  RCP<MV> Y, Z, S, T;
452  S = MVT::Clone( *R_, numRHS_ );
453  T = MVT::Clone( *R_, numRHS_ );
454  if (lp_->isLeftPrec() || lp_->isRightPrec()) {
455  Y = MVT::Clone( *R_, numRHS_ );
456  Z = MVT::Clone( *R_, numRHS_ );
457  }
458  else {
459  Y = P_;
460  Z = S;
461  }
462 
463  // Get the current solution std::vector.
464  Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
465 
467  // Iterate until the status test tells us to stop.
468  //
469  while (stest_->checkStatus(this) != Passed && !breakdown_) {
470 
471  // Increment the iteration
472  iter_++;
473 
474  // rho_new = <R_, Rhat_>
475  MVT::MvDot(*R_,*Rhat_,rho_new);
476 
477  // beta = ( rho_new / rho_old ) (alpha / omega )
478  // TODO: None of these loops are currently threaded
479  for(i=0; i<numRHS_; i++) {
480  // Catch breakdown in rho_old here, since
481  // it is just rho_new from the previous iteration.
482  if (SCT::magnitude(rho_new[i]) < MT::sfmin())
483  breakdown_ = true;
484 
485  beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
486  }
487 
488  // p = r + beta (p - omega v)
489  // TODO: Is it safe to call MvAddMv with A or B = mv?
490  // TODO: Not all of these things have to be part of the state
491  axpy(one, *P_, omega_, *V_, *P_, true); // p = p - omega v
492  axpy(one, *R_, beta, *P_, *P_); // p = r + beta (p - omega v)
493 
494  // y = K\p, unless K does not exist
495  // TODO: There may be a more efficient way to apply the preconditioners
496  if(lp_->isLeftPrec()) {
497  if(lp_->isRightPrec()) {
498  if(leftPrecVec == Teuchos::null) {
499  leftPrecVec = MVT::Clone( *R_, numRHS_ );
500  }
501  lp_->applyLeftPrec(*P_,*leftPrecVec);
502  lp_->applyRightPrec(*leftPrecVec,*Y);
503  }
504  else {
505  lp_->applyLeftPrec(*P_,*Y);
506  }
507  }
508  else if(lp_->isRightPrec()) {
509  lp_->applyRightPrec(*P_,*Y);
510  }
511 
512  // v = Ay
513  lp_->applyOp(*Y,*V_);
514 
515  // alpha = rho_new / <Rhat, V>
516  MVT::MvDot(*V_,*Rhat_,rhatV);
517  for(i=0; i<numRHS_; i++) {
518  if (SCT::magnitude(rhatV[i]) < MT::sfmin())
519  {
520  breakdown_ = true;
521  return;
522  }
523  else
524  alpha_[i] = rho_new[i] / rhatV[i];
525  }
526 
527  // s = r - alpha v
528  axpy(one, *R_, alpha_, *V_, *S, true);
529 
530  // z = K\s, unless K does not exist
531  if(lp_->isLeftPrec()) {
532  if(lp_->isRightPrec()) {
533  if(leftPrecVec == Teuchos::null) {
534  leftPrecVec = MVT::Clone( *R_, numRHS_ );
535  }
536  lp_->applyLeftPrec(*S,*leftPrecVec);
537  lp_->applyRightPrec(*leftPrecVec,*Z);
538  }
539  else {
540  lp_->applyLeftPrec(*S,*Z);
541  }
542  }
543  else if(lp_->isRightPrec()) {
544  lp_->applyRightPrec(*S,*Z);
545  }
546 
547  // t = Az
548  lp_->applyOp(*Z,*T);
549 
550  // omega = <K1\t,K1\s> / <K1\t,K1\t>
551  if(lp_->isLeftPrec()) {
552  if(leftPrecVec == Teuchos::null) {
553  leftPrecVec = MVT::Clone( *R_, numRHS_ );
554  }
555  if(leftPrecVec2 == Teuchos::null) {
556  leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
557  }
558  lp_->applyLeftPrec(*T,*leftPrecVec2);
559  MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
560  MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
561  }
562  else {
563  MVT::MvDot(*T,*T,tT);
564  MVT::MvDot(*S,*T,tS);
565  }
566  for(i=0; i<numRHS_; i++) {
567  if (SCT::magnitude(tT[i]) < MT::sfmin())
568  {
569  omega_[i] = SCT::zero();
570  breakdown_ = true;
571  }
572  else
573  omega_[i] = tS[i] / tT[i];
574  }
575 
576  // x = x + alpha y + omega z
577  axpy(one, *X, alpha_, *Y, *X); // x = x + alpha y
578  axpy(one, *X, omega_, *Z, *X); // x = x + alpha y + omega z
579 
580  // r = s - omega t
581  axpy(one, *S, omega_, *T, *R_, true);
582 
583  // Update rho_old
584  rho_old_ = rho_new;
585  } // end while (sTest_->checkStatus(this) != Passed)
586  }
587 
588 
590  // Iterate until the status test informs us we should stop.
591  template <class ScalarType, class MV, class OP>
592  void BiCGStabIter<ScalarType,MV,OP>::axpy(const ScalarType alpha, const MV & A,
593  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus)
594  {
595  Teuchos::RCP<const MV> A1, B1;
596  Teuchos::RCP<MV> mv1;
597  std::vector<int> index(1);
598 
599  for(int i=0; i<numRHS_; i++) {
600  index[0] = i;
601  A1 = MVT::CloneView(A,index);
602  B1 = MVT::CloneView(B,index);
603  mv1 = MVT::CloneViewNonConst(mv,index);
604  if(minus) {
605  MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
606  }
607  else {
608  MVT::MvAddMv(alpha,*A1,beta[i],*B1,*mv1);
609  }
610  }
611  }
612 
613 } // end Belos namespace
614 
615 #endif /* BELOS_BICGSTAB_ITER_HPP */
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...
This class implements the pseudo-block BiCGStab iteration, where the basic BiCGStab algorithm is perf...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > R
The current residual.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > Rhat
The initial residual.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Structure to contain pointers to BiCGStabIteration state variables.
Declaration of basic traits for the multivector type.
void initializeBiCGStab(BiCGStabIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void setBlockSize(int blockSize)
Set the blocksize.
void resetNumIters(int iter=0)
Reset the iteration count.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
std::vector< ScalarType > alpha
Traits class which defines basic operations on multivectors.
std::vector< ScalarType > omega
BiCGStabIter(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)
BiCGStabIter constructor with linear problem, solver utilities, and parameter list of solver options...
A linear system to solve, and its associated information.
Teuchos::RCP< const MV > V
A * M * the first decent direction vector.
Class which describes the linear problem to be solved by the iterative solver.
virtual ~BiCGStabIter()
Destructor.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< const MV > P
The first decent direction vector.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void iterate()
This method performs BiCGStab iterations on each linear system until the status test indicates the ne...
std::vector< ScalarType > rho_old
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
SCT::magnitudeType MagnitudeType
MultiVecTraits< ScalarType, MV > MVT
bool breakdownDetected()
Has breakdown been detected in any linear system.
Belos header file which uses auto-configuration information to include necessary C++ headers...
BiCGStabIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated on Wed Dec 11 2024 09:24:06 for Belos by doxygen 1.8.5