Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockCGIter.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_PSEUDO_BLOCK_CG_ITER_HPP
43 #define BELOS_PSEUDO_BLOCK_CG_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 
80  template<class ScalarType, class MV, class OP>
81  class PseudoBlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
82 
83  public:
84 
85  //
86  // Convenience typedefs
87  //
92 
94 
95 
102  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
104  Teuchos::ParameterList &params );
105 
107  virtual ~PseudoBlockCGIter() {};
109 
110 
112 
113 
127  void iterate();
128 
150 
154  void initialize()
155  {
157  initializeCG(empty);
158  }
159 
169  state.R = R_;
170  state.P = P_;
171  state.AP = AP_;
172  state.Z = Z_;
173  return state;
174  }
175 
177 
178 
180 
181 
183  int getNumIters() const { return iter_; }
184 
186  void resetNumIters( int iter = 0 ) { iter_ = iter; }
187 
190  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
191 
193 
196 
198 
200 
201 
203  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
204 
206  int getBlockSize() const { return 1; }
207 
209  void setBlockSize(int blockSize) {
210  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
211  "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
212  }
213 
215  bool isInitialized() { return initialized_; }
216 
218 
220  void setDoCondEst(bool val) {
222  }
223 
226  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
227  // getDiag() didn't actually throw for me in that case, but why
228  // not be cautious?
229  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
230  if (static_cast<size_type> (iter_) >= diag_.size ()) {
231  return diag_ ();
232  } else {
233  return diag_ (0, iter_);
234  }
235  }
236 
239  // NOTE (mfh 30 Jul 2015) The implementation as I found it
240  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
241  // debug mode) when the maximum number of iterations has been
242  // reached, because iter_ == offdiag_.size() in that case. The
243  // new logic fixes this.
244  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
245  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
246  return offdiag_ ();
247  } else {
248  return offdiag_ (0, iter_);
249  }
250  }
251 
252  private:
253 
254  //
255  // Classes inputed through constructor that define the linear problem to be solved.
256  //
260 
261  //
262  // Algorithmic parameters
263  //
264  // numRHS_ is the current number of linear systems being solved.
265  int numRHS_;
266 
267  //
268  // Current solver state
269  //
270  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
271  // is capable of running; _initialize is controlled by the initialize() member method
272  // For the implications of the state of initialized_, please see documentation for initialize()
274 
275  // Current number of iterations performed.
276  int iter_;
277 
278  // Assert that the matrix is positive definite
280 
281  // Tridiagonal system for condition estimation (if needed)
283  ScalarType pAp_old_, beta_old_, rHz_old2_; // Put scalars here so that estimate is correct for multiple RHS, when deflation occurs.
286 
287  //
288  // State Storage
289  //
290  // Residual
292  //
293  // Preconditioned residual
295  //
296  // Direction vector
298  //
299  // Operator applied to direction vector
301 
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  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
318  numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
319  doCondEst_(false)
320  {
321  }
322 
323 
325  // Initialize this iteration object
326  template <class ScalarType, class MV, class OP>
328  {
329  // Check if there is any mltivector to clone from.
330  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
331  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
332  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
333  "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
334 
335  // Get the multivector that is not null.
336  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
337 
338  // Get the number of right-hand sides we're solving for now.
339  int numRHS = MVT::GetNumberVecs(*tmp);
340  numRHS_ = numRHS;
341 
342  // Initialize the state storage
343  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
344  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
345  R_ = MVT::Clone( *tmp, numRHS_ );
346  Z_ = MVT::Clone( *tmp, numRHS_ );
347  P_ = MVT::Clone( *tmp, numRHS_ );
348  AP_ = MVT::Clone( *tmp, numRHS_ );
349  }
350 
351  // Tracking information for condition number estimation
352  if(numEntriesForCondEst_ > 0) {
353  diag_.resize(numEntriesForCondEst_);
354  offdiag_.resize(numEntriesForCondEst_-1);
355  }
356 
357  // NOTE: In CGIter R_, the initial residual, is required!!!
358  //
359  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
360 
361  // Create convenience variables for zero and one.
362  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
364 
365  if (!Teuchos::is_null(newstate.R)) {
366 
367  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
368  std::invalid_argument, errstr );
369  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
370  std::invalid_argument, errstr );
371 
372  // Copy basis vectors from newstate into V
373  if (newstate.R != R_) {
374  // copy over the initial residual (unpreconditioned).
375  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
376  }
377 
378  // Compute initial direction vectors
379  // Initially, they are set to the preconditioned residuals
380  //
381  if ( lp_->getLeftPrec() != Teuchos::null ) {
382  lp_->applyLeftPrec( *R_, *Z_ );
383  if ( lp_->getRightPrec() != Teuchos::null ) {
384  Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
385  lp_->applyRightPrec( *Z_, *tmp1 );
386  Z_ = tmp1;
387  }
388  }
389  else if ( lp_->getRightPrec() != Teuchos::null ) {
390  lp_->applyRightPrec( *R_, *Z_ );
391  }
392  else {
393  Z_ = R_;
394  }
395  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
396  }
397  else {
398 
399  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
400  "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
401  }
402 
403  // The solver is initialized
404  initialized_ = true;
405  }
406 
407 
409  // Iterate until the status test informs us we should stop.
410  template <class ScalarType, class MV, class OP>
412  {
413  //
414  // Allocate/initialize data structures
415  //
416  if (initialized_ == false) {
417  initialize();
418  }
419 
420  // Allocate memory for scalars.
421  int i=0;
422  std::vector<int> index(1);
423  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ ), beta( numRHS_ );
424  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ );
425 
426  // Create convenience variables for zero and one.
427  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
429 
430  // Get the current solution std::vector.
431  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
432 
433  // Compute first <r,z> a.k.a. rHz
434  MVT::MvDot( *R_, *Z_, rHz );
435 
436  if ( assertPositiveDefiniteness_ )
437  for (i=0; i<numRHS_; ++i)
438  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
440  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
441 
443  // Iterate until the status test tells us to stop.
444  //
445  while (stest_->checkStatus(this) != Passed) {
446 
447  // Increment the iteration
448  iter_++;
449 
450  // Multiply the current direction std::vector by A and store in AP_
451  lp_->applyOp( *P_, *AP_ );
452 
453  // Compute alpha := <R_,Z_> / <P_,AP_>
454  MVT::MvDot( *P_, *AP_, pAp );
455 
456  for (i=0; i<numRHS_; ++i) {
457  if ( assertPositiveDefiniteness_ )
458  // Check that pAp[i] is a positive number!
459  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
461  "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
462 
463  alpha(i,i) = rHz[i] / pAp[i];
464  }
465 
466  //
467  // Update the solution std::vector x := x + alpha * P_
468  //
469  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
470  lp_->updateSolution();// what does this do?
471  //
472  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
473  //
474  for (i=0; i<numRHS_; ++i) {
475  rHz_old[i] = rHz[i];
476  }
477  //
478  // Compute the new residual R_ := R_ - alpha * AP_
479  //
480  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
481  //
482  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
483  // and the new direction std::vector p.
484  //
485  if ( lp_->getLeftPrec() != Teuchos::null ) {
486  lp_->applyLeftPrec( *R_, *Z_ );
487  if ( lp_->getRightPrec() != Teuchos::null ) {
488  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
489  lp_->applyRightPrec( *Z_, *tmp );
490  Z_ = tmp;
491  }
492  }
493  else if ( lp_->getRightPrec() != Teuchos::null ) {
494  lp_->applyRightPrec( *R_, *Z_ );
495  }
496  else {
497  Z_ = R_;
498  }
499  //
500  MVT::MvDot( *R_, *Z_, rHz );
501  if ( assertPositiveDefiniteness_ )
502  for (i=0; i<numRHS_; ++i)
503  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
505  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
506  //
507  // Update the search directions.
508  for (i=0; i<numRHS_; ++i) {
509  beta[i] = rHz[i] / rHz_old[i];
510  index[0] = i;
511  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
512  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
513  MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
514  }
515 
516  // Condition estimate (if needed)
517  if (doCondEst_) {
518  if (iter_ > 1) {
519  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old_ * beta_old_ * pAp_old_ + pAp[0]) / rHz_old[0]);
520  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old_ * pAp_old_ / (sqrt( rHz_old[0] * rHz_old2_)));
521  }
522  else {
523  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
524  }
525  rHz_old2_ = rHz_old[0];
526  beta_old_ = beta[0];
527  pAp_old_ = pAp[0];
528  }
529 
530 
531  //
532  } // end while (sTest_->checkStatus(this) != Passed)
533  }
534 
535 } // end Belos namespace
536 
537 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
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.
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.
virtual ~PseudoBlockCGIter()
Destructor.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Structure to contain pointers to CGIteration state variables.
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)
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::ArrayRCP< MagnitudeType > offdiag_
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize.
size_type size() const
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.
OperatorTraits< ScalarType, MV, OP > OPT
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void resetNumIters(int iter=0)
Reset the iteration count.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
A linear system to solve, and its associated information.
int getNumIters() const
Get the current iteration count.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::ScalarTraits< ScalarType > SCT
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
PseudoBlockCGIter(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)
PseudoBlockCGIter constructor with linear problem, solver utilities, and parameter list of solver opt...
TransListIter iter
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Teuchos::RCP< OutputManager< ScalarType > > om_
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ArrayRCP< MagnitudeType > diag_