Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Belos_PseudoBlockCGIter_MP_Vector.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_MP_VECTOR_HPP
43 #define BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
44 
54 #include "BelosPseudoBlockCGIter.hpp"
55 
71 namespace Belos {
72 
73  template<class Storage, class MV, class OP>
74  class PseudoBlockCGIter<Sacado::MP::Vector<Storage>, MV, OP> :
75  virtual public CGIteration<Sacado::MP::Vector<Storage>, MV, OP> {
76 
77  public:
78 
79  //
80  // Convenience typedefs
81  //
83  typedef MultiVecTraits<ScalarType,MV> MVT;
84  typedef OperatorTraits<ScalarType,MV,OP> OPT;
85  typedef Teuchos::ScalarTraits<ScalarType> SCT;
86  typedef typename SCT::magnitudeType MagnitudeType;
87  typedef Teuchos::ScalarTraits<typename Storage::value_type> SVT;
88 
90 
91 
97  PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
98  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
99  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
100  Teuchos::ParameterList &params );
101 
103  virtual ~PseudoBlockCGIter() {};
105 
106 
108 
109 
123  void iterate();
124 
145  void initializeCG(CGIterationState<ScalarType,MV>& newstate);
146 
150  void initialize()
151  {
152  CGIterationState<ScalarType,MV> empty;
153  initializeCG(empty);
154  }
155 
163  CGIterationState<ScalarType,MV> getState() const {
164  CGIterationState<ScalarType,MV> state;
165  state.R = R_;
166  state.P = P_;
167  state.AP = AP_;
168  state.Z = Z_;
169  return state;
170  }
171 
173 
174 
176 
177 
179  int getNumIters() const { return iter_; }
180 
182  void resetNumIters( int iter = 0 ) { iter_ = iter; }
183 
186  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
187 
189 
191  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
192 
194 
196 
197 
199  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
200 
202  int getBlockSize() const { return 1; }
203 
205  void setBlockSize(int blockSize) {
206  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
207  "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
208  }
209 
211  bool isInitialized() { return initialized_; }
212 
214 
216  void setDoCondEst(bool val){doCondEst_=val;}
217 
219  Teuchos::ArrayView<MagnitudeType> getDiag() {
220  // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
221  // getDiag() didn't actually throw for me in that case, but why
222  // not be cautious?
223  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
224  if (static_cast<size_type> (iter_) >= diag_.size ()) {
225  return diag_ ();
226  } else {
227  return diag_ (0, iter_);
228  }
229  }
230 
232  Teuchos::ArrayView<MagnitudeType> getOffDiag() {
233  // NOTE (mfh 30 Jul 2015) The implementation as I found it
234  // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
235  // debug mode) when the maximum number of iterations has been
236  // reached, because iter_ == offdiag_.size() in that case. The
237  // new logic fixes this.
238  typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
239  if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
240  return offdiag_ ();
241  } else {
242  return offdiag_ (0, iter_);
243  }
244  }
245 
246  private:
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  // Algorithmic parameters
257  //
258  // numRHS_ is the current number of linear systems being solved.
259  int numRHS_;
260 
261  //
262  // Current solver state
263  //
264  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
265  // is capable of running; _initialize is controlled by the initialize() member method
266  // For the implications of the state of initialized_, please see documentation for initialize()
268 
269  // Current number of iterations performed.
270  int iter_;
271 
272  // Assert that the matrix is positive definite
274 
275  // Tridiagonal system for condition estimation (if needed)
276  Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
279 
280  //
281  // State Storage
282  //
283  // Residual
284  Teuchos::RCP<MV> R_;
285  //
286  // Preconditioned residual
287  Teuchos::RCP<MV> Z_;
288  //
289  // Direction vector
290  Teuchos::RCP<MV> P_;
291  //
292  // Operator applied to direction vector
293  Teuchos::RCP<MV> AP_;
294 
295  // Tolerance for ensemble breakdown
296  typename SVT::magnitudeType breakDownTol_;
297 
298  };
299 
301  // Constructor.
302  template<class StorageType, class MV, class OP>
304  PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
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  breakDownTol_(params.get("Ensemble Breakdown Tolerance", 0.0))
317  {
318  }
319 
320 
322  // Initialize this iteration object
323  template <class StorageType, class MV, class OP>
324  void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
325  initializeCG(CGIterationState<ScalarType,MV>& newstate)
326  {
327  // Check if there is any mltivector to clone from.
328  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
329  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
330  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
331  "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
332 
333  // Get the multivector that is not null.
334  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
335 
336  // Get the number of right-hand sides we're solving for now.
337  int numRHS = MVT::GetNumberVecs(*tmp);
338  numRHS_ = numRHS;
339 
340  // Initialize the state storage
341  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
342  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
343  R_ = MVT::Clone( *tmp, numRHS_ );
344  Z_ = MVT::Clone( *tmp, numRHS_ );
345  P_ = MVT::Clone( *tmp, numRHS_ );
346  AP_ = MVT::Clone( *tmp, numRHS_ );
347  }
348 
349  // Tracking information for condition number estimation
350  if(numEntriesForCondEst_ > 0) {
351  diag_.resize(numEntriesForCondEst_);
352  offdiag_.resize(numEntriesForCondEst_-1);
353  }
354 
355  // NOTE: In CGIter R_, the initial residual, is required!!!
356  //
357  std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
358 
359  // Create convenience variables for zero and one.
360  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
361  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
362 
363  if (!Teuchos::is_null(newstate.R)) {
364 
365  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
366  std::invalid_argument, errstr );
367  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
368  std::invalid_argument, errstr );
369 
370  // Copy basis vectors from newstate into V
371  if (newstate.R != R_) {
372  // copy over the initial residual (unpreconditioned).
373  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
374  }
375 
376  // Compute initial direction vectors
377  // Initially, they are set to the preconditioned residuals
378  //
379  if ( lp_->getLeftPrec() != Teuchos::null ) {
380  lp_->applyLeftPrec( *R_, *Z_ );
381  if ( lp_->getRightPrec() != Teuchos::null ) {
382  Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
383  lp_->applyRightPrec( *Z_, *tmp1 );
384  Z_ = tmp1;
385  }
386  }
387  else if ( lp_->getRightPrec() != Teuchos::null ) {
388  lp_->applyRightPrec( *R_, *Z_ );
389  }
390  else {
391  Z_ = R_;
392  }
393  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
394  }
395  else {
396 
397  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
398  "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
399  }
400 
401  // The solver is initialized
402  initialized_ = true;
403  }
404 
405 
407  // Iterate until the status test informs us we should stop.
408  template <class StorageType, class MV, class OP>
409  void PseudoBlockCGIter<Sacado::MP::Vector<StorageType>,MV,OP>::
410  iterate()
411  {
412  //
413  // Allocate/initialize data structures
414  //
415  if (initialized_ == false) {
416  initialize();
417  }
418 
419  // Allocate memory for scalars.
420  int i=0;
421  std::vector<int> index(1);
422  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
423  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
424 
425  // Create convenience variables for zero and one.
426  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
427  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
428 
429  // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
430  ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
431 
432  // Get the current solution std::vector.
433  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
434 
435  // Compute first <r,z> a.k.a. rHz
436  MVT::MvDot( *R_, *Z_, rHz );
437 
438  if ( assertPositiveDefiniteness_ )
439  for (i=0; i<numRHS_; ++i)
440  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
441  CGIterateFailure,
442  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
443 
445  // Iterate until the status test tells us to stop.
446  //
447  while (stest_->checkStatus(this) != Passed) {
448 
449  // Increment the iteration
450  iter_++;
451 
452  // Multiply the current direction std::vector by A and store in AP_
453  lp_->applyOp( *P_, *AP_ );
454 
455  // Compute alpha := <R_,Z_> / <P_,AP_>
456  MVT::MvDot( *P_, *AP_, pAp );
457 
458  for (i=0; i<numRHS_; ++i) {
459  if ( assertPositiveDefiniteness_ )
460  // Check that pAp[i] is a positive number!
461  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
462  CGIterateFailure,
463  "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
464 
465  const int ensemble_size = pAp[i].size();
466  for (int k=0; k<ensemble_size; ++k) {
467  if (SVT::magnitude(pAp[i].coeff(k)) <= breakDownTol_)
468  alpha(i,i).fastAccessCoeff(k) = 0.0;
469  else
470  alpha(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / pAp[i].coeff(k);
471  }
472  }
473 
474  //
475  // Update the solution std::vector x := x + alpha * P_
476  //
477  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
478  lp_->updateSolution();
479  //
480  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
481  //
482  for (i=0; i<numRHS_; ++i) {
483  rHz_old[i] = rHz[i];
484  }
485  //
486  // Compute the new residual R_ := R_ - alpha * AP_
487  //
488  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
489  //
490  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
491  // and the new direction std::vector p.
492  //
493  if ( lp_->getLeftPrec() != Teuchos::null ) {
494  lp_->applyLeftPrec( *R_, *Z_ );
495  if ( lp_->getRightPrec() != Teuchos::null ) {
496  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
497  lp_->applyRightPrec( *Z_, *tmp );
498  Z_ = tmp;
499  }
500  }
501  else if ( lp_->getRightPrec() != Teuchos::null ) {
502  lp_->applyRightPrec( *R_, *Z_ );
503  }
504  else {
505  Z_ = R_;
506  }
507  //
508  MVT::MvDot( *R_, *Z_, rHz );
509  if ( assertPositiveDefiniteness_ )
510  for (i=0; i<numRHS_; ++i)
511  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
512  CGIterateFailure,
513  "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
514  //
515  // Update the search directions.
516  for (i=0; i<numRHS_; ++i) {
517  const int ensemble_size = rHz[i].size();
518  for (int k=0; k<ensemble_size; ++k) {
519  if (SVT::magnitude(rHz[i].coeff(k)) <= breakDownTol_)
520  beta(i,i).fastAccessCoeff(k) = 0.0;
521  else
522  beta(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / rHz_old[i].coeff(k);
523  }
524  index[0] = i;
525  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
526  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
527  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
528  }
529 
530  // Condition estimate (if needed)
531  if(doCondEst_ > 0) {
532  if(iter_ > 1 ) {
533  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
534  offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
535  }
536  else {
537  diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
538  }
539  rHz_old2 = rHz_old[0];
540  beta_old = beta(0,0);
541  pAp_old = pAp[0];
542  }
543 
544 
545  //
546  } // end while (sTest_->checkStatus(this) != Passed)
547  }
548 
549 } // end Belos namespace
550 
551 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
expr val()
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void setDoCondEst(bool val)
Sets whether or not to store 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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.