Belos Package Browser (Single Doxygen Collection)  Development
 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 //
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_STOCHASTIC_CG_ITER_HPP
43 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.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"
64 #include "Teuchos_ScalarTraits.hpp"
66 #include "Teuchos_TimeMonitor.hpp"
67 
83 namespace Belos {
84 
85  template<class ScalarType, class MV, class OP>
86  class PseudoBlockStochasticCGIter : virtual public StochasticCGIteration<ScalarType,MV,OP> {
87 
88  public:
89 
90  //
91  // Convenience typedefs
92  //
97 
99 
100 
107  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
109  Teuchos::ParameterList &params );
110 
114 
115 
117 
118 
132  void iterate();
133 
155 
159  void initialize()
160  {
162  initializeCG(empty);
163  }
164 
174  state.R = R_;
175  state.P = P_;
176  state.AP = AP_;
177  state.Z = Z_;
178  state.Y = Y_;
179  return state;
180  }
181 
183 
184 
186 
187 
189  int getNumIters() const { return iter_; }
190 
192  void resetNumIters( int iter = 0 ) { iter_ = iter; }
193 
196  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
197 
199 
202 
205 
207 
209 
210 
212  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
213 
215  int getBlockSize() const { return 1; }
216 
218  void setBlockSize(int blockSize) {
219  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
220  "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
221  }
222 
224  bool isInitialized() { return initialized_; }
225 
227 
228  private:
229 
232  // Do all of the calculations with doubles, because that is what the Odeh and Evans 1974 constants are for.
233  // Then cast to ScalarType.
234 
235  const double p0 = -0.322232431088;
236  const double p1 = -1.0;
237  const double p2 = -0.342242088547;
238  const double p3 = -0.204231210245e-1;
239  const double p4 = -0.453642210148e-4;
240  const double q0 = 0.993484626060e-1;
241  const double q1 = 0.588581570495;
242  const double q2 = 0.531103462366;
243  const double q3 = 0.103537752850;
244  const double q4 = 0.38560700634e-2;
245  double r,p,q,y,z;
246 
247  // Return a vector with random entries that are synchronized across processors.
248  Teuchos::randomSyncedMatrix( randvec_ );
249 
250  for (int i=0; i<numRHS_; i++)
251  {
252  // Get a random number (-1,1) and rescale to (0,1).
253  r=0.5*randvec_[i] + 1.0;
254 
255  // Odeh and Evans algorithm (as modified by Park & Geyer)
256  if(r < 0.5) y=std::sqrt(-2.0 * log(r));
257  else y=std::sqrt(-2.0 * log(1.0 - r));
258 
259  p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
260  q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
261 
262  if(r < 0.5) z = (p / q) - y;
263  else z = y - (p / q);
264 
265  randvec_[i] = Teuchos::as<ScalarType,double>(z);
266  }
267 
268  return randvec_;
269  }
270 
271  //
272  // Classes inputed through constructor that define the linear problem to be solved.
273  //
277 
278  //
279  // Algorithmic parameters
280  //
281  // numRHS_ is the current number of linear systems being solved.
282  int numRHS_;
283 
284  //
285  // Current solver state
286  //
287  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
288  // is capable of running; _initialize is controlled by the initialize() member method
289  // For the implications of the state of initialized_, please see documentation for initialize()
291 
292  // Current number of iterations performed.
293  int iter_;
294 
295  // Current number of iterations performed.
297 
298  //
299  // State Storage
300  //
301  // Residual
303  //
304  // Preconditioned residual
306  //
307  // Direction vector
309  //
310  // Operator applied to direction vector
312  //
313  // Stochastic recurrence vector
315  //
316  // Stochastic variable storage (for normal() method)
318 
319  };
320 
322  // Constructor.
323  template<class ScalarType, class MV, class OP>
325  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
327  Teuchos::ParameterList &params ):
328  lp_(problem),
329  om_(printer),
330  stest_(tester),
331  numRHS_(0),
332  initialized_(false),
333  iter_(0),
334  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
335  {
336  }
337 
338 
340  // Initialize this iteration object
341  template <class ScalarType, class MV, class OP>
343  {
344  // Check if there is any multivector to clone from.
345  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
346  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
347  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
348  "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
349 
350  // Get the multivector that is not null.
351  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
352 
353  // Get the number of right-hand sides we're solving for now.
354  int numRHS = MVT::GetNumberVecs(*tmp);
355  numRHS_ = numRHS;
356 
357  // Initialize the state storage
358  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
359  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
360  R_ = MVT::Clone( *tmp, numRHS_ );
361  Z_ = MVT::Clone( *tmp, numRHS_ );
362  P_ = MVT::Clone( *tmp, numRHS_ );
363  AP_ = MVT::Clone( *tmp, numRHS_ );
364  Y_ = MVT::Clone( *tmp, numRHS_ );
365  }
366 
367  // Initialize the random vector container with zeros.
368  randvec_.size( numRHS_ );
369 
370  // NOTE: In StochasticCGIter R_, the initial residual, is required!!!
371  //
372  std::string errstr("Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
373 
374  // Create convenience variables for zero and one.
375  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
377 
378  if (!Teuchos::is_null(newstate.R)) {
379 
380  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
381  std::invalid_argument, errstr );
382  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
383  std::invalid_argument, errstr );
384 
385  // Copy basis vectors from newstate into V
386  if (newstate.R != R_) {
387  // copy over the initial residual (unpreconditioned).
388  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
389  }
390 
391  // Compute initial direction vectors
392  // Initially, they are set to the preconditioned residuals
393 
394  if ( lp_->getLeftPrec() != Teuchos::null ) {
395  lp_->applyLeftPrec( *R_, *Z_ );
396  if ( lp_->getRightPrec() != Teuchos::null ) {
397  Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
398  lp_->applyRightPrec( *Z_, *tmp2 );
399  Z_ = tmp2;
400  }
401  }
402  else if ( lp_->getRightPrec() != Teuchos::null ) {
403  lp_->applyRightPrec( *R_, *Z_ );
404  }
405  else {
406  Z_ = R_;
407  }
408  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
409  }
410  else {
411 
412  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
413  "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
414  }
415 
416  // The solver is initialized
417  initialized_ = true;
418  }
419 
420 
422  // Iterate until the status test informs us we should stop.
423  template <class ScalarType, class MV, class OP>
425  {
426  //
427  // Allocate/initialize data structures
428  //
429  if (initialized_ == false) {
430  initialize();
431  }
432 
433  // Allocate memory for scalars.
434  int i=0;
435  std::vector<int> index(1);
436  std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
437  Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(numRHS_,numRHS_);
438 
439  // Create convenience variables for zero and one.
440  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
442 
443  // Get the current solution std::vector.
444  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
445 
446  // Compute first <r,z> a.k.a. rHz
447  MVT::MvDot( *R_, *Z_, rHz );
448 
449  if ( assertPositiveDefiniteness_ )
450  for (i=0; i<numRHS_; ++i)
451  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
453  "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
454 
456  // Iterate until the status test tells us to stop.
457  //
458  while (stest_->checkStatus(this) != Passed) {
459 
460  // Increment the iteration
461  iter_++;
462 
463  // Multiply the current direction std::vector by A and store in AP_
464  lp_->applyOp( *P_, *AP_ );
465 
466  // Compute alpha := <R_,Z_> / <P_,AP_>
467  MVT::MvDot( *P_, *AP_, pAp );
468 
470 
471  for (i=0; i<numRHS_; ++i) {
472  if ( assertPositiveDefiniteness_ )
473  // Check that pAp[i] is a positive number!
474  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
476  "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
477 
478  alpha(i,i) = rHz[i] / pAp[i];
479 
480  // Compute the scaling parameter for the stochastic vector
481  zeta(i,i) = z[i] / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
482  }
483 
484  //
485  // Update the solution std::vector x := x + alpha * P_
486  //
487  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
488  lp_->updateSolution();
489 
490  // Updates the stochastic vector y := y + zeta * P_
491  MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
492 
493  //
494  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
495  //
496  for (i=0; i<numRHS_; ++i) {
497  rHz_old[i] = rHz[i];
498  }
499  //
500  // Compute the new residual R_ := R_ - alpha * AP_
501  //
502  MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
503  //
504  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
505  // and the new direction std::vector p.
506  //
507  if ( lp_->getLeftPrec() != Teuchos::null ) {
508  lp_->applyLeftPrec( *R_, *Z_ );
509  if ( lp_->getRightPrec() != Teuchos::null ) {
510  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
511  lp_->applyRightPrec( *Z_, *tmp );
512  Z_ = tmp;
513  }
514  }
515  else if ( lp_->getRightPrec() != Teuchos::null ) {
516  lp_->applyRightPrec( *R_, *Z_ );
517  }
518  else {
519  Z_ = R_;
520  }
521  //
522  MVT::MvDot( *R_, *Z_, rHz );
523  if ( assertPositiveDefiniteness_ )
524  for (i=0; i<numRHS_; ++i)
525  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
527  "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
528  //
529  // Update the search directions.
530  for (i=0; i<numRHS_; ++i) {
531  beta(i,i) = rHz[i] / rHz_old[i];
532  index[0] = i;
533  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
534  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
535  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
536  }
537  //
538  } // end while (sTest_->checkStatus(this) != Passed)
539  }
540 
541 } // end Belos namespace
542 
543 #endif /* BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP */
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
const Teuchos::RCP< OutputManager< ScalarType > > om_
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)
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.
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
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.
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
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
Teuchos::SerialDenseVector< int, ScalarType > & normal()
Wrapper for Normal(0,1) random variables.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
TransListIter iter
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.
Teuchos::SerialDenseVector< int, ScalarType > randvec_
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Structure to contain pointers to CGIteration state variables.