Belos Version of the Day
Loading...
Searching...
No Matches
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
55#include "BelosStatusTest.hpp"
59
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
65
76namespace Belos {
77
78template<class ScalarType, class MV, class OP>
79class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
80
81 public:
82
83 //
84 // Convenience typedefs
85 //
88 typedef Teuchos::ScalarTraits<ScalarType> SCT;
89 typedef typename SCT::magnitudeType MagnitudeType;
90
92
93
99 CGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
100 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
101 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
102 const Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > &convTester,
103 Teuchos::ParameterList &params );
104
106 virtual ~CGIter() {};
108
109
111
112
125 void iterate();
126
142
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 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
183 (*norms)[0] = std::sqrt(Teuchos::ScalarTraits<ScalarType>::magnitude(rHr_));
184 return Teuchos::null;
185 } else
186 return R_;
187 }
188
190
192 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
193
195
197
198
200 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
201
203 int getBlockSize() const { return 1; }
204
206 void setBlockSize(int blockSize) {
207 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
208 "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
209 }
210
212 bool isInitialized() { return initialized_; }
213
214
216 void setDoCondEst(bool val) {
217 if (numEntriesForCondEst_) doCondEst_=val;
218 }
219
221 Teuchos::ArrayView<MagnitudeType> getDiag() {
222 // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
223 // getDiag() didn't actually throw for me in that case, but why
224 // not be cautious?
225 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
226 if (static_cast<size_type> (iter_) >= diag_.size ()) {
227 return diag_ ();
228 } else {
229 return diag_ (0, iter_);
230 }
231 }
232
234 Teuchos::ArrayView<MagnitudeType> getOffDiag() {
235 // NOTE (mfh 30 Jul 2015) The implementation as I found it
236 // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
237 // debug mode) when the maximum number of iterations has been
238 // reached, because iter_ == offdiag_.size() in that case. The
239 // new logic fixes this.
240 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
241 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
242 return offdiag_ ();
243 } else {
244 return offdiag_ (0, iter_);
245 }
246 }
247
249
250 private:
251
252 //
253 // Internal methods
254 //
256 void setStateSize();
257
258 //
259 // Classes inputed through constructor that define the linear problem to be solved.
260 //
261 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
262 const Teuchos::RCP<OutputManager<ScalarType> > om_;
263 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
264 const Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
265
266 //
267 // Current solver state
268 //
269 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
270 // is capable of running; _initialize is controlled by the initialize() member method
271 // For the implications of the state of initialized_, please see documentation for initialize()
272 bool initialized_;
273
274 // stateStorageInitialized_ specifies that the state storage has been initialized.
275 // This initialization may be postponed if the linear problem was generated without
276 // the right-hand side or solution vectors.
277 bool stateStorageInitialized_;
278
279 // Current number of iterations performed.
280 int iter_;
281
282 // Should the allreduce for convergence detection be merged with one of the inner products?
283 bool foldConvergenceDetectionIntoAllreduce_;
284
285 // <r,r>
286 ScalarType rHr_;
287
288 // Assert that the matrix is positive definite
289 bool assertPositiveDefiniteness_;
290
291 // Tridiagonal system for condition estimation (if needed)
292 Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
293 int numEntriesForCondEst_;
294 bool doCondEst_;
295
296
297
298 //
299 // State Storage
300 //
301 // Residual
302 Teuchos::RCP<MV> R_;
303 //
304 // Preconditioned residual
305 Teuchos::RCP<MV> Z_;
306 //
307 // Direction vector
308 Teuchos::RCP<MV> P_;
309 //
310 // Operator applied to direction vector
311 Teuchos::RCP<MV> AP_;
312
313 Teuchos::RCP<MV> S_;
314
315};
316
318 // Constructor.
319 template<class ScalarType, class MV, class OP>
321 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
322 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
323 const Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > &convTester,
324 Teuchos::ParameterList &params ):
325 lp_(problem),
326 om_(printer),
327 stest_(tester),
328 convTest_(convTester),
329 initialized_(false),
330 stateStorageInitialized_(false),
331 iter_(0),
332 assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
333 numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
334 doCondEst_(false)
335 {
336 foldConvergenceDetectionIntoAllreduce_ = params.get<bool>("Fold Convergence Detection Into Allreduce",false);
337 }
338
340 // Setup the state storage.
341 template <class ScalarType, class MV, class OP>
343 {
344 if (!stateStorageInitialized_) {
345
346 // Check if there is any multivector to clone from.
347 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
348 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ = false;
351 return;
352 }
353 else {
354
355 // Initialize the state storage
356 // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
357 if (R_ == Teuchos::null) {
358 // Get the multivector that is not null.
359 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
360 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
361 "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
362 S_ = MVT::Clone( *tmp, 2 );
363 std::vector<int> index(1,0);
364 index[0] = 0;
365 R_ = MVT::CloneViewNonConst( *S_, index );
366 index[0] = 1;
367 Z_ = MVT::CloneViewNonConst( *S_, index );
368 P_ = MVT::Clone( *tmp, 1 );
369 AP_ = MVT::Clone( *tmp, 1 );
370
371 }
372
373 // Tracking information for condition number estimation
374 if(numEntriesForCondEst_ > 0) {
375 diag_.resize(numEntriesForCondEst_);
376 offdiag_.resize(numEntriesForCondEst_-1);
377 }
378
379 // State storage has now been initialized.
380 stateStorageInitialized_ = true;
381 }
382 }
383 }
384
385
387 // Initialize this iteration object
388 template <class ScalarType, class MV, class OP>
390 {
391 // Initialize the state storage if it isn't already.
392 if (!stateStorageInitialized_)
393 setStateSize();
394
395 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
396 "Belos::CGIter::initialize(): Cannot initialize state storage!");
397
398 // NOTE: In CGIter R_, the initial residual, is required!!!
399 //
400 std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
401
402 if (newstate.R != Teuchos::null) {
403
404 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
405 std::invalid_argument, errstr );
406 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
407 std::invalid_argument, errstr );
408
409 // Copy basis vectors from newstate into V
410 if (newstate.R != R_) {
411 // copy over the initial residual (unpreconditioned).
412 MVT::Assign( *newstate.R, *R_ );
413 }
414
415 // Compute initial direction vectors
416 // Initially, they are set to the preconditioned residuals
417 //
418 if ( lp_->getLeftPrec() != Teuchos::null ) {
419 lp_->applyLeftPrec( *R_, *Z_ );
420 if ( lp_->getRightPrec() != Teuchos::null ) {
421 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_ );
422 lp_->applyRightPrec( *tmp, *Z_ );
423 }
424 }
425 else if ( lp_->getRightPrec() != Teuchos::null ) {
426 lp_->applyRightPrec( *R_, *Z_ );
427 }
428 else {
429 MVT::Assign( *R_, *Z_ );
430 }
431 MVT::Assign( *Z_, *P_ );
432 }
433 else {
434
435 TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
436 "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
437 }
438
439 // The solver is initialized
440 initialized_ = true;
441 }
442
443
445 // Iterate until the status test informs us we should stop.
446 template <class ScalarType, class MV, class OP>
448 {
449 //
450 // Allocate/initialize data structures
451 //
452 if (initialized_ == false) {
453 initialize();
454 }
455
456 // Allocate memory for scalars.
457 std::vector<ScalarType> alpha(1);
458 std::vector<ScalarType> beta(1);
459 std::vector<ScalarType> rHz(1), rHz_old(1), pAp(1);
460 Teuchos::SerialDenseMatrix<int,ScalarType> rHs( 1, 2 );
461
462 // Create convenience variables for zero and one.
463 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
464 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
465
466 // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
467 ScalarType pAp_old = one, beta_old = one, rHz_old2 = one;
468
469 // Get the current solution vector.
470 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
471
472 // Check that the current solution vector only has one column.
473 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
474 "Belos::CGIter::iterate(): current linear system has more than one vector!" );
475
476 // Compute first <r,z> a.k.a. rHz
477 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
478 MVT::MvTransMv( one, *R_, *S_, rHs );
479 rHr_ = rHs(0,0);
480 rHz[0] = rHs(0,1);
481 } else
482 MVT::MvDot( *R_, *Z_, rHz );
483
485 // Iterate until the status test tells us to stop.
486 //
487 while (stest_->checkStatus(this) != Passed) {
488
489 // Increment the iteration
490 iter_++;
491
492 // Multiply the current direction vector by A and store in AP_
493 lp_->applyOp( *P_, *AP_ );
494
495 // Compute alpha := <R_,Z_> / <P_,AP_>
496 MVT::MvDot( *P_, *AP_, pAp );
497 alpha[0] = rHz[0] / pAp[0];
498
499 // Check that alpha is a positive number!
500 if(assertPositiveDefiniteness_) {
501 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha[0]) <= zero, CGPositiveDefiniteFailure,
502 "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
503 }
504
505 //
506 // Update the solution vector x := x + alpha * P_
507 //
508 MVT::MvAddMv( one, *cur_soln_vec, alpha[0], *P_, *cur_soln_vec );
509 lp_->updateSolution();
510 //
511 // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
512 //
513 rHz_old[0] = rHz[0];
514 //
515 // Compute the new residual R_ := R_ - alpha * AP_
516 //
517 MVT::MvAddMv( one, *R_, -alpha[0], *AP_, *R_ );
518 //
519 // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
520 // and the new direction vector p.
521 //
522 if ( lp_->getLeftPrec() != Teuchos::null ) {
523 lp_->applyLeftPrec( *R_, *Z_ );
524 if ( lp_->getRightPrec() != Teuchos::null ) {
525 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Z_);
526 lp_->applyRightPrec( *tmp, *Z_ );
527 }
528 }
529 else if ( lp_->getRightPrec() != Teuchos::null ) {
530 lp_->applyRightPrec( *R_, *Z_ );
531 }
532 else {
533 MVT::Assign( *R_, *Z_ );
534 }
535 //
536 if (foldConvergenceDetectionIntoAllreduce_ && convTest_->getResNormType() == Belos::TwoNorm) {
537 MVT::MvTransMv( one, *R_, *S_, rHs );
538 rHr_ = rHs(0,0);
539 rHz[0] = rHs(0,1);
540 } else
541 MVT::MvDot( *R_, *Z_, rHz );
542 //
543 beta[0] = rHz[0] / rHz_old[0];
544 //
545 MVT::MvAddMv( one, *Z_, beta[0], *P_, *P_ );
546
547 // Condition estimate (if needed)
548 if (doCondEst_) {
549 if (iter_ > 1) {
550 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
551 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
552 }
553 else {
554 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
555 }
556 rHz_old2 = rHz_old[0];
557 beta_old = beta[0];
558 pAp_old = pAp[0];
559 }
560
561 } // end while (sTest_->checkStatus(this) != Passed)
562 }
563
564} // end Belos namespace
565
566#endif /* BELOS_CG_ITER_HPP */
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:79
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
SCT::magnitudeType MagnitudeType
Definition: BelosCGIter.hpp:89
OperatorTraits< ScalarType, MV, OP > OPT
Definition: BelosCGIter.hpp:87
virtual ~CGIter()
Destructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::ScalarTraits< ScalarType > SCT
Definition: BelosCGIter.hpp:88
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
void resetNumIters(int iter=0)
Reset the iteration count.
MultiVecTraits< ScalarType, MV > MVT
Definition: BelosCGIter.hpp:86
CGIter(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)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
An implementation of StatusTestResNorm using a family of residual norms.
A pure virtual class for defining the status tests for the Belos iterative solvers.
@ TwoNorm
Definition: BelosTypes.hpp:98
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.

Generated for Belos by doxygen 1.9.6