Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockStochasticCGIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
11 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosMatOrthoManager.hpp"
23 #include "BelosOutputManager.hpp"
24 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
31 #include "Teuchos_ScalarTraits.hpp"
33 #include "Teuchos_TimeMonitor.hpp"
34 
50 namespace Belos {
51 
52  template<class ScalarType, class MV, class OP>
53  class PseudoBlockStochasticCGIter : virtual public StochasticCGIteration<ScalarType,MV,OP> {
54 
55  public:
56 
57  //
58  // Convenience typedefs
59  //
64 
66 
67 
74  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
76  Teuchos::ParameterList &params );
77 
81 
82 
84 
85 
99  void iterate();
100 
122 
126  void initialize()
127  {
129  initializeCG(empty);
130  }
131 
141  state.R = R_;
142  state.P = P_;
143  state.AP = AP_;
144  state.Z = Z_;
145  state.Y = Y_;
146  return state;
147  }
148 
150 
151 
153 
154 
156  int getNumIters() const { return iter_; }
157 
159  void resetNumIters( int iter = 0 ) { iter_ = iter; }
160 
163  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
164 
166 
168  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
169 
171  Teuchos::RCP<MV> getStochasticVector() const { return Y_; }
172 
174 
176 
177 
179  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
180 
182  int getBlockSize() const { return 1; }
183 
185  void setBlockSize(int blockSize) {
186  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
187  "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
188  }
189 
191  bool isInitialized() { return initialized_; }
192 
194 
195  private:
196 
199  // Do all of the calculations with doubles, because that is what the Odeh and Evans 1974 constants are for.
200  // Then cast to ScalarType.
201 
202  const double p0 = -0.322232431088;
203  const double p1 = -1.0;
204  const double p2 = -0.342242088547;
205  const double p3 = -0.204231210245e-1;
206  const double p4 = -0.453642210148e-4;
207  const double q0 = 0.993484626060e-1;
208  const double q1 = 0.588581570495;
209  const double q2 = 0.531103462366;
210  const double q3 = 0.103537752850;
211  const double q4 = 0.38560700634e-2;
212  double r,p,q,y,z;
213 
214  // Return a vector with random entries that are synchronized across processors.
215  Teuchos::randomSyncedMatrix( randvec_ );
216 
217  for (int i=0; i<numRHS_; i++)
218  {
219  // Get a random number (-1,1) and rescale to (0,1).
220  r=0.5*SCT::real(randvec_[i]) + 1.0;
221 
222  // Odeh and Evans algorithm (as modified by Park & Geyer)
223  if(r < 0.5) y=std::sqrt(-2.0 * log(r));
224  else y=std::sqrt(-2.0 * log(1.0 - r));
225 
226  p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
227  q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
228 
229  if(r < 0.5) z = (p / q) - y;
230  else z = y - (p / q);
231 
232  randvec_[i] = Teuchos::as<ScalarType,double>(z);
233  }
234 
235  return randvec_;
236  }
237 
238  //
239  // Classes inputed through constructor that define the linear problem to be solved.
240  //
244 
245  //
246  // Algorithmic parameters
247  //
248  // numRHS_ is the current number of linear systems being solved.
249  int numRHS_;
250 
251  //
252  // Current solver state
253  //
254  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
255  // is capable of running; _initialize is controlled by the initialize() member method
256  // For the implications of the state of initialized_, please see documentation for initialize()
257  bool initialized_;
258 
259  // Current number of iterations performed.
260  int iter_;
261 
262  // Current number of iterations performed.
263  bool assertPositiveDefiniteness_;
264 
265  //
266  // State Storage
267  //
268  // Residual
269  Teuchos::RCP<MV> R_;
270  //
271  // Preconditioned residual
272  Teuchos::RCP<MV> Z_;
273  //
274  // Direction vector
275  Teuchos::RCP<MV> P_;
276  //
277  // Operator applied to direction vector
278  Teuchos::RCP<MV> AP_;
279  //
280  // Stochastic recurrence vector
281  Teuchos::RCP<MV> Y_;
282  //
283  // Stochastic variable storage (for normal() method)
285 
286  };
287 
289  // Constructor.
290  template<class ScalarType, class MV, class OP>
292  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
294  Teuchos::ParameterList &params ):
295  lp_(problem),
296  om_(printer),
297  stest_(tester),
298  numRHS_(0),
299  initialized_(false),
300  iter_(0),
301  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
302  {
303  }
304 
305 
307  // Initialize this iteration object
308  template <class ScalarType, class MV, class OP>
310  {
311  // Check if there is any multivector to clone from.
312  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
313  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
314  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
315  "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
316 
317  // Get the multivector that is not null.
318  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
319 
320  // Get the number of right-hand sides we're solving for now.
321  int numRHS = MVT::GetNumberVecs(*tmp);
322  numRHS_ = numRHS;
323 
324  // Initialize the state storage
325  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
326  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
327  R_ = MVT::Clone( *tmp, numRHS_ );
328  Z_ = MVT::Clone( *tmp, numRHS_ );
329  P_ = MVT::Clone( *tmp, numRHS_ );
330  AP_ = MVT::Clone( *tmp, numRHS_ );
331  Y_ = MVT::Clone( *tmp, numRHS_ );
332  }
333 
334  // Initialize the random vector container with zeros.
335  randvec_.size( numRHS_ );
336 
337  // NOTE: In StochasticCGIter R_, the initial residual, is required!!!
338  //
339  std::string errstr("Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
340 
341  // Create convenience variables for zero and one.
342  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
344 
345  if (!Teuchos::is_null(newstate.R)) {
346 
347  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
348  std::invalid_argument, errstr );
349  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
350  std::invalid_argument, errstr );
351 
352  // Copy basis vectors from newstate into V
353  if (newstate.R != R_) {
354  // copy over the initial residual (unpreconditioned).
355  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
356  }
357 
358  // Compute initial direction vectors
359  // Initially, they are set to the preconditioned residuals
360 
361  if ( lp_->getLeftPrec() != Teuchos::null ) {
362  lp_->applyLeftPrec( *R_, *Z_ );
363  if ( lp_->getRightPrec() != Teuchos::null ) {
364  Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
365  lp_->applyRightPrec( *Z_, *tmp2 );
366  Z_ = tmp2;
367  }
368  }
369  else if ( lp_->getRightPrec() != Teuchos::null ) {
370  lp_->applyRightPrec( *R_, *Z_ );
371  }
372  else {
373  Z_ = R_;
374  }
375  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
376  }
377  else {
378 
379  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
380  "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
381  }
382 
383  // The solver is initialized
384  initialized_ = true;
385  }
386 
387 
389  // Iterate until the status test informs us we should stop.
390  template <class ScalarType, class MV, class OP>
392  {
393  //
394  // Allocate/initialize data structures
395  //
396  if (initialized_ == false) {
397  initialize();
398  }
399 
400  // Allocate memory for scalars.
401  int i=0;
402  std::vector<int> index(1);
403  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
404  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(numRHS_,numRHS_);
405 
406  // Create convenience variables for zero and one.
407  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
409 
410  // Get the current solution std::vector.
411  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
412 
413  // Compute first <r,z> a.k.a. rHz
414  MVT::MvDot( *R_, *Z_, rHz );
415 
416  if ( assertPositiveDefiniteness_ )
417  for (i=0; i<numRHS_; ++i)
418  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
420  "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
421 
423  // Iterate until the status test tells us to stop.
424  //
425  while (stest_->checkStatus(this) != Passed) {
426 
427  // Increment the iteration
428  iter_++;
429 
430  // Multiply the current direction std::vector by A and store in AP_
431  lp_->applyOp( *P_, *AP_ );
432 
433  // Compute alpha := <R_,Z_> / <P_,AP_>
434  MVT::MvDot( *P_, *AP_, pAp );
435 
437 
438  for (i=0; i<numRHS_; ++i) {
439  if ( assertPositiveDefiniteness_ )
440  // Check that pAp[i] is a positive number!
441  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
443  "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
444 
445  alpha(i,i) = rHz[i] / pAp[i];
446 
447  // Compute the scaling parameter for the stochastic vector
448  zeta(i,i) = z[i] / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
449  }
450 
451  //
452  // Update the solution std::vector x := x + alpha * P_
453  //
454  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
455  lp_->updateSolution();
456 
457  // Updates the stochastic vector y := y + zeta * P_
458  MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
459 
460  //
461  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
462  //
463  for (i=0; i<numRHS_; ++i) {
464  rHz_old[i] = rHz[i];
465  }
466  //
467  // Compute the new residual R_ := R_ - alpha * AP_
468  //
469  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
470  //
471  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
472  // and the new direction std::vector p.
473  //
474  if ( lp_->getLeftPrec() != Teuchos::null ) {
475  lp_->applyLeftPrec( *R_, *Z_ );
476  if ( lp_->getRightPrec() != Teuchos::null ) {
477  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
478  lp_->applyRightPrec( *Z_, *tmp );
479  Z_ = tmp;
480  }
481  }
482  else if ( lp_->getRightPrec() != Teuchos::null ) {
483  lp_->applyRightPrec( *R_, *Z_ );
484  }
485  else {
486  Z_ = R_;
487  }
488  //
489  MVT::MvDot( *R_, *Z_, rHz );
490  if ( assertPositiveDefiniteness_ )
491  for (i=0; i<numRHS_; ++i)
492  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
494  "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
495  //
496  // Update the search directions.
497  for (i=0; i<numRHS_; ++i) {
498  beta(i,i) = rHz[i] / rHz_old[i];
499  index[0] = i;
500  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
501  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
502  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
503  }
504  //
505  } // end while (sTest_->checkStatus(this) != Passed)
506  }
507 
508 } // end Belos namespace
509 
510 #endif /* BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP */
Pure virtual base class which augments the basic interface for a stochastic conjugate gradient linear...
Collection of types and exceptions used within the Belos solvers.
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 initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
void setBlockSize(int blockSize)
Set the blocksize.
bool is_null(const std::shared_ptr< T > &p)
Pure virtual base class for defining the status testing capabilities of Belos.
static T squareroot(T x)
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(ScalarTypea)
Declaration of basic traits for the multivector type.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Traits class which defines basic operations on multivectors.
CGPositiveDefiniteFailure is thrown when the the CG &#39;alpha = p^H*A*P&#39; value is less than zero...
void resetNumIters(int iter=0)
Reset the iteration count.
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
OperatorTraits< ScalarType, MV, OP > OPT
void iterate()
This method performs stochastic CG iterations on each linear system until the status test indicates t...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
A linear system to solve, and its associated information.
bool isInitialized()
States whether the solver has been initialized or not.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
PseudoBlockStochasticCGIter(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)
PseudoBlockStochasticCGIter constructor with linear problem, solver utilities, and parameter list of ...
void initializeCG(StochasticCGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > P
The current decent direction vector.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
This class implements the stochastic pseudo-block CG iteration, where the basic stochastic CG algorit...
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 > Y
The current stochastic recurrence vector.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Structure to contain pointers to CGIteration state variables.

Generated on Fri Nov 22 2024 09:23:06 for Belos by doxygen 1.8.5