Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosCGSingleRedIter.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_CG_SINGLE_RED_ITER_HPP
43 #define BELOS_CG_SINGLE_RED_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosCGIteration.hpp"
52 
53 #include "BelosLinearProblem.hpp"
54 #include "BelosOutputManager.hpp"
55 #include "BelosStatusTest.hpp"
57 #include "BelosOperatorTraits.hpp"
58 #include "BelosMultiVecTraits.hpp"
59 
62 #include "Teuchos_ScalarTraits.hpp"
64 #include "Teuchos_TimeMonitor.hpp"
65 
76 namespace Belos {
77 
78 template<class ScalarType, class MV, class OP>
79 class CGSingleRedIter : virtual public CGIteration<ScalarType,MV,OP> {
80 
81  public:
82 
83  //
84  // Convenience typedefs
85  //
90 
92 
93 
100  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
103  Teuchos::ParameterList &params );
104 
106  virtual ~CGSingleRedIter() {};
108 
109 
111 
112 
125  void iterate();
126 
142 
146  void initialize()
147  {
149  initializeCG(empty);
150  }
151 
160  state.R = R_;
161  state.P = P_;
162  state.AP = AP_;
163  state.Z = Z_;
164  return state;
165  }
166 
168 
169 
171 
172 
174  int getNumIters() const { return iter_; }
175 
177  void resetNumIters( int iter = 0 ) { iter_ = iter; }
178 
181  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
182 
184 
187 
189 
191 
192 
194  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
195 
197  int getBlockSize() const { return 1; }
198 
200  void setBlockSize(int blockSize) {
201  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
202  "Belos::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
203  }
204 
206  bool isInitialized() { return initialized_; }
207 
209  void setDoCondEst(bool /* val */){/*ignored*/}
210 
214  return temp;
215  }
216 
220  return temp;
221  }
222 
224 
225  private:
226 
227  //
228  // Internal methods
229  //
231  void setStateSize();
232 
233  //
234  // Classes inputed through constructor that define the linear problem to be solved.
235  //
240 
241  //
242  // Current solver state
243  //
244  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
245  // is capable of running; _initialize is controlled by the initialize() member method
246  // For the implications of the state of initialized_, please see documentation for initialize()
248 
249  // stateStorageInitialized_ specifies that the state storage has been initialized.
250  // This initialization may be postponed if the linear problem was generated without
251  // the right-hand side or solution vectors.
253 
254  // Current number of iterations performed.
255  int iter_;
256 
257  // Should the allreduce for convergence detection be merged with the inner products?
259 
260  // <r,z>
261  ScalarType rHz_;
262  // <r,r>
263  ScalarType rHr_;
264 
265  //
266  // State Storage
267  //
268  // Residual
270  //
271  // Preconditioned residual
273  //
274  // Operator applied to preconditioned residual
276  //
277  // This is the additional storage needed for single-reduction CG (Chronopoulos/Gear)
278  // R_ views the first vector and AZ_ views the second vector
280 
281  // This is the additional storage needed for merging convergence detection and inner products
282  // (Option "Fold Convergence Detection Into Allreduce")
283  // Z_ views the first vector and R2_ views the second
286  //
287  // Direction vector
289  //
290  // Operator applied to direction vector (S)
292 
293 };
294 
296  // Constructor.
297  template<class ScalarType, class MV, class OP>
299  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
302  Teuchos::ParameterList &params ):
303  lp_(problem),
304  om_(printer),
305  stest_(tester),
306  convTest_(convTester),
307  initialized_(false),
308  stateStorageInitialized_(false),
309  iter_(0)
310  {
311  foldConvergenceDetectionIntoAllreduce_ = params.get<bool>("Fold Convergence Detection Into Allreduce",false);
312  }
313 
315  // Setup the state storage.
316  template <class ScalarType, class MV, class OP>
318  {
319  if (!stateStorageInitialized_) {
320 
321  // Check if there is any multivector to clone from.
322  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
323  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
324  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
325  stateStorageInitialized_ = false;
326  return;
327  }
328  else {
329 
330  // Initialize the state storage
331  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
332  if (R_ == Teuchos::null) {
333  // Get the multivector that is not null.
334  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
335  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
336  "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
337  S_ = MVT::Clone( *tmp, 2 );
338  if (foldConvergenceDetectionIntoAllreduce_) {
339  T_ = MVT::Clone( *tmp, 2 );
340  // Z_ will view the first column of T_, R2_ will view the second.
341  std::vector<int> index(1,0);
342  Z_ = MVT::CloneViewNonConst( *T_, index );
343  index[0] = 1;
344  R2_ = MVT::CloneViewNonConst( *T_, index );
345  }
346  else
347  Z_ = MVT::Clone( *tmp, 1 );
348  P_ = MVT::Clone( *tmp, 1 );
349  AP_ = MVT::Clone( *tmp, 1 );
350 
351  // R_ will view the first column of S_, AZ_ will view the second.
352  std::vector<int> index(1,0);
353  R_ = MVT::CloneViewNonConst( *S_, index );
354  index[0] = 1;
355  AZ_ = MVT::CloneViewNonConst( *S_, index );
356  }
357 
358  // State storage has now been initialized.
359  stateStorageInitialized_ = true;
360  }
361  }
362  }
363 
364 
366  // Initialize this iteration object
367  template <class ScalarType, class MV, class OP>
369  {
370  // Initialize the state storage if it isn't already.
371  if (!stateStorageInitialized_)
372  setStateSize();
373 
374  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
375  "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
376 
377  // NOTE: In CGSingleRedIter R_, the initial residual, is required!!!
378  //
379  std::string errstr("Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
380 
381  if (newstate.R != Teuchos::null) {
382 
383  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
384  std::invalid_argument, errstr );
385  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
386  std::invalid_argument, errstr );
387 
388  // Copy basis vectors from newstate into V
389  if (newstate.R != R_) {
390  // copy over the initial residual (unpreconditioned).
391  MVT::Assign( *newstate.R, *R_ );
392  }
393 
394  // Compute initial direction vectors
395  // Initially, they are set to the preconditioned residuals
396  //
397  if ( lp_->getLeftPrec() != Teuchos::null ) {
398  lp_->applyLeftPrec( *R_, *Z_ );
399  if ( lp_->getRightPrec() != Teuchos::null ) {
400  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
401  lp_->applyRightPrec( *Z_, *tmp );
402  MVT::Assign( *tmp, *Z_ );
403  }
404  }
405  else if ( lp_->getRightPrec() != Teuchos::null ) {
406  lp_->applyRightPrec( *R_, *Z_ );
407  }
408  else {
409  MVT::Assign( *R_, *Z_ );
410  }
411  MVT::Assign( *Z_, *P_ );
412 
413  // Multiply the current preconditioned residual vector by A and store in AZ_
414  lp_->applyOp( *Z_, *AZ_ );
415 
416  // Logically, AP_ = AZ_
417  MVT::Assign( *AZ_, *AP_ );
418  }
419  else {
420 
421  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
422  "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
423  }
424 
425  // The solver is initialized
426  initialized_ = true;
427  }
428 
429 
431  // Get the norms of the residuals native to the solver.
432  template <class ScalarType, class MV, class OP>
434  CGSingleRedIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
435  if (convTest_->getResNormType() == Belos::PreconditionerNorm) {
436  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHz_));
437  return Teuchos::null;
438  } else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
439  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
440  return Teuchos::null;
441  } else
442  return R_;
443  }
444 
445 
447  // Iterate until the status test informs us we should stop.
448  template <class ScalarType, class MV, class OP>
450  {
451  //
452  // Allocate/initialize data structures
453  //
454  if (initialized_ == false) {
455  initialize();
456  }
457 
458  // Allocate memory for scalars.
461  ScalarType rHz_old, alpha, beta, delta;
462 
463  // Create convenience variables for zero and one.
464  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
466 
467  // Get the current solution vector.
468  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
469 
470  // Check that the current solution vector only has one column.
471  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
472  "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
473 
474  if (foldConvergenceDetectionIntoAllreduce_) {
475  // Compute first <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R2_> combined (also computes unneeded <AZ_,R2_>)
476  MVT::Assign( *R_, *R2_ );
477  MVT::MvTransMv( one, *S_, *T_, sHt );
478  rHz_ = sHt(0,0);
479  delta = sHt(1,0);
480  rHr_ = sHt(1,1);
481  } else {
482  // Compute first <s,z> a.k.a. <r,z> and <Az,z> combined
483  MVT::MvTransMv( one, *S_, *Z_, sHz );
484  rHz_ = sHz(0,0);
485  delta = sHz(1,0);
486  }
488  (stest_->checkStatus(this) == Passed))
489  return;
490  alpha = rHz_ / delta;
491 
492  // Check that alpha is a positive number!
493  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGIterateFailure,
494  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
495 
497  // Iterate until the status test tells us to stop.
498  //
499  if (foldConvergenceDetectionIntoAllreduce_) {
501  // Iterate until the status test tells us to stop.
502  //
503  while (1) {
504 
505  // Update the solution vector x := x + alpha * P_
506  //
507  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
508  //
509  // Compute the new residual R_ := R_ - alpha * AP_
510  //
511  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
512  //
513  // Apply preconditioner to new residual to update Z_
514  //
515  if ( lp_->getLeftPrec() != Teuchos::null ) {
516  lp_->applyLeftPrec( *R_, *Z_ );
517  if ( lp_->getRightPrec() != Teuchos::null ) {
518  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
519  lp_->applyRightPrec( *Z_, *tmp );
520  MVT::Assign( *tmp, *Z_ );
521  }
522  }
523  else if ( lp_->getRightPrec() != Teuchos::null ) {
524  lp_->applyRightPrec( *R_, *Z_ );
525  }
526  else {
527  MVT::Assign( *R_, *Z_ );
528  }
529  //
530  // Multiply the current preconditioned residual vector by A and store in AZ_
531  lp_->applyOp( *Z_, *AZ_ );
532  //
533  // Compute <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R2_> combined (also computes unneeded <AZ_,R2_>)
534  MVT::Assign( *R_, *R2_ );
535  MVT::MvTransMv( one, *S_, *T_, sHt );
536  //
537  // Update scalars.
538  rHz_old = rHz_;
539  rHz_ = sHt(0,0);
540  delta = sHt(1,0);
541  rHr_ = sHt(1,1);
542  //
543  // Check the status test, now that the solution and residual have been updated
544  //
545  if (stest_->checkStatus(this) == Passed) {
546  break;
547  }
548  // Increment the iteration
549  iter_++;
550  //
551  beta = rHz_ / rHz_old;
552  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
553  //
554  // Check that alpha is a positive number!
555  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGIterateFailure,
556  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
557  //
558  // Update the direction vector P_ := Z_ + beta * P_
559  //
560  MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
561  //
562  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
563  // NOTE: This increases the number of vector updates by 1, in exchange for
564  // reducing the collectives from 2 to 1.
565  //
566  MVT::MvAddMv( one, *AZ_, beta, *AP_, *AP_ );
567  //
568  } // end while (1)
569  } else {
571  // Iterate until the status test tells us to stop.
572  //
573  while (1) {
574 
575  // Update the solution vector x := x + alpha * P_
576  //
577  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
578  //
579  // Compute the new residual R_ := R_ - alpha * AP_
580  //
581  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
582  //
583  // Check the status test, now that the solution and residual have been updated
584  //
585  if (stest_->checkStatus(this) == Passed) {
586  break;
587  }
588  // Increment the iteration
589  iter_++;
590  //
591  // Apply preconditioner to new residual to update Z_
592  //
593  if ( lp_->getLeftPrec() != Teuchos::null ) {
594  lp_->applyLeftPrec( *R_, *Z_ );
595  if ( lp_->getRightPrec() != Teuchos::null ) {
596  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
597  lp_->applyRightPrec( *Z_, *tmp );
598  Z_ = tmp;
599  }
600  }
601  else if ( lp_->getRightPrec() != Teuchos::null ) {
602  lp_->applyRightPrec( *R_, *Z_ );
603  }
604  else {
605  Z_ = R_;
606  }
607  //
608  // Multiply the current preconditioned residual vector by A and store in AZ_
609  lp_->applyOp( *Z_, *AZ_ );
610  //
611  // Compute <S_,Z_> a.k.a. <R_,Z_> and <AZ_,Z_> combined
612  MVT::MvTransMv( one, *S_, *Z_, sHz );
613  //
614  // Update scalars.
615  rHz_old = rHz_;
616  rHz_ = sHz(0,0);
617  delta = sHz(1,0);
618  //
619  beta = rHz_ / rHz_old;
620  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
621  //
622  // Check that alpha is a positive number!
623  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGIterateFailure,
624  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
625  //
626  // Update the direction vector P_ := Z_ + beta * P_
627  //
628  MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
629  //
630  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
631  // NOTE: This increases the number of vector updates by 1, in exchange for
632  // reducing the collectives from 2 to 1.
633  //
634  MVT::MvAddMv( one, *AZ_, beta, *AP_, *AP_ );
635  //
636  } // end while (1)
637  }
638  }
639 
640 } // end Belos namespace
641 
642 #endif /* BELOS_CG_SINGLE_RED_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void setStateSize()
Method for initalizing the state storage needed by CG.
int getNumIters() const
Get the current iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
T & get(ParameterList &l, const std::string &name)
Structure to contain pointers to CGIteration state variables.
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
An implementation of StatusTestResNorm using a family of residual norms.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
MultiVecTraits< ScalarType, MV > MVT
Traits class which defines basic operations on multivectors.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
bool isInitialized()
States whether the solver has been initialized or not.
virtual ~CGSingleRedIter()
Destructor.
SCT::magnitudeType MagnitudeType
const Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
TransListIter iter
CGSingleRedIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > &convTester, Teuchos::ParameterList &params)
CGSingleRedIter constructor with linear problem, solver utilities, and parameter list of solver optio...
const Teuchos::RCP< OutputManager< ScalarType > > om_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Class which defines basic traits for the operator type.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Teuchos::ScalarTraits< ScalarType > SCT
Belos header file which uses auto-configuration information to include necessary C++ headers...
void resetNumIters(int iter=0)
Reset the iteration count.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > P
The current decent direction vector.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void setDoCondEst(bool)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation (NOT_IMPLEMENTED)