Belos  Version of the Day
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"
61 #include "Teuchos_SerialDenseMatrix.hpp"
62 #include "Teuchos_SerialDenseVector.hpp"
63 #include "Teuchos_ScalarTraits.hpp"
64 #include "Teuchos_ParameterList.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  //
90  typedef Teuchos::ScalarTraits<ScalarType> SCT;
91  typedef typename SCT::magnitudeType MagnitudeType;
92 
94 
95 
101  PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
102  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
103  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
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 
195  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
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){doCondEst_=val;}
221 
223  Teuchos::ArrayView<MagnitudeType> getDiag() {
224  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
225  // getDiag() didn't actually throw for me in that case, but why
226  // not be cautious?
227  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
228  if (static_cast<size_type> (iter_) >= diag_.size ()) {
229  return diag_ ();
230  } else {
231  return diag_ (0, iter_);
232  }
233  }
234 
236  Teuchos::ArrayView<MagnitudeType> getOffDiag() {
237  // NOTE (mfh 30 Jul 2015) The implementation as I found it
238  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
239  // debug mode) when the maximum number of iterations has been
240  // reached, because iter_ == offdiag_.size() in that case. The
241  // new logic fixes this.
242  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
243  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
244  return offdiag_ ();
245  } else {
246  return offdiag_ (0, iter_);
247  }
248  }
249 
250  private:
251 
252  //
253  // Classes inputed through constructor that define the linear problem to be solved.
254  //
255  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
256  const Teuchos::RCP<OutputManager<ScalarType> > om_;
257  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
258 
259  //
260  // Algorithmic parameters
261  //
262  // numRHS_ is the current number of linear systems being solved.
263  int numRHS_;
264 
265  //
266  // Current solver state
267  //
268  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
269  // is capable of running; _initialize is controlled by the initialize() member method
270  // For the implications of the state of initialized_, please see documentation for initialize()
271  bool initialized_;
272 
273  // Current number of iterations performed.
274  int iter_;
275 
276  // Assert that the matrix is positive definite
277  bool assertPositiveDefiniteness_;
278 
279  // Tridiagonal system for condition estimation (if needed)
280  Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
281  int numEntriesForCondEst_;
282  bool doCondEst_;
283 
284  //
285  // State Storage
286  //
287  // Residual
288  Teuchos::RCP<MV> R_;
289  //
290  // Preconditioned residual
291  Teuchos::RCP<MV> Z_;
292  //
293  // Direction vector
294  Teuchos::RCP<MV> P_;
295  //
296  // Operator applied to direction vector
297  Teuchos::RCP<MV> AP_;
298 
299  };
300 
302  // Constructor.
303  template<class ScalarType, class MV, class OP>
305  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
306  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
307  Teuchos::ParameterList &params ):
308  lp_(problem),
309  om_(printer),
310  stest_(tester),
311  numRHS_(0),
312  initialized_(false),
313  iter_(0),
314  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
315  numEntriesForCondEst_(params.get("Max Size For Condest",0) )
316  {
317  }
318 
319 
321  // Initialize this iteration object
322  template <class ScalarType, class MV, class OP>
324  {
325  // Check if there is any mltivector to clone from.
326  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
327  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
328  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
329  "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
330 
331  // Get the multivector that is not null.
332  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
333 
334  // Get the number of right-hand sides we're solving for now.
335  int numRHS = MVT::GetNumberVecs(*tmp);
336  numRHS_ = numRHS;
337 
338  // Initialize the state storage
339  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
340  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
341  R_ = MVT::Clone( *tmp, numRHS_ );
342  Z_ = MVT::Clone( *tmp, numRHS_ );
343  P_ = MVT::Clone( *tmp, numRHS_ );
344  AP_ = MVT::Clone( *tmp, numRHS_ );
345  }
346 
347  // Tracking information for condition number estimation
348  if(numEntriesForCondEst_ > 0) {
349  diag_.resize(numEntriesForCondEst_);
350  offdiag_.resize(numEntriesForCondEst_-1);
351  }
352 
353  // NOTE: In CGIter R_, the initial residual, is required!!!
354  //
355  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
356 
357  // Create convenience variables for zero and one.
358  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
359  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
360 
361  if (!Teuchos::is_null(newstate.R)) {
362 
363  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
364  std::invalid_argument, errstr );
365  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
366  std::invalid_argument, errstr );
367 
368  // Copy basis vectors from newstate into V
369  if (newstate.R != R_) {
370  // copy over the initial residual (unpreconditioned).
371  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
372  }
373 
374  // Compute initial direction vectors
375  // Initially, they are set to the preconditioned residuals
376  //
377  if ( lp_->getLeftPrec() != Teuchos::null ) {
378  lp_->applyLeftPrec( *R_, *Z_ );
379  if ( lp_->getRightPrec() != Teuchos::null ) {
380  Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
381  lp_->applyRightPrec( *Z_, *tmp1 );
382  Z_ = tmp1;
383  }
384  }
385  else if ( lp_->getRightPrec() != Teuchos::null ) {
386  lp_->applyRightPrec( *R_, *Z_ );
387  }
388  else {
389  Z_ = R_;
390  }
391  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
392  }
393  else {
394 
395  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
396  "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
397  }
398 
399  // The solver is initialized
400  initialized_ = true;
401  }
402 
403 
405  // Iterate until the status test informs us we should stop.
406  template <class ScalarType, class MV, class OP>
408  {
409  //
410  // Allocate/initialize data structures
411  //
412  if (initialized_ == false) {
413  initialize();
414  }
415 
416  // Allocate memory for scalars.
417  int i=0;
418  std::vector<int> index(1);
419  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
420  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
421 
422  // Create convenience variables for zero and one.
423  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
424  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
425 
426  // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
427  ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
428 
429  // Get the current solution std::vector.
430  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
431 
432  // Compute first <r,z> a.k.a. rHz
433  MVT::MvDot( *R_, *Z_, rHz );
434 
435  if ( assertPositiveDefiniteness_ )
436  for (i=0; i<numRHS_; ++i)
437  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
439  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
440 
442  // Iterate until the status test tells us to stop.
443  //
444  while (stest_->checkStatus(this) != Passed) {
445 
446  // Increment the iteration
447  iter_++;
448 
449  // Multiply the current direction std::vector by A and store in AP_
450  lp_->applyOp( *P_, *AP_ );
451 
452  // Compute alpha := <R_,Z_> / <P_,AP_>
453  MVT::MvDot( *P_, *AP_, pAp );
454 
455  for (i=0; i<numRHS_; ++i) {
456  if ( assertPositiveDefiniteness_ )
457  // Check that pAp[i] is a positive number!
458  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
460  "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
461 
462  alpha(i,i) = rHz[i] / pAp[i];
463  }
464 
465  //
466  // Update the solution std::vector x := x + alpha * P_
467  //
468  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
469  lp_->updateSolution();
470  //
471  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
472  //
473  for (i=0; i<numRHS_; ++i) {
474  rHz_old[i] = rHz[i];
475  }
476  //
477  // Compute the new residual R_ := R_ - alpha * AP_
478  //
479  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
480  //
481  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
482  // and the new direction std::vector p.
483  //
484  if ( lp_->getLeftPrec() != Teuchos::null ) {
485  lp_->applyLeftPrec( *R_, *Z_ );
486  if ( lp_->getRightPrec() != Teuchos::null ) {
487  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
488  lp_->applyRightPrec( *Z_, *tmp );
489  Z_ = tmp;
490  }
491  }
492  else if ( lp_->getRightPrec() != Teuchos::null ) {
493  lp_->applyRightPrec( *R_, *Z_ );
494  }
495  else {
496  Z_ = R_;
497  }
498  //
499  MVT::MvDot( *R_, *Z_, rHz );
500  if ( assertPositiveDefiniteness_ )
501  for (i=0; i<numRHS_; ++i)
502  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
504  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
505  //
506  // Update the search directions.
507  for (i=0; i<numRHS_; ++i) {
508  beta(i,i) = rHz[i] / rHz_old[i];
509  index[0] = i;
510  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
511  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
512  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
513  }
514 
515  // Condition estimate (if needed)
516  if(doCondEst_ > 0) {
517  if(iter_ > 1 ) {
518  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
519  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
520  }
521  else {
522  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
523  }
524  rHz_old2 = rHz_old[0];
525  beta_old = beta(0,0);
526  pAp_old = pAp[0];
527  }
528 
529 
530  //
531  } // end while (sTest_->checkStatus(this) != Passed)
532  }
533 
534 } // end Belos namespace
535 
536 #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.
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.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
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.
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.
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 getNumIters() const
Get the current iteration count.
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.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
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...
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated for Belos by doxygen 1.8.14