Belos Version of the Day
Loading...
Searching...
No Matches
BelosPseudoBlockCGIter.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_HPP
43#define BELOS_PSEUDO_BLOCK_CG_ITER_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51#include "BelosCGIteration.hpp"
52
56#include "BelosStatusTest.hpp"
59
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
65
77namespace Belos {
78
79 template<class ScalarType, class MV, class OP>
80 class PseudoBlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
81
82 public:
83
84 //
85 // Convenience typedefs
86 //
89 typedef Teuchos::ScalarTraits<ScalarType> SCT;
90 typedef typename SCT::magnitudeType MagnitudeType;
91
93
94
100 PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
101 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
102 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
103 Teuchos::ParameterList &params );
104
106 virtual ~PseudoBlockCGIter() {};
108
109
111
112
126 void iterate();
127
149
154 {
156 initializeCG(empty);
157 }
158
168 state.R = R_;
169 state.P = P_;
170 state.AP = AP_;
171 state.Z = Z_;
172 return state;
173 }
174
176
177
179
180
182 int getNumIters() const { return iter_; }
183
185 void resetNumIters( int iter = 0 ) { iter_ = iter; }
186
189 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
190
192
194 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
195
197
199
200
202 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
203
205 int getBlockSize() const { return 1; }
206
208 void setBlockSize(int blockSize) {
209 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
210 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
211 }
212
214 bool isInitialized() { return initialized_; }
215
217
219 void setDoCondEst(bool val) {
220 if (numEntriesForCondEst_) doCondEst_=val;
221 }
222
224 Teuchos::ArrayView<MagnitudeType> getDiag() {
225 // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
226 // getDiag() didn't actually throw for me in that case, but why
227 // not be cautious?
228 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
229 if (static_cast<size_type> (iter_) >= diag_.size ()) {
230 return diag_ ();
231 } else {
232 return diag_ (0, iter_);
233 }
234 }
235
237 Teuchos::ArrayView<MagnitudeType> getOffDiag() {
238 // NOTE (mfh 30 Jul 2015) The implementation as I found it
239 // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
240 // debug mode) when the maximum number of iterations has been
241 // reached, because iter_ == offdiag_.size() in that case. The
242 // new logic fixes this.
243 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
244 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
245 return offdiag_ ();
246 } else {
247 return offdiag_ (0, iter_);
248 }
249 }
250
251 private:
252
253 //
254 // Classes inputed through constructor that define the linear problem to be solved.
255 //
256 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
257 const Teuchos::RCP<OutputManager<ScalarType> > om_;
258 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
259
260 //
261 // Algorithmic parameters
262 //
263 // numRHS_ is the current number of linear systems being solved.
264 int numRHS_;
265
266 //
267 // Current solver state
268 //
269 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
270 // is capable of running; _initialize is controlled by the initialize() member method
271 // For the implications of the state of initialized_, please see documentation for initialize()
272 bool initialized_;
273
274 // Current number of iterations performed.
275 int iter_;
276
277 // Assert that the matrix is positive definite
278 bool assertPositiveDefiniteness_;
279
280 // Tridiagonal system for condition estimation (if needed)
281 Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
282 ScalarType pAp_old_, beta_old_, rHz_old2_; // Put scalars here so that estimate is correct for multiple RHS, when deflation occurs.
283 int numEntriesForCondEst_;
284 bool doCondEst_;
285
286 //
287 // State Storage
288 //
289 // Residual
290 Teuchos::RCP<MV> R_;
291 //
292 // Preconditioned residual
293 Teuchos::RCP<MV> Z_;
294 //
295 // Direction vector
296 Teuchos::RCP<MV> P_;
297 //
298 // Operator applied to direction vector
299 Teuchos::RCP<MV> AP_;
300
301 };
302
304 // Constructor.
305 template<class ScalarType, class MV, class OP>
307 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
308 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
309 Teuchos::ParameterList &params ):
310 lp_(problem),
311 om_(printer),
312 stest_(tester),
313 numRHS_(0),
314 initialized_(false),
315 iter_(0),
316 assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
317 numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
318 doCondEst_(false)
319 {
320 }
321
322
324 // Initialize this iteration object
325 template <class ScalarType, class MV, class OP>
327 {
328 // Check if there is any mltivector to clone from.
329 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
330 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
331 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
332 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
333
334 // Get the multivector that is not null.
335 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
336
337 // Get the number of right-hand sides we're solving for now.
338 int numRHS = MVT::GetNumberVecs(*tmp);
339 numRHS_ = numRHS;
340
341 // Initialize the state storage
342 // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
343 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
344 R_ = MVT::Clone( *tmp, numRHS_ );
345 Z_ = MVT::Clone( *tmp, numRHS_ );
346 P_ = MVT::Clone( *tmp, numRHS_ );
347 AP_ = MVT::Clone( *tmp, numRHS_ );
348 }
349
350 // Tracking information for condition number estimation
351 if(numEntriesForCondEst_ > 0) {
352 diag_.resize(numEntriesForCondEst_);
353 offdiag_.resize(numEntriesForCondEst_-1);
354 }
355
356 // NOTE: In CGIter R_, the initial residual, is required!!!
357 //
358 std::string errstr("Belos::BlockPseudoCGIter::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> tmp1 = MVT::Clone( *Z_, numRHS_ );
384 lp_->applyRightPrec( *Z_, *tmp1 );
385 Z_ = tmp1;
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::CGIter::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_ ), beta( numRHS_ );
423 Teuchos::SerialDenseMatrix<int, ScalarType> alpha( 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::PseudoBlockCGIter::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::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
461
462 alpha(i,i) = rHz[i] / pAp[i];
463 }
464
465 //
466 // Update the solution std::vector x := x + alpha * P_
467 //
468 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
469 lp_->updateSolution();// what does this do?
470 //
471 // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
472 //
473 for (i=0; i<numRHS_; ++i) {
474 rHz_old[i] = rHz[i];
475 }
476 //
477 // Compute the new residual R_ := R_ - alpha * AP_
478 //
479 MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
480 //
481 // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
482 // and the new direction std::vector p.
483 //
484 if ( lp_->getLeftPrec() != Teuchos::null ) {
485 lp_->applyLeftPrec( *R_, *Z_ );
486 if ( lp_->getRightPrec() != Teuchos::null ) {
487 Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
488 lp_->applyRightPrec( *Z_, *tmp );
489 Z_ = tmp;
490 }
491 }
492 else if ( lp_->getRightPrec() != Teuchos::null ) {
493 lp_->applyRightPrec( *R_, *Z_ );
494 }
495 else {
496 Z_ = R_;
497 }
498 //
499 MVT::MvDot( *R_, *Z_, rHz );
500 if ( assertPositiveDefiniteness_ )
501 for (i=0; i<numRHS_; ++i)
502 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
504 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
505 //
506 // Update the search directions.
507 for (i=0; i<numRHS_; ++i) {
508 beta[i] = rHz[i] / rHz_old[i];
509 index[0] = i;
510 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
511 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
512 MVT::MvAddMv( one, *Z_i, beta[i], *P_i, *P_i );
513 }
514
515 // Condition estimate (if needed)
516 if (doCondEst_ && (iter_ - 1) < diag_.size()) {
517 if (iter_ > 1) {
518 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old_ * beta_old_ * pAp_old_ + pAp[0]) / rHz_old[0]);
519 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old_ * pAp_old_ / (sqrt( rHz_old[0] * rHz_old2_)));
520 }
521 else {
522 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
523 }
524 rHz_old2_ = rHz_old[0];
525 beta_old_ = beta[0];
526 pAp_old_ = pAp[0];
527 }
528
529
530 //
531 } // end while (sTest_->checkStatus(this) != Passed)
532 }
533
534} // end Belos namespace
535
536#endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
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.
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 pseudo-block CG iteration, where the basic CG algorithm is performed on all...
void resetNumIters(int iter=0)
Reset the iteration count.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
PseudoBlockCGIter(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)
PseudoBlockCGIter constructor with linear problem, solver utilities, and parameter list of solver opt...
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
virtual ~PseudoBlockCGIter()
Destructor.
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
MultiVecTraits< ScalarType, MV > MVT
void setBlockSize(int blockSize)
Set the blocksize.
OperatorTraits< ScalarType, MV, OP > OPT
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 > P
The current decent direction vector.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.

Generated for Belos by doxygen 1.9.6