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_BLAS.hpp"
61 #include "Teuchos_SerialDenseMatrix.hpp"
62 #include "Teuchos_SerialDenseVector.hpp"
63 #include "Teuchos_ScalarTraits.hpp"
64 #include "Teuchos_ParameterList.hpp"
65 #include "Teuchos_TimeMonitor.hpp"
66 
78 namespace Belos {
79 
81 
82 
87  template <class ScalarType, class MV>
89 
91  Teuchos::RCP<const MV> R;
92 
94  Teuchos::RCP<const MV> Rhat;
95 
97  Teuchos::RCP<const MV> P;
98 
100  Teuchos::RCP<const MV> V;
101 
102  std::vector<ScalarType> rho_old, alpha, omega;
103 
104  BiCGStabIterationState() : R(Teuchos::null), Rhat(Teuchos::null),
105  P(Teuchos::null), V(Teuchos::null)
106  {
107  rho_old.clear();
108  alpha.clear();
109  omega.clear();
110  }
111  };
112 
113  template<class ScalarType, class MV, class OP>
114  class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP> {
115 
116  public:
117 
118  //
119  // Convenience typedefs
120  //
123  typedef Teuchos::ScalarTraits<ScalarType> SCT;
124  typedef typename SCT::magnitudeType MagnitudeType;
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 
238 
239 
241  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
242 
244  int getBlockSize() const { return 1; }
245 
247  void setBlockSize(int blockSize) {
248  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
249  "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
250  }
251 
253  bool isInitialized() { return initialized_; }
254 
256 
257  private:
258 
259  void axpy(const ScalarType alpha, const MV & A,
260  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus=false);
261 
262  //
263  // Classes inputed through constructor that define the linear problem to be solved.
264  //
265  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
266  const Teuchos::RCP<OutputManager<ScalarType> > om_;
267  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
268 
269  //
270  // Algorithmic parameters
271  //
272  // numRHS_ is the current number of linear systems being solved.
273  int numRHS_;
274 
275  //
276  // Current solver state
277  //
278  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
279  // is capable of running; _initialize is controlled by the initialize() member method
280  // For the implications of the state of initialized_, please see documentation for initialize()
281  bool initialized_;
282 
283  // Current number of iterations performed.
284  int iter_;
285 
286  //
287  // State Storage
288  //
289  // Initial residual
290  Teuchos::RCP<MV> Rhat_;
291  //
292  // Residual
293  Teuchos::RCP<MV> R_;
294  //
295  // Direction vector 1
296  Teuchos::RCP<MV> P_;
297  //
298  // Operator applied to preconditioned direction vector 1
299  Teuchos::RCP<MV> V_;
300  //
301  std::vector<ScalarType> rho_old_, alpha_, omega_;
302  };
303 
305  // Constructor.
306  template<class ScalarType, class MV, class OP>
308  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
309  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
310  Teuchos::ParameterList &params ):
311  lp_(problem),
312  om_(printer),
313  stest_(tester),
314  numRHS_(0),
315  initialized_(false),
316  iter_(0)
317  {
318  }
319 
320 
322  // Initialize this iteration object
323  template <class ScalarType, class MV, class OP>
325  {
326  // Check if there is any multivector to clone from.
327  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
328  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
329  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
330  "Belos::BiCGStabIter::initialize(): Cannot initialize state storage!");
331 
332  // Get the multivector that is not null.
333  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
334 
335  // Get the number of right-hand sides we're solving for now.
336  int numRHS = MVT::GetNumberVecs(*tmp);
337  numRHS_ = numRHS;
338 
339  // Initialize the state storage
340  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
341  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
342  R_ = MVT::Clone( *tmp, numRHS_ );
343  Rhat_ = MVT::Clone( *tmp, numRHS_ );
344  P_ = MVT::Clone( *tmp, numRHS_ );
345  V_ = MVT::Clone( *tmp, numRHS_ );
346 
347  rho_old_.resize(numRHS_);
348  alpha_.resize(numRHS_);
349  omega_.resize(numRHS_);
350  }
351 
352  // NOTE: In BiCGStabIter R_, the initial residual, is required!!!
353  //
354  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
355 
356  // Create convenience variable for one.
357  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
358 
359  if (!Teuchos::is_null(newstate.R)) {
360 
361  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
362  std::invalid_argument, errstr );
363  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
364  std::invalid_argument, errstr );
365 
366  // Copy residual vectors from newstate into R
367  if (newstate.R != R_) {
368  // Assigned by the new state
369  MVT::Assign(*newstate.R, *R_);
370  }
371  else {
372  // Computed
373  lp_->computeCurrResVec(R_.get());
374  }
375 
376  // Set Rhat
377  if (!Teuchos::is_null(newstate.Rhat) && newstate.Rhat != Rhat_) {
378  // Assigned by the new state
379  MVT::Assign(*newstate.Rhat, *Rhat_);
380  }
381  else {
382  // Set to be the initial residual
383  MVT::Assign(*lp_->getInitResVec(), *Rhat_);
384  }
385 
386  // Set V
387  if (!Teuchos::is_null(newstate.V) && newstate.V != V_) {
388  // Assigned by the new state
389  MVT::Assign(*newstate.V, *V_);
390  }
391  else {
392  // Initial V = 0
393  MVT::MvInit(*V_);
394  }
395 
396  // Set P
397  if (!Teuchos::is_null(newstate.P) && newstate.P != P_) {
398  // Assigned by the new state
399  MVT::Assign(*newstate.P, *P_);
400  }
401  else {
402  // Initial P = 0
403  MVT::MvInit(*P_);
404  }
405 
406  // Set rho_old
407  if (newstate.rho_old.size () == static_cast<size_t> (numRHS_)) {
408  // Assigned by the new state
409  rho_old_ = newstate.rho_old;
410  }
411  else {
412  // Initial rho = 1
413  rho_old_.assign(numRHS_,one);
414  }
415 
416  // Set alpha
417  if (newstate.alpha.size() == static_cast<size_t> (numRHS_)) {
418  // Assigned by the new state
419  alpha_ = newstate.alpha;
420  }
421  else {
422  // Initial rho = 1
423  alpha_.assign(numRHS_,one);
424  }
425 
426  // Set omega
427  if (newstate.omega.size() == static_cast<size_t> (numRHS_)) {
428  // Assigned by the new state
429  omega_ = newstate.omega;
430  }
431  else {
432  // Initial rho = 1
433  omega_.assign(numRHS_,one);
434  }
435 
436  }
437  else {
438 
439  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
440  "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
441  }
442 
443  // The solver is initialized
444  initialized_ = true;
445  }
446 
447 
449  // Iterate until the status test informs us we should stop.
450  template <class ScalarType, class MV, class OP>
452  {
453  using Teuchos::RCP;
454 
455  //
456  // Allocate/initialize data structures
457  //
458  if (initialized_ == false) {
459  initialize();
460  }
461 
462  // Allocate memory for scalars.
463  int i=0;
464  std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
465  std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
466 
467  // Create convenience variable for one.
468  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
469 
470  // TODO: We may currently be using more space than is required
471  RCP<MV> leftPrecVec, leftPrecVec2;
472 
473  RCP<MV> Y, Z, S, T;
474  S = MVT::Clone( *R_, numRHS_ );
475  T = MVT::Clone( *R_, numRHS_ );
476  if (lp_->isLeftPrec() || lp_->isRightPrec()) {
477  Y = MVT::Clone( *R_, numRHS_ );
478  Z = MVT::Clone( *R_, numRHS_ );
479  }
480  else {
481  Y = P_;
482  Z = S;
483  }
484 
485  // Get the current solution std::vector.
486  Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
487 
489  // Iterate until the status test tells us to stop.
490  //
491  while (stest_->checkStatus(this) != Passed) {
492 // std::cout << std::endl;
493 
494  // std::vector<ScalarType> tempResids(numRHS_);
495  // MVT::MvNorm(*R_,tempResids);
496 // for(i=0; i<numRHS_; i++)
497 // std::cout << "r[" << i << "] = " << tempResids[i] << std::endl;
498 
499  // Increment the iteration
500  iter_++;
501 
502  // rho_new = <R_, Rhat_>
503  MVT::MvDot(*R_,*Rhat_,rho_new);
504 
505 // for(i=0; i<numRHS_; i++) {
506 // std::cout << "rho[" << i << "] = " << rho_new[i] << std::endl;
507 // }
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  beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
513 // std::cout << "beta[" << i << "] = " << beta[i] << std::endl;
514  }
515 
516  // p = r + beta (p - omega v)
517  // TODO: Is it safe to call MvAddMv with A or B = mv?
518  // TODO: Not all of these things have to be part of the state
519  axpy(one, *P_, omega_, *V_, *P_, true); // p = p - omega v
520  axpy(one, *R_, beta, *P_, *P_); // p = r + beta (p - omega v)
521 
522  // y = K\p, unless K does not exist
523  // TODO: There may be a more efficient way to apply the preconditioners
524  if(lp_->isLeftPrec()) {
525  if(lp_->isRightPrec()) {
526  if(leftPrecVec == Teuchos::null) {
527  leftPrecVec = MVT::Clone( *R_, numRHS_ );
528  }
529  lp_->applyLeftPrec(*P_,*leftPrecVec);
530  lp_->applyRightPrec(*leftPrecVec,*Y);
531  }
532  else {
533  lp_->applyLeftPrec(*P_,*Y);
534  }
535  }
536  else if(lp_->isRightPrec()) {
537  lp_->applyRightPrec(*P_,*Y);
538  }
539 
540  // v = Ay
541  lp_->applyOp(*Y,*V_);
542 
543  // alpha = rho_new / <Rhat, V>
544  MVT::MvDot(*V_,*Rhat_,rhatV);
545  for(i=0; i<numRHS_; i++) {
546  alpha_[i] = rho_new[i] / rhatV[i];
547  }
548 
549 // for(i=0; i<numRHS_; i++) {
550 // std::cout << "alpha[" << i << "] = " << alpha_[i] << std::endl;
551 // }
552 
553  // s = r - alpha v
554  axpy(one, *R_, alpha_, *V_, *S, true);
555 
556  // z = K\s, unless K does not exist
557  if(lp_->isLeftPrec()) {
558  if(lp_->isRightPrec()) {
559  if(leftPrecVec == Teuchos::null) {
560  leftPrecVec = MVT::Clone( *R_, numRHS_ );
561  }
562  lp_->applyLeftPrec(*S,*leftPrecVec);
563  lp_->applyRightPrec(*leftPrecVec,*Z);
564  }
565  else {
566  lp_->applyLeftPrec(*S,*Z);
567  }
568  }
569  else if(lp_->isRightPrec()) {
570  lp_->applyRightPrec(*S,*Z);
571  }
572 
573  // t = Az
574  lp_->applyOp(*Z,*T);
575 
576 // std::cout << "T:\n";
577 // T->Print(std::cout);
578 
579  // omega = <K1\t,K1\s> / <K1\t,K1\t>
580  if(lp_->isLeftPrec()) {
581  if(leftPrecVec == Teuchos::null) {
582  leftPrecVec = MVT::Clone( *R_, numRHS_ );
583  }
584  if(leftPrecVec2 == Teuchos::null) {
585  leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
586  }
587  lp_->applyLeftPrec(*T,*leftPrecVec2);
588  MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
589  MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
590  }
591  else {
592  MVT::MvDot(*T,*T,tT);
593  MVT::MvDot(*T,*S,tS);
594  }
595  for(i=0; i<numRHS_; i++) {
596  omega_[i] = tS[i] / tT[i];
597  }
598 
599 // for(i=0; i<numRHS_; i++) {
600 // std::cout << "omega[" << i << "] = " << omega_[i] << " = " << tS[i] << " / " << tT[i] << std::endl;
601 // }
602 
603  // x = x + alpha y + omega z
604  axpy(one, *X, alpha_, *Y, *X); // x = x + alpha y
605  axpy(one, *X, omega_, *Z, *X); // x = x + alpha y + omega z
606 
607  // r = s - omega t
608  axpy(one, *S, omega_, *T, *R_, true);
609 
610  // Update rho_old
611  rho_old_ = rho_new;
612  } // end while (sTest_->checkStatus(this) != Passed)
613  }
614 
615 
617  // Iterate until the status test informs us we should stop.
618  template <class ScalarType, class MV, class OP>
619  void BiCGStabIter<ScalarType,MV,OP>::axpy(const ScalarType alpha, const MV & A,
620  const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus)
621  {
622  Teuchos::RCP<const MV> A1, B1;
623  Teuchos::RCP<MV> mv1;
624  std::vector<int> index(1);
625 
626  for(int i=0; i<numRHS_; i++) {
627  index[0] = i;
628  A1 = MVT::CloneView(A,index);
629  B1 = MVT::CloneView(B,index);
630  mv1 = MVT::CloneViewNonConst(mv,index);
631  if(minus) {
632  MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
633  }
634  else {
635  MVT::MvAddMv(alpha,*A1,beta[i],*B1,*mv1);
636  }
637  }
638  }
639 
640 } // end Belos namespace
641 
642 #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.
Class which defines basic traits for the operator type.
std::vector< ScalarType > alpha
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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::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
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