Belos Version of the Day
Loading...
Searching...
No Matches
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
56#include "BelosStatusTest.hpp"
59
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_SerialDenseHelpers.hpp"
63#include "Teuchos_ScalarTraits.hpp"
64#include "Teuchos_ParameterList.hpp"
65#include "Teuchos_TimeMonitor.hpp"
66
82namespace 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
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 Teuchos::SerialDenseVector<int,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 // Return a vector with random entries that are synchronized across processors.
247 Teuchos::randomSyncedMatrix( randvec_ );
248
249 for (int i=0; i<numRHS_; i++)
250 {
251 // Get a random number (-1,1) and rescale to (0,1).
252 r=0.5*SCT::real(randvec_[i]) + 1.0;
253
254 // Odeh and Evans algorithm (as modified by Park & Geyer)
255 if(r < 0.5) y=std::sqrt(-2.0 * log(r));
256 else y=std::sqrt(-2.0 * log(1.0 - r));
257
258 p = p0 + y * (p1 + y* (p2 + y * (p3 + y * p4)));
259 q = q0 + y * (q1 + y* (q2 + y * (q3 + y * q4)));
260
261 if(r < 0.5) z = (p / q) - y;
262 else z = y - (p / q);
263
264 randvec_[i] = Teuchos::as<ScalarType,double>(z);
265 }
266
267 return randvec_;
268 }
269
270 //
271 // Classes inputed through constructor that define the linear problem to be solved.
272 //
273 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
274 const Teuchos::RCP<OutputManager<ScalarType> > om_;
275 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
276
277 //
278 // Algorithmic parameters
279 //
280 // numRHS_ is the current number of linear systems being solved.
281 int numRHS_;
282
283 //
284 // Current solver state
285 //
286 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
287 // is capable of running; _initialize is controlled by the initialize() member method
288 // For the implications of the state of initialized_, please see documentation for initialize()
289 bool initialized_;
290
291 // Current number of iterations performed.
292 int iter_;
293
294 // Current number of iterations performed.
295 bool assertPositiveDefiniteness_;
296
297 //
298 // State Storage
299 //
300 // Residual
301 Teuchos::RCP<MV> R_;
302 //
303 // Preconditioned residual
304 Teuchos::RCP<MV> Z_;
305 //
306 // Direction vector
307 Teuchos::RCP<MV> P_;
308 //
309 // Operator applied to direction vector
310 Teuchos::RCP<MV> AP_;
311 //
312 // Stochastic recurrence vector
313 Teuchos::RCP<MV> Y_;
314 //
315 // Stochastic variable storage (for normal() method)
316 Teuchos::SerialDenseVector<int,ScalarType> randvec_;
317
318 };
319
321 // Constructor.
322 template<class ScalarType, class MV, class OP>
324 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
325 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
326 Teuchos::ParameterList &params ):
327 lp_(problem),
328 om_(printer),
329 stest_(tester),
330 numRHS_(0),
331 initialized_(false),
332 iter_(0),
333 assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
334 {
335 }
336
337
339 // Initialize this iteration object
340 template <class ScalarType, class MV, class OP>
342 {
343 // Check if there is any multivector to clone from.
344 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
345 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
346 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
347 "Belos::PseudoBlockStochasticCGIter::initialize(): Cannot initialize state storage!");
348
349 // Get the multivector that is not null.
350 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
351
352 // Get the number of right-hand sides we're solving for now.
353 int numRHS = MVT::GetNumberVecs(*tmp);
354 numRHS_ = numRHS;
355
356 // Initialize the state storage
357 // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
358 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
359 R_ = MVT::Clone( *tmp, numRHS_ );
360 Z_ = MVT::Clone( *tmp, numRHS_ );
361 P_ = MVT::Clone( *tmp, numRHS_ );
362 AP_ = MVT::Clone( *tmp, numRHS_ );
363 Y_ = MVT::Clone( *tmp, numRHS_ );
364 }
365
366 // Initialize the random vector container with zeros.
367 randvec_.size( numRHS_ );
368
369 // NOTE: In StochasticCGIter R_, the initial residual, is required!!!
370 //
371 std::string errstr("Belos::BlockPseudoStochasticCGIter::initialize(): Specified multivectors must have a consistent length and width.");
372
373 // Create convenience variables for zero and one.
374 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
375 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
376
377 if (!Teuchos::is_null(newstate.R)) {
378
379 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
380 std::invalid_argument, errstr );
381 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
382 std::invalid_argument, errstr );
383
384 // Copy basis vectors from newstate into V
385 if (newstate.R != R_) {
386 // copy over the initial residual (unpreconditioned).
387 MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
388 }
389
390 // Compute initial direction vectors
391 // Initially, they are set to the preconditioned residuals
392
393 if ( lp_->getLeftPrec() != Teuchos::null ) {
394 lp_->applyLeftPrec( *R_, *Z_ );
395 if ( lp_->getRightPrec() != Teuchos::null ) {
396 Teuchos::RCP<MV> tmp2 = MVT::Clone( *Z_, numRHS_ );
397 lp_->applyRightPrec( *Z_, *tmp2 );
398 Z_ = tmp2;
399 }
400 }
401 else if ( lp_->getRightPrec() != Teuchos::null ) {
402 lp_->applyRightPrec( *R_, *Z_ );
403 }
404 else {
405 Z_ = R_;
406 }
407 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
408 }
409 else {
410
411 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
412 "Belos::StochasticCGIter::initialize(): CGStateIterState does not have initial residual.");
413 }
414
415 // The solver is initialized
416 initialized_ = true;
417 }
418
419
421 // Iterate until the status test informs us we should stop.
422 template <class ScalarType, class MV, class OP>
424 {
425 //
426 // Allocate/initialize data structures
427 //
428 if (initialized_ == false) {
429 initialize();
430 }
431
432 // Allocate memory for scalars.
433 int i=0;
434 std::vector<int> index(1);
435 std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
436 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ ), zeta(numRHS_,numRHS_);
437
438 // Create convenience variables for zero and one.
439 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
441
442 // Get the current solution std::vector.
443 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
444
445 // Compute first <r,z> a.k.a. rHz
446 MVT::MvDot( *R_, *Z_, rHz );
447
448 if ( assertPositiveDefiniteness_ )
449 for (i=0; i<numRHS_; ++i)
450 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
452 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
453
455 // Iterate until the status test tells us to stop.
456 //
457 while (stest_->checkStatus(this) != Passed) {
458
459 // Increment the iteration
460 iter_++;
461
462 // Multiply the current direction std::vector by A and store in AP_
463 lp_->applyOp( *P_, *AP_ );
464
465 // Compute alpha := <R_,Z_> / <P_,AP_>
466 MVT::MvDot( *P_, *AP_, pAp );
467
468 Teuchos::SerialDenseVector<int,ScalarType>& z = normal();
469
470 for (i=0; i<numRHS_; ++i) {
471 if ( assertPositiveDefiniteness_ )
472 // Check that pAp[i] is a positive number!
473 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
475 "Belos::PseudoBlockStochasticCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
476
477 alpha(i,i) = rHz[i] / pAp[i];
478
479 // Compute the scaling parameter for the stochastic vector
480 zeta(i,i) = z[i] / Teuchos::ScalarTraits<ScalarType>::squareroot(pAp[i]);
481 }
482
483 //
484 // Update the solution std::vector x := x + alpha * P_
485 //
486 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
487 lp_->updateSolution();
488
489 // Updates the stochastic vector y := y + zeta * P_
490 MVT::MvTimesMatAddMv( one, *P_, zeta, one, *Y_);
491
492 //
493 // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
494 //
495 for (i=0; i<numRHS_; ++i) {
496 rHz_old[i] = rHz[i];
497 }
498 //
499 // Compute the new residual R_ := R_ - alpha * AP_
500 //
501 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
502 //
503 // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
504 // and the new direction std::vector p.
505 //
506 if ( lp_->getLeftPrec() != Teuchos::null ) {
507 lp_->applyLeftPrec( *R_, *Z_ );
508 if ( lp_->getRightPrec() != Teuchos::null ) {
509 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
510 lp_->applyRightPrec( *Z_, *tmp );
511 Z_ = tmp;
512 }
513 }
514 else if ( lp_->getRightPrec() != Teuchos::null ) {
515 lp_->applyRightPrec( *R_, *Z_ );
516 }
517 else {
518 Z_ = R_;
519 }
520 //
521 MVT::MvDot( *R_, *Z_, rHz );
522 if ( assertPositiveDefiniteness_ )
523 for (i=0; i<numRHS_; ++i)
524 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
526 "Belos::PseudoBlockStochasticCGIter::iterate(): negative value for r^H*M*r encountered!" );
527 //
528 // Update the search directions.
529 for (i=0; i<numRHS_; ++i) {
530 beta(i,i) = rHz[i] / rHz_old[i];
531 index[0] = i;
532 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
533 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
534 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
535 }
536 //
537 } // end while (sTest_->checkStatus(this) != Passed)
538 }
539
540} // end Belos namespace
541
542#endif /* BELOS_PSEUDO_BLOCK_STOCHASTIC_CG_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
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.
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the stochastic pseudo-block CG iteration, where the basic stochastic CG algorit...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
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 ...
bool isInitialized()
States whether the solver has been initialized or not.
int getNumIters() const
Get the current iteration count.
void initializeCG(StochasticCGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) 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.
StochasticCGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getStochasticVector() const
Get the stochastic vector
void resetNumIters(int iter=0)
Reset the iteration count.
void iterate()
This method performs stochastic CG iterations on each linear system until the status test indicates t...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void setBlockSize(int blockSize)
Set the blocksize.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > Y
The current stochastic recurrence vector.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > R
The current residual.

Generated for Belos by doxygen 1.10.0