Belos  Version of the Day
 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 
186  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
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()
247  bool initialized_;
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.
252  bool stateStorageInitialized_;
253 
254  // Current number of iterations performed.
255  int iter_;
256 
257  // Should the allreduce for convergence detection be merged with the inner products?
258  bool foldConvergenceDetectionIntoAllreduce_;
259 
260  // <r,z>
261  ScalarType rHz_;
262  // <r,r>
263  ScalarType rHr_;
264 
265  //
266  // State Storage
267  //
268  // Residual
269  Teuchos::RCP<MV> R_;
270  // Preconditioned residual
271  Teuchos::RCP<MV> Z_;
272  // Operator applied to preconditioned residual
273  Teuchos::RCP<MV> AZ_;
274  // Direction vector
275  Teuchos::RCP<MV> P_;
276  // Operator applied to direction vector
277  Teuchos::RCP<MV> AP_;
278  //
279  // W_ = (R_, AZ_, Z_)
280  Teuchos::RCP<MV> W_;
281  // S_ = (R_, AZ_)
282  Teuchos::RCP<MV> S_;
283  // T_ = (Z_, R_)
284  Teuchos::RCP<MV> T_;
285  // U_ = (AZ_, Z_)
286  Teuchos::RCP<MV> U_;
287  // V_ = (AP_, P_)
288  Teuchos::RCP<MV> V_;
289 
290 };
291 
293  // Constructor.
294  template<class ScalarType, class MV, class OP>
296  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
299  Teuchos::ParameterList &params ):
300  lp_(problem),
301  om_(printer),
302  stest_(tester),
303  convTest_(convTester),
304  initialized_(false),
305  stateStorageInitialized_(false),
306  iter_(0)
307  {
308  foldConvergenceDetectionIntoAllreduce_ = params.get<bool>("Fold Convergence Detection Into Allreduce",false);
309  }
310 
312  // Setup the state storage.
313  template <class ScalarType, class MV, class OP>
315  {
316  if (!stateStorageInitialized_) {
317 
318  // Check if there is any multivector to clone from.
319  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
320  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
321  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
322  stateStorageInitialized_ = false;
323  return;
324  }
325  else {
326 
327  // Initialize the state storage
328  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
329  if (R_ == Teuchos::null) {
330  // Get the multivector that is not null.
331  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
332  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
333  "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
334 
335  // W_ = (R_, AZ_, Z_)
336  W_ = MVT::Clone( *tmp, 3 );
337  std::vector<int> index2(2,0);
338  std::vector<int> index(1,0);
339 
340  // S_ = (R_, AZ_)
341  index2[0] = 0;
342  index2[1] = 1;
343  S_ = MVT::CloneViewNonConst( *W_, index2 );
344 
345  // U_ = (AZ_, Z_)
346  index2[0] = 1;
347  index2[1] = 2;
348  U_ = MVT::CloneViewNonConst( *W_, index2 );
349 
350  index[0] = 0;
351  R_ = MVT::CloneViewNonConst( *W_, index );
352  index[0] = 1;
353  AZ_ = MVT::CloneViewNonConst( *W_, index );
354  index[0] = 2;
355  Z_ = MVT::CloneViewNonConst( *W_, index );
356 
357  // T_ = (R_, Z_)
358  index2[0] = 0;
359  index2[1] = 2;
360  T_ = MVT::CloneViewNonConst( *W_, index2 );
361 
362  // V_ = (AP_, P_)
363  V_ = MVT::Clone( *tmp, 2 );
364  index[0] = 0;
365  AP_ = MVT::CloneViewNonConst( *V_, index );
366  index[0] = 1;
367  P_ = MVT::CloneViewNonConst( *V_, index );
368 
369  }
370 
371  // State storage has now been initialized.
372  stateStorageInitialized_ = true;
373  }
374  }
375  }
376 
377 
379  // Initialize this iteration object
380  template <class ScalarType, class MV, class OP>
382  {
383  // Initialize the state storage if it isn't already.
384  if (!stateStorageInitialized_)
385  setStateSize();
386 
387  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
388  "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
389 
390  // NOTE: In CGSingleRedIter R_, the initial residual, is required!!!
391  //
392  std::string errstr("Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
393 
394  if (newstate.R != Teuchos::null) {
395 
396  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
397  std::invalid_argument, errstr );
398  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
399  std::invalid_argument, errstr );
400 
401  // Copy basis vectors from newstate into V
402  if (newstate.R != R_) {
403  // copy over the initial residual (unpreconditioned).
404  MVT::Assign( *newstate.R, *R_ );
405  }
406 
407  // Compute initial direction vectors
408  // Initially, they are set to the preconditioned residuals
409  //
410  if ( lp_->getLeftPrec() != Teuchos::null ) {
411  lp_->applyLeftPrec( *R_, *Z_ );
412  if ( lp_->getRightPrec() != Teuchos::null ) {
413  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
414  lp_->applyRightPrec( *Z_, *tmp );
415  MVT::Assign( *tmp, *Z_ );
416  }
417  }
418  else if ( lp_->getRightPrec() != Teuchos::null ) {
419  lp_->applyRightPrec( *R_, *Z_ );
420  }
421  else {
422  MVT::Assign( *R_, *Z_ );
423  }
424 
425  // Multiply the current preconditioned residual vector by A and store in AZ_
426  lp_->applyOp( *Z_, *AZ_ );
427 
428  // P_ := Z_
429  // Logically, AP_ := AZ_
430  MVT::Assign( *U_, *V_);
431  }
432  else {
433 
434  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
435  "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
436  }
437 
438  // The solver is initialized
439  initialized_ = true;
440  }
441 
442 
444  // Get the norms of the residuals native to the solver.
445  template <class ScalarType, class MV, class OP>
447  CGSingleRedIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
448  if (convTest_->getResNormType() == Belos::PreconditionerNorm) {
449  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHz_));
450  return Teuchos::null;
451  } else if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
452  (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
453  return Teuchos::null;
454  } else
455  return R_;
456  }
457 
458 
460  // Iterate until the status test informs us we should stop.
461  template <class ScalarType, class MV, class OP>
463  {
464  //
465  // Allocate/initialize data structures
466  //
467  if (initialized_ == false) {
468  initialize();
469  }
470 
471  // Allocate memory for scalars.
474  ScalarType rHz_old, alpha, beta, delta;
475 
476  // Create convenience variables for zero and one.
477  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
479 
480  // Get the current solution vector.
481  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
482 
483  // Check that the current solution vector only has one column.
484  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
485  "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
486 
487  if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
488  // Compute first <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R_> combined (also computes unneeded <AZ_,R_>)
489  MVT::MvTransMv( one, *S_, *T_, sHt );
490  rHz_ = sHt(0,1);
491  delta = sHt(1,1);
492  rHr_ = sHt(0,0);
493  } else {
494  // Compute first <s,z> a.k.a. <r,z> and <Az,z> combined
495  MVT::MvTransMv( one, *S_, *Z_, sHz );
496  rHz_ = sHz(0,0);
497  delta = sHz(1,0);
498  }
500  (stest_->checkStatus(this) == Passed))
501  return;
502  alpha = rHz_ / delta;
503 
504  // Check that alpha is a positive number!
505  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
506  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
507 
509  // Iterate until the status test tells us to stop.
510  //
511  if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
513  // Iterate until the status test tells us to stop.
514  //
515  while (1) {
516 
517  // Update the solution vector x := x + alpha * P_
518  //
519  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
520  //
521  // Compute the new residual R_ := R_ - alpha * AP_
522  //
523  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
524  //
525  // Apply preconditioner to new residual to update Z_
526  //
527  if ( lp_->getLeftPrec() != Teuchos::null ) {
528  lp_->applyLeftPrec( *R_, *Z_ );
529  if ( lp_->getRightPrec() != Teuchos::null ) {
530  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
531  lp_->applyRightPrec( *tmp, *Z_ );
532  }
533  }
534  else if ( lp_->getRightPrec() != Teuchos::null ) {
535  lp_->applyRightPrec( *R_, *Z_ );
536  }
537  else {
538  MVT::Assign( *R_, *Z_ );
539  }
540  //
541  // Multiply the current preconditioned residual vector by A and store in AZ_
542  lp_->applyOp( *Z_, *AZ_ );
543  //
544  // Compute <S_,T_> a.k.a. <R_,Z_>, <AZ_,Z_> and <R_,R_> combined (also computes unneeded <AZ_,R_>)
545  MVT::MvTransMv( one, *S_, *T_, sHt );
546  //
547  // Update scalars.
548  rHz_old = rHz_;
549  rHz_ = sHt(0,1);
550  delta = sHt(1,1);
551  rHr_ = sHt(0,0);
552 
553  // Increment the iteration
554  iter_++;
555  //
556  // Check the status test, now that the solution and residual have been updated
557  //
558  if (stest_->checkStatus(this) == Passed) {
559  break;
560  }
561  //
562  beta = rHz_ / rHz_old;
563  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
564  //
565  // Check that alpha is a positive number!
566  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
567  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
568 
569  //
570  // Update the direction vector P_ := Z_ + beta * P_
571  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
572  // Hence: V_ = (AP_, P_) := (AZ_, Z_) + beta (AP_, P_) = U_ + beta * V_
573  //
574  MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
575 
576  } // end while (1)
577  } else {
579  // Iterate until the status test tells us to stop.
580  //
581  while (1) {
582 
583  // Update the solution vector x := x + alpha * P_
584  //
585  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
586  //
587  // Compute the new residual R_ := R_ - alpha * AP_
588  //
589  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
590 
591  // Increment the iteration
592  iter_++;
593  //
594  // Check the status test, now that the solution and residual have been updated
595  //
596  if (stest_->checkStatus(this) == Passed) {
597  break;
598  }
599  //
600  // Apply preconditioner to new residual to update Z_
601  //
602  if ( lp_->getLeftPrec() != Teuchos::null ) {
603  lp_->applyLeftPrec( *R_, *Z_ );
604  if ( lp_->getRightPrec() != Teuchos::null ) {
605  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
606  lp_->applyRightPrec( *tmp, *Z_ );
607  }
608  }
609  else if ( lp_->getRightPrec() != Teuchos::null ) {
610  lp_->applyRightPrec( *R_, *Z_ );
611  }
612  else {
613  MVT::Assign( *R_, *Z_ );
614  }
615  //
616  // Multiply the current preconditioned residual vector by A and store in AZ_
617  lp_->applyOp( *Z_, *AZ_ );
618  //
619  // Compute <S_,Z_> a.k.a. <R_,Z_> and <AZ_,Z_> combined
620  MVT::MvTransMv( one, *S_, *Z_, sHz );
621  //
622  // Update scalars.
623  rHz_old = rHz_;
624  rHz_ = sHz(0,0);
625  delta = sHz(1,0);
626  //
627  beta = rHz_ / rHz_old;
628  alpha = rHz_ / (delta - (beta*rHz_ / alpha));
629  //
630  // Check that alpha is a positive number!
631  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGPositiveDefiniteFailure,
632  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
633 
634  //
635  // Update the direction vector P_ := Z_ + beta * P_
636  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
637  // Hence: V_ = (AP_, P_) := (AZ_, Z_) + beta (AP_, P_) = U_ + beta * V_
638  //
639  MVT::MvAddMv( one, *U_, beta, *V_, *V_ );
640 
641  } // end while (1)
642  }
643  }
644 
645 } // end Belos namespace
646 
647 #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.
int getNumIters() const
Get the current iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
T & get(const std::string &name, T def_value)
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.
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.
CGPositiveDefiniteFailure is thrown when the the CG &#39;alpha = p^H*A*P&#39; value is less than zero...
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
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.
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...
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Class which defines basic traits for the operator type.
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)

Generated on Thu Mar 28 2024 09:24:28 for Belos by doxygen 1.8.5