Belos Package Browser (Single Doxygen Collection)  Development
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"
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 
82 namespace Belos {
83 
84  template<class ScalarType, class MV, class OP>
85  class PseudoBlockStochasticCGIter : virtual public StochasticCGIteration<ScalarType,MV,OP> {
86 
87  public:
88 
89  //
90  // Convenience typedefs
91  //
94  typedef Teuchos::ScalarTraits<ScalarType> SCT;
95  typedef typename SCT::magnitudeType MagnitudeType;
96 
98 
99 
105  PseudoBlockStochasticCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
106  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
107  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
108  Teuchos::ParameterList &params );
109 
113 
114 
116 
117 
131  void iterate();
132 
154 
158  void initialize()
159  {
161  initializeCG(empty);
162  }
163 
173  state.R = R_;
174  state.P = P_;
175  state.AP = AP_;
176  state.Z = Z_;
177  state.Y = Y_;
178  return state;
179  }
180 
182 
183 
185 
186 
188  int getNumIters() const { return iter_; }
189 
191  void resetNumIters( int iter = 0 ) { iter_ = iter; }
192 
195  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
196 
198 
200  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
201 
203  Teuchos::RCP<MV> getStochasticVector() const { return Y_; }
204 
206 
208 
209 
211  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
212 
214  int getBlockSize() const { return 1; }
215 
217  void setBlockSize(int blockSize) {
218  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
219  "Belos::PseudoBlockStochasticCGIter::setBlockSize(): Cannot use a block size that is not one.");
220  }
221 
223  bool isInitialized() { return initialized_; }
224 
226 
227  private:
228 
230  inline ScalarType normal() {
231  // Do all of the calculations with doubles, because that is what the Odeh and Evans 1974 constants are for.
232  // Then cast to ScalarType.
233 
234  const double p0 = -0.322232431088;
235  const double p1 = -1.0;
236  const double p2 = -0.342242088547;
237  const double p3 = -0.204231210245e-1;
238  const double p4 = -0.453642210148e-4;
239  const double q0 = 0.993484626060e-1;
240  const double q1 = 0.588581570495;
241  const double q2 = 0.531103462366;
242  const double q3 = 0.103537752850;
243  const double q4 = 0.38560700634e-2;
244  double r,p,q,y,z;
245 
246  // Get a random number (-1,1) and rescale to (0,1).
247  r=0.5*(Teuchos::ScalarTraits<double>::random() + 1.0);
248 
249  // Odeh and Evans algorithm (as modified by Park & Geyer)
250  if(r < 0.5) y=std::sqrt(-2.0 * log(r));
251  else y=std::sqrt(-2.0 * log(1.0 - r));
252 
253  p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
254  q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
255 
256  if(r < 0.5) z = (p / q) - y;
257  else z = y - (p / q);
258 
259  return Teuchos::as<ScalarType,double>(z);
260  }
261 
262  //
263  // Classes inputed through constructor that define the linear problem to be solved.
264  //
265  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
266  const Teuchos::RCP<OutputManager<ScalarType> > om_;
267  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
268 
269  //
270  // Algorithmic parameters
271  //
272  // numRHS_ is the current number of linear systems being solved.
273  int numRHS_;
274 
275  //
276  // Current solver state
277  //
278  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
279  // is capable of running; _initialize is controlled by the initialize() member method
280  // For the implications of the state of initialized_, please see documentation for initialize()
282 
283  // Current number of iterations performed.
284  int iter_;
285 
286  // Current number of iterations performed.
288 
289  //
290  // State Storage
291  //
292  // Residual
293  Teuchos::RCP<MV> R_;
294  //
295  // Preconditioned residual
296  Teuchos::RCP<MV> Z_;
297  //
298  // Direction vector
299  Teuchos::RCP<MV> P_;
300  //
301  // Operator applied to direction vector
302  Teuchos::RCP<MV> AP_;
303  //
304  // Stochastic recurrence vector
305  Teuchos::RCP<MV> Y_;
306 
307 
308  };
309 
311  // Constructor.
312  template<class ScalarType, class MV, class OP>
314  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
315  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
316  Teuchos::ParameterList &params ):
317  lp_(problem),
318  om_(printer),
319  stest_(tester),
320  numRHS_(0),
321  initialized_(false),
322  iter_(0),
323  assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
324  {
325  }
326 
327 
329  // Initialize this iteration object
330  template <class ScalarType, class MV, class OP>
332  {
333  // Check if there is any multivector to clone from.
334  Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
335  Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
336  TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
337  "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
338 
339  // Get the multivector that is not null.
340  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
341 
342  // Get the number of right-hand sides we're solving for now.
343  int numRHS = MVT::GetNumberVecs(*tmp);
344  numRHS_ = numRHS;
345 
346  // Initialize the state storage
347  // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
348  if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
349  R_ = MVT::Clone( *tmp, numRHS_ );
350  Z_ = MVT::Clone( *tmp, numRHS_ );
351  P_ = MVT::Clone( *tmp, numRHS_ );
352  AP_ = MVT::Clone( *tmp, numRHS_ );
353  Y_ = MVT::Clone( *tmp, numRHS_ );
354  }
355 
356  // NOTE: In StochasticCGIter R_, the initial residual, is required!!!
357  //
358  std::string errstr("Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
359 
360  // Create convenience variables for zero and one.
361  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
362  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
363 
364  if (!Teuchos::is_null(newstate.R)) {
365 
366  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
367  std::invalid_argument, errstr );
368  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
369  std::invalid_argument, errstr );
370 
371  // Copy basis vectors from newstate into V
372  if (newstate.R != R_) {
373  // copy over the initial residual (unpreconditioned).
374  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
375  }
376 
377  // Compute initial direction vectors
378  // Initially, they are set to the preconditioned residuals
379 
380  if ( lp_->getLeftPrec() != Teuchos::null ) {
381  lp_->applyLeftPrec( *R_, *Z_ );
382  if ( lp_->getRightPrec() != Teuchos::null ) {
383  Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
384  lp_->applyRightPrec( *Z_, *tmp2 );
385  Z_ = tmp2;
386  }
387  }
388  else if ( lp_->getRightPrec() != Teuchos::null ) {
389  lp_->applyRightPrec( *R_, *Z_ );
390  }
391  else {
392  Z_ = R_;
393  }
394  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
395  }
396  else {
397 
398  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
399  "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
400  }
401 
402  // The solver is initialized
403  initialized_ = true;
404  }
405 
406 
408  // Iterate until the status test informs us we should stop.
409  template <class ScalarType, class MV, class OP>
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_ ), zeta(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  // 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::PseudoBlockStochasticCGIter::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::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
461 
462  alpha(i,i) = rHz[i] / pAp[i];
463 
464 
465  // Compute the scaling parameter for the stochastic vector
466  ScalarType z = normal();
467  zeta(i,i) = z / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
468  }
469 
470  //
471  // Update the solution std::vector x := x + alpha * P_
472  //
473  MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
474  lp_->updateSolution();
475 
476  // Updates the stochastic vector y := y + zeta * P_
477  MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
478 
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,
513  "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
514  //
515  // Update the search directions.
516  for (i=0; i<numRHS_; ++i) {
517  beta(i,i) = rHz[i] / rHz_old[i];
518  index[0] = i;
519  Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
520  Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
521  MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
522  }
523  //
524  } // end while (sTest_->checkStatus(this) != Passed)
525  }
526 
527 } // end Belos namespace
528 
529 #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.
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
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...
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.
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.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector.
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.
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.
ScalarType normal()
Wrapper for Normal(0,1) random variables.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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...
int getNumIters() const
Get the current iteration count.
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.