Belos  Version of the Day
BelosCGIter.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_ITER_HPP
43 #define BELOS_CG_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"
56 #include "BelosOperatorTraits.hpp"
57 #include "BelosMultiVecTraits.hpp"
58 
59 #include "Teuchos_SerialDenseMatrix.hpp"
60 #include "Teuchos_SerialDenseVector.hpp"
61 #include "Teuchos_ScalarTraits.hpp"
62 #include "Teuchos_ParameterList.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
75 namespace Belos {
76 
77 template<class ScalarType, class MV, class OP>
78 class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
79 
80  public:
81 
82  //
83  // Convenience typedefs
84  //
87  typedef Teuchos::ScalarTraits<ScalarType> SCT;
88  typedef typename SCT::magnitudeType MagnitudeType;
89 
91 
92 
98  CGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
99  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
100  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
101  Teuchos::ParameterList &params );
102 
104  virtual ~CGIter() {};
106 
107 
109 
110 
123  void iterate();
124 
140 
144  void initialize()
145  {
147  initializeCG(empty);
148  }
149 
158  state.R = R_;
159  state.P = P_;
160  state.AP = AP_;
161  state.Z = Z_;
162  return state;
163  }
164 
166 
167 
169 
170 
172  int getNumIters() const { return iter_; }
173 
175  void resetNumIters( int iter = 0 ) { iter_ = iter; }
176 
179  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
180 
182 
184  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
185 
187 
189 
190 
192  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
193 
195  int getBlockSize() const { return 1; }
196 
198  void setBlockSize(int blockSize) {
199  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
200  "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
201  }
202 
204  bool isInitialized() { return initialized_; }
205 
206 
208  void setDoCondEst(bool val){doCondEst_=val;}
209 
211  Teuchos::ArrayView<MagnitudeType> getDiag() {
212  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
213  // getDiag() didn't actually throw for me in that case, but why
214  // not be cautious?
215  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
216  if (static_cast<size_type> (iter_) >= diag_.size ()) {
217  return diag_ ();
218  } else {
219  return diag_ (0, iter_);
220  }
221  }
222 
224  Teuchos::ArrayView<MagnitudeType> getOffDiag() {
225  // NOTE (mfh 30 Jul 2015) The implementation as I found it
226  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
227  // debug mode) when the maximum number of iterations has been
228  // reached, because iter_ == offdiag_.size() in that case. The
229  // new logic fixes this.
230  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
231  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
232  return offdiag_ ();
233  } else {
234  return offdiag_ (0, iter_);
235  }
236  }
237 
239 
240  private:
241 
242  //
243  // Internal methods
244  //
246  void setStateSize();
247 
248  //
249  // Classes inputed through constructor that define the linear problem to be solved.
250  //
251  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
252  const Teuchos::RCP<OutputManager<ScalarType> > om_;
253  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
254 
255  //
256  // Current solver state
257  //
258  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
259  // is capable of running; _initialize is controlled by the initialize() member method
260  // For the implications of the state of initialized_, please see documentation for initialize()
261  bool initialized_;
262 
263  // stateStorageInitialized_ specifies that the state storage has been initialized.
264  // This initialization may be postponed if the linear problem was generated without
265  // the right-hand side or solution vectors.
266  bool stateStorageInitialized_;
267 
268  // Current number of iterations performed.
269  int iter_;
270 
271  // Assert that the matrix is positive definite
272  bool assertPositiveDefiniteness_;
273 
274  // Tridiagonal system for condition estimation (if needed)
275  Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
276  int numEntriesForCondEst_;
277  bool doCondEst_;
278 
279 
280 
281  //
282  // State Storage
283  //
284  // Residual
285  Teuchos::RCP<MV> R_;
286  //
287  // Preconditioned residual
288  Teuchos::RCP<MV> Z_;
289  //
290  // Direction vector
291  Teuchos::RCP<MV> P_;
292  //
293  // Operator applied to direction vector
294  Teuchos::RCP<MV> AP_;
295 
296 };
297 
299  // Constructor.
300  template<class ScalarType, class MV, class OP>
302  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
303  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
304  Teuchos::ParameterList &params ):
305  lp_(problem),
306  om_(printer),
307  stest_(tester),
308  initialized_(false),
309  stateStorageInitialized_(false),
310  iter_(0),
311  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
312  numEntriesForCondEst_(params.get("Max Size For Condest",0) )
313  {
314  }
315 
317  // Setup the state storage.
318  template <class ScalarType, class MV, class OP>
320  {
321  if (!stateStorageInitialized_) {
322 
323  // Check if there is any multivector to clone from.
324  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
325  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
326  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
327  stateStorageInitialized_ = false;
328  return;
329  }
330  else {
331 
332  // Initialize the state storage
333  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
334  if (R_ == Teuchos::null) {
335  // Get the multivector that is not null.
336  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
337  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
338  "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
339  R_ = MVT::Clone( *tmp, 1 );
340  Z_ = MVT::Clone( *tmp, 1 );
341  P_ = MVT::Clone( *tmp, 1 );
342  AP_ = MVT::Clone( *tmp, 1 );
343  }
344 
345  // Tracking information for condition number estimation
346  if(numEntriesForCondEst_ > 0) {
347  diag_.resize(numEntriesForCondEst_);
348  offdiag_.resize(numEntriesForCondEst_-1);
349  }
350 
351  // State storage has now been initialized.
352  stateStorageInitialized_ = true;
353  }
354  }
355  }
356 
357 
359  // Initialize this iteration object
360  template <class ScalarType, class MV, class OP>
362  {
363  // Initialize the state storage if it isn't already.
364  if (!stateStorageInitialized_)
365  setStateSize();
366 
367  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
368  "Belos::CGIter::initialize(): Cannot initialize state storage!");
369 
370  // NOTE: In CGIter R_, the initial residual, is required!!!
371  //
372  std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
373 
374  // Create convenience variables for zero and one.
375  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
376  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
377 
378  if (newstate.R != Teuchos::null) {
379 
380  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
381  std::invalid_argument, errstr );
382  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
383  std::invalid_argument, errstr );
384 
385  // Copy basis vectors from newstate into V
386  if (newstate.R != R_) {
387  // copy over the initial residual (unpreconditioned).
388  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
389  }
390 
391  // Compute initial direction vectors
392  // Initially, they are set to the preconditioned residuals
393  //
394  if ( lp_->getLeftPrec() != Teuchos::null ) {
395  lp_->applyLeftPrec( *R_, *Z_ );
396  if ( lp_->getRightPrec() != Teuchos::null ) {
397  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
398  lp_->applyRightPrec( *Z_, *tmp );
399  Z_ = tmp;
400  }
401  }
402  else if ( lp_->getRightPrec() != Teuchos::null ) {
403  lp_->applyRightPrec( *R_, *Z_ );
404  }
405  else {
406  Z_ = R_;
407  }
408  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
409  }
410  else {
411 
412  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
413  "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
414  }
415 
416  // The solver is initialized
417  initialized_ = true;
418  }
419 
420 
422  // Iterate until the status test informs us we should stop.
423  template <class ScalarType, class MV, class OP>
425  {
426  //
427  // Allocate/initialize data structures
428  //
429  if (initialized_ == false) {
430  initialize();
431  }
432 
433  // Allocate memory for scalars.
434  Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
435  Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
436  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
437 
438  // Create convenience variables for zero and one.
439  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
441 
442  // Get the current solution vector.
443  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
444 
445  // Check that the current solution vector only has one column.
446  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
447  "Belos::CGIter::iterate(): current linear system has more than one vector!" );
448 
449  // Compute first <r,z> a.k.a. rHz
450  MVT::MvTransMv( one, *R_, *Z_, rHz );
451 
453  // Iterate until the status test tells us to stop.
454  //
455  while (stest_->checkStatus(this) != Passed) {
456 
457  // Increment the iteration
458  iter_++;
459 
460  // Multiply the current direction vector by A and store in AP_
461  lp_->applyOp( *P_, *AP_ );
462 
463  // Compute alpha := <R_,Z_> / <P_,AP_>
464  MVT::MvTransMv( one, *P_, *AP_, pAp );
465  alpha(0,0) = rHz(0,0) / pAp(0,0);
466 
467  // Check that alpha is a positive number!
468  if(assertPositiveDefiniteness_)
469  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
470  "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
471  //
472  // Update the solution vector x := x + alpha * P_
473  //
474  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
475  lp_->updateSolution();
476  //
477  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
478  //
479  rHz_old(0,0) = rHz(0,0);
480  //
481  // Compute the new residual R_ := R_ - alpha * AP_
482  //
483  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
484  //
485  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
486  // and the new direction vector p.
487  //
488  if ( lp_->getLeftPrec() != Teuchos::null ) {
489  lp_->applyLeftPrec( *R_, *Z_ );
490  if ( lp_->getRightPrec() != Teuchos::null ) {
491  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
492  lp_->applyRightPrec( *Z_, *tmp );
493  Z_ = tmp;
494  }
495  }
496  else if ( lp_->getRightPrec() != Teuchos::null ) {
497  lp_->applyRightPrec( *R_, *Z_ );
498  }
499  else {
500  Z_ = R_;
501  }
502  //
503  MVT::MvTransMv( one, *R_, *Z_, rHz );
504  //
505  beta(0,0) = rHz(0,0) / rHz_old(0,0);
506  //
507  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
508 
509  } // end while (sTest_->checkStatus(this) != Passed)
510  }
511 
512 } // end Belos namespace
513 
514 #endif /* BELOS_CG_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::ScalarTraits< ScalarType > SCT
Definition: BelosCGIter.hpp:87
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.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
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.
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
Definition: BelosCGIter.hpp:85
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIter(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)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options...
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.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
Definition: BelosCGIter.hpp:86
Traits class which defines basic operations on multivectors.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
int getNumIters() const
Get the current iteration count.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
virtual ~CGIter()
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 > *norms) const
Get the norms of the residuals native to the solver.
SCT::magnitudeType MagnitudeType
Definition: BelosCGIter.hpp:88
void resetNumIters(int iter=0)
Reset the iteration count.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > P
The current decent direction vector.
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.

Generated for Belos by doxygen 1.8.14