Belos  Version of the Day
BelosBiCGStabIter.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_BICGSTAB_ITER_HPP
43 #define BELOS_BICGSTAB_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosCGIteration.hpp"
52 
53 #include "BelosLinearProblem.hpp"
54 #include "BelosMatOrthoManager.hpp"
55 #include "BelosOutputManager.hpp"
56 #include "BelosStatusTest.hpp"
57 #include "BelosOperatorTraits.hpp"
58 #include "BelosMultiVecTraits.hpp"
59 
60 #include "Teuchos_SerialDenseMatrix.hpp"
61 #include "Teuchos_SerialDenseVector.hpp"
62 #include "Teuchos_ScalarTraits.hpp"
63 #include "Teuchos_ParameterList.hpp"
64 #include "Teuchos_TimeMonitor.hpp"
65 
77 namespace Belos {
78 
80 
81 
86  template <class ScalarType, class MV>
88 
90  Teuchos::RCP<const MV> R;
91 
93  Teuchos::RCP<const MV> Rhat;
94 
96  Teuchos::RCP<const MV> P;
97 
99  Teuchos::RCP<const MV> V;
100 
101  std::vector<ScalarType> rho_old, alpha, omega;
102 
103  BiCGStabIterationState() : R(Teuchos::null), Rhat(Teuchos::null),
104  P(Teuchos::null), V(Teuchos::null)
105  {
106  rho_old.clear();
107  alpha.clear();
108  omega.clear();
109  }
110  };
111 
112  template<class ScalarType, class MV, class OP>
113  class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP> {
114 
115  public:
116 
117  //
118  // Convenience typedefs
119  //
122  typedef Teuchos::ScalarTraits<ScalarType> SCT;
123  typedef typename SCT::magnitudeType MagnitudeType;
124  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
125 
127 
128 
134  BiCGStabIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
135  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
136  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
137  Teuchos::ParameterList &params );
138 
140  virtual ~BiCGStabIter() {};
142 
143 
145 
146 
160  void iterate();
161 
183 
187  void initialize()
188  {
190  initializeBiCGStab(empty);
191  }
192 
202  state.R = R_;
203  state.Rhat = Rhat_;
204  state.P = P_;
205  state.V = V_;
206  state.rho_old = rho_old_;
207  state.alpha = alpha_;
208  state.omega = omega_;
209  return state;
210  }
211 
213 
214 
216 
217 
219  int getNumIters() const { return iter_; }
220 
222  void resetNumIters( int iter = 0 ) { iter_ = iter; }
223 
226  // amk TODO: are the residuals actually being set? What is a native residual?
227  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
228 
230 
232  // amk TODO: what is this supposed to be doing?
233  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
234 
236  bool breakdownDetected() { return breakdown_; }
237 
239 
241 
242 
244  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
245 
247  int getBlockSize() const { return 1; }
248 
250  void setBlockSize(int blockSize) {
251  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
252  "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
253  }
254 
256  bool isInitialized() { return initialized_; }
257 
259 
260  private:
261 
262  void axpy(const ScalarType alpha, const MV & A,
263  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus=false);
264 
265  //
266  // Classes inputed through constructor that define the linear problem to be solved.
267  //
268  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
269  const Teuchos::RCP<OutputManager<ScalarType> > om_;
270  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
271 
272  //
273  // Algorithmic parameters
274  //
275  // numRHS_ is the current number of linear systems being solved.
276  int numRHS_;
277 
278  //
279  // Current solver state
280  //
281  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
282  // is capable of running; _initialize is controlled by the initialize() member method
283  // For the implications of the state of initialized_, please see documentation for initialize()
284  bool initialized_;
285 
286  // Breakdown has been observed for at least one of the linear systems
287  bool breakdown_;
288 
289  // Current number of iterations performed.
290  int iter_;
291 
292  //
293  // State Storage
294  //
295  // Initial residual
296  Teuchos::RCP<MV> Rhat_;
297  //
298  // Residual
299  Teuchos::RCP<MV> R_;
300  //
301  // Direction vector 1
302  Teuchos::RCP<MV> P_;
303  //
304  // Operator applied to preconditioned direction vector 1
305  Teuchos::RCP<MV> V_;
306  //
307  std::vector<ScalarType> rho_old_, alpha_, omega_;
308  };
309 
311  // Constructor.
312  template<class ScalarType, class MV, class OP>
314  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
315  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
316  Teuchos::ParameterList &/* params */ ):
317  lp_(problem),
318  om_(printer),
319  stest_(tester),
320  numRHS_(0),
321  initialized_(false),
322  breakdown_(false),
323  iter_(0)
324  {
325  }
326 
327 
329  // Initialize this iteration object
330  template <class ScalarType, class MV, class OP>
332  {
333  // Check if there is any multivector to clone from.
334  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
335  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
336  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
337  "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
338 
339  // Get the multivector that is not null.
340  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
341 
342  // Get the number of right-hand sides we're solving for now.
343  int numRHS = MVT::GetNumberVecs(*tmp);
344  numRHS_ = numRHS;
345 
346  // Initialize the state storage
347  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
348  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
349  R_ = MVT::Clone( *tmp, numRHS_ );
350  Rhat_ = MVT::Clone( *tmp, numRHS_ );
351  P_ = MVT::Clone( *tmp, numRHS_ );
352  V_ = MVT::Clone( *tmp, numRHS_ );
353 
354  rho_old_.resize(numRHS_);
355  alpha_.resize(numRHS_);
356  omega_.resize(numRHS_);
357  }
358 
359  // Reset breakdown to false before initializing iteration
360  breakdown_ = false;
361 
362  // NOTE: In BiCGStabIter R_, the initial residual, is required!!!
363  //
364  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
365 
366  // Create convenience variable for one.
367  const ScalarType one = SCT::one();
368 
369  if (!Teuchos::is_null(newstate.R)) {
370 
371  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
372  std::invalid_argument, errstr );
373  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
374  std::invalid_argument, errstr );
375 
376  // Copy residual vectors from newstate into R
377  if (newstate.R != R_) {
378  // Assigned by the new state
379  MVT::Assign(*newstate.R, *R_);
380  }
381  else {
382  // Computed
383  lp_->computeCurrResVec(R_.get());
384  }
385 
386  // Set Rhat
387  if (!Teuchos::is_null(newstate.Rhat) && newstate.Rhat != Rhat_) {
388  // Assigned by the new state
389  MVT::Assign(*newstate.Rhat, *Rhat_);
390  }
391  else {
392  // Set to be the initial residual
393  MVT::Assign(*R_, *Rhat_);
394  }
395 
396  // Set V
397  if (!Teuchos::is_null(newstate.V) && newstate.V != V_) {
398  // Assigned by the new state
399  MVT::Assign(*newstate.V, *V_);
400  }
401  else {
402  // Initial V = 0
403  MVT::MvInit(*V_);
404  }
405 
406  // Set P
407  if (!Teuchos::is_null(newstate.P) && newstate.P != P_) {
408  // Assigned by the new state
409  MVT::Assign(*newstate.P, *P_);
410  }
411  else {
412  // Initial P = 0
413  MVT::MvInit(*P_);
414  }
415 
416  // Set rho_old
417  if (newstate.rho_old.size () == static_cast<size_t> (numRHS_)) {
418  // Assigned by the new state
419  rho_old_ = newstate.rho_old;
420  }
421  else {
422  // Initial rho = 1
423  rho_old_.assign(numRHS_,one);
424  }
425 
426  // Set alpha
427  if (newstate.alpha.size() == static_cast<size_t> (numRHS_)) {
428  // Assigned by the new state
429  alpha_ = newstate.alpha;
430  }
431  else {
432  // Initial rho = 1
433  alpha_.assign(numRHS_,one);
434  }
435 
436  // Set omega
437  if (newstate.omega.size() == static_cast<size_t> (numRHS_)) {
438  // Assigned by the new state
439  omega_ = newstate.omega;
440  }
441  else {
442  // Initial rho = 1
443  omega_.assign(numRHS_,one);
444  }
445 
446  }
447  else {
448 
449  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
450  "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
451  }
452 
453  // The solver is initialized
454  initialized_ = true;
455  }
456 
457 
459  // Iterate until the status test informs us we should stop.
460  template <class ScalarType, class MV, class OP>
462  {
463  using Teuchos::RCP;
464 
465  //
466  // Allocate/initialize data structures
467  //
468  if (initialized_ == false) {
469  initialize();
470  }
471 
472  // Allocate memory for scalars.
473  int i=0;
474  std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
475  std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
476 
477  // Create convenience variable for one.
478  const ScalarType one = SCT::one();
479 
480  // TODO: We may currently be using more space than is required
481  RCP<MV> leftPrecVec, leftPrecVec2;
482 
483  RCP<MV> Y, Z, S, T;
484  S = MVT::Clone( *R_, numRHS_ );
485  T = MVT::Clone( *R_, numRHS_ );
486  if (lp_->isLeftPrec() || lp_->isRightPrec()) {
487  Y = MVT::Clone( *R_, numRHS_ );
488  Z = MVT::Clone( *R_, numRHS_ );
489  }
490  else {
491  Y = P_;
492  Z = S;
493  }
494 
495  // Get the current solution std::vector.
496  Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
497 
499  // Iterate until the status test tells us to stop.
500  //
501  while (stest_->checkStatus(this) != Passed && !breakdown_) {
502 
503  // Increment the iteration
504  iter_++;
505 
506  // rho_new = <R_, Rhat_>
507  MVT::MvDot(*R_,*Rhat_,rho_new);
508 
509  // beta = ( rho_new / rho_old ) (alpha / omega )
510  // TODO: None of these loops are currently threaded
511  for(i=0; i<numRHS_; i++) {
512  // Catch breakdown in rho_old here, since
513  // it is just rho_new from the previous iteration.
514  if (SCT::magnitude(rho_new[i]) < MT::sfmin())
515  breakdown_ = true;
516 
517  beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
518  }
519 
520  // p = r + beta (p - omega v)
521  // TODO: Is it safe to call MvAddMv with A or B = mv?
522  // TODO: Not all of these things have to be part of the state
523  axpy(one, *P_, omega_, *V_, *P_, true); // p = p - omega v
524  axpy(one, *R_, beta, *P_, *P_); // p = r + beta (p - omega v)
525 
526  // y = K\p, unless K does not exist
527  // TODO: There may be a more efficient way to apply the preconditioners
528  if(lp_->isLeftPrec()) {
529  if(lp_->isRightPrec()) {
530  if(leftPrecVec == Teuchos::null) {
531  leftPrecVec = MVT::Clone( *R_, numRHS_ );
532  }
533  lp_->applyLeftPrec(*P_,*leftPrecVec);
534  lp_->applyRightPrec(*leftPrecVec,*Y);
535  }
536  else {
537  lp_->applyLeftPrec(*P_,*Y);
538  }
539  }
540  else if(lp_->isRightPrec()) {
541  lp_->applyRightPrec(*P_,*Y);
542  }
543 
544  // v = Ay
545  lp_->applyOp(*Y,*V_);
546 
547  // alpha = rho_new / <Rhat, V>
548  MVT::MvDot(*V_,*Rhat_,rhatV);
549  for(i=0; i<numRHS_; i++) {
550  if (SCT::magnitude(rhatV[i]) < MT::sfmin())
551  {
552  breakdown_ = true;
553  return;
554  }
555  else
556  alpha_[i] = rho_new[i] / rhatV[i];
557  }
558 
559  // s = r - alpha v
560  axpy(one, *R_, alpha_, *V_, *S, true);
561 
562  // z = K\s, unless K does not exist
563  if(lp_->isLeftPrec()) {
564  if(lp_->isRightPrec()) {
565  if(leftPrecVec == Teuchos::null) {
566  leftPrecVec = MVT::Clone( *R_, numRHS_ );
567  }
568  lp_->applyLeftPrec(*S,*leftPrecVec);
569  lp_->applyRightPrec(*leftPrecVec,*Z);
570  }
571  else {
572  lp_->applyLeftPrec(*S,*Z);
573  }
574  }
575  else if(lp_->isRightPrec()) {
576  lp_->applyRightPrec(*S,*Z);
577  }
578 
579  // t = Az
580  lp_->applyOp(*Z,*T);
581 
582  // omega = <K1\t,K1\s> / <K1\t,K1\t>
583  if(lp_->isLeftPrec()) {
584  if(leftPrecVec == Teuchos::null) {
585  leftPrecVec = MVT::Clone( *R_, numRHS_ );
586  }
587  if(leftPrecVec2 == Teuchos::null) {
588  leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
589  }
590  lp_->applyLeftPrec(*T,*leftPrecVec2);
591  MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
592  MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
593  }
594  else {
595  MVT::MvDot(*T,*T,tT);
596  MVT::MvDot(*T,*S,tS);
597  }
598  for(i=0; i<numRHS_; i++) {
599  if (SCT::magnitude(tT[i]) < MT::sfmin())
600  {
601  omega_[i] = SCT::zero();
602  breakdown_ = true;
603  }
604  else
605  omega_[i] = tS[i] / tT[i];
606  }
607 
608  // x = x + alpha y + omega z
609  axpy(one, *X, alpha_, *Y, *X); // x = x + alpha y
610  axpy(one, *X, omega_, *Z, *X); // x = x + alpha y + omega z
611 
612  // r = s - omega t
613  axpy(one, *S, omega_, *T, *R_, true);
614 
615  // Update rho_old
616  rho_old_ = rho_new;
617  } // end while (sTest_->checkStatus(this) != Passed)
618  }
619 
620 
622  // Iterate until the status test informs us we should stop.
623  template <class ScalarType, class MV, class OP>
624  void BiCGStabIter<ScalarType,MV,OP>::axpy(const ScalarType alpha, const MV & A,
625  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus)
626  {
627  Teuchos::RCP<const MV> A1, B1;
628  Teuchos::RCP<MV> mv1;
629  std::vector<int> index(1);
630 
631  for(int i=0; i<numRHS_; i++) {
632  index[0] = i;
633  A1 = MVT::CloneView(A,index);
634  B1 = MVT::CloneView(B,index);
635  mv1 = MVT::CloneViewNonConst(mv,index);
636  if(minus) {
637  MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
638  }
639  else {
640  MVT::MvAddMv(alpha,*A1,beta[i],*B1,*mv1);
641  }
642  }
643  }
644 
645 } // end Belos namespace
646 
647 #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.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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...
Pure virtual base class for defining the status testing capabilities of Belos.
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.
int getNumIters() const
Get the current iteration count.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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 > P
The first decent direction vector.
BiCGStabIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
SCT::magnitudeType MagnitudeType
MultiVecTraits< ScalarType, MV > MVT
bool breakdownDetected()
Has breakdown been detected in any linear system.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Belos header file which uses auto-configuration information to include necessary C++ headers...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated for Belos by doxygen 1.8.14