Belos Package Browser (Single Doxygen Collection)  Development
 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 //
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"
63 #include "Teuchos_ScalarTraits.hpp"
65 #include "Teuchos_TimeMonitor.hpp"
66 
78 namespace Belos {
79 
81 
82 
87  template <class ScalarType, class MV>
89 
92 
95 
98 
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  //
125 
127 
128 
135  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
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?
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  //
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()
282 
283  // Current number of iterations performed.
284  int iter_;
285 
286  //
287  // State Storage
288  //
289  // Initial residual
291  //
292  // Residual
294  //
295  // Direction vector 1
297  //
298  // Operator applied to preconditioned direction vector 1
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,
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(*R_, *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  {
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.
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.
std::vector< ScalarType > rho_old_
Class which defines basic traits for the operator type.
std::vector< ScalarType > alpha
const Teuchos::RCP< OutputManager< ScalarType > > om_
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > Rhat_
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...
void axpy(const ScalarType alpha, const MV &A, const std::vector< ScalarType > beta, const MV &B, MV &mv, bool minus=false)
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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 > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< MV > V_
TransListIter iter
Teuchos::RCP< const MV > P
The first decent direction vector.
Teuchos::RCP< MV > R_
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.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
int getNumIters() const
Get the current iteration count.
std::vector< ScalarType > alpha_
SCT::magnitudeType MagnitudeType
MultiVecTraits< ScalarType, MV > MVT
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.
Teuchos::RCP< MV > P_
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
std::vector< ScalarType > omega_