Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Belos_PseudoBlockCGIter_MP_Vector.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_MP_VECTOR_HPP
43#define BELOS_PSEUDO_BLOCK_CG_ITER_MP_VECTOR_HPP
44
54#include "BelosPseudoBlockCGIter.hpp"
55
71namespace Belos {
72
73 template<class Storage, class MV, class OP>
74 class PseudoBlockCGIter<Sacado::MP::Vector<Storage>, MV, OP> :
75 virtual public CGIteration<Sacado::MP::Vector<Storage>, MV, OP> {
76
77 public:
78
79 //
80 // Convenience typedefs
81 //
83 typedef MultiVecTraits<ScalarType,MV> MVT;
84 typedef OperatorTraits<ScalarType,MV,OP> OPT;
85 typedef Teuchos::ScalarTraits<ScalarType> SCT;
86 typedef typename SCT::magnitudeType MagnitudeType;
87 typedef Teuchos::ScalarTraits<typename Storage::value_type> SVT;
88
90
91
97 PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
98 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
99 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
100 Teuchos::ParameterList &params );
101
103 virtual ~PseudoBlockCGIter() {};
105
106
108
109
123 void iterate();
124
145 void initializeCG(CGIterationState<ScalarType,MV>& newstate);
146
151 {
152 CGIterationState<ScalarType,MV> empty;
153 initializeCG(empty);
154 }
155
163 CGIterationState<ScalarType,MV> getState() const {
164 CGIterationState<ScalarType,MV> state;
165 state.R = R_;
166 state.P = P_;
167 state.AP = AP_;
168 state.Z = Z_;
169 return state;
170 }
171
173
174
176
177
179 int getNumIters() const { return iter_; }
180
182 void resetNumIters( int iter = 0 ) { iter_ = iter; }
183
186 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
187
189
191 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
192
194
196
197
199 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
200
202 int getBlockSize() const { return 1; }
203
205 void setBlockSize(int blockSize) {
206 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
207 "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
208 }
209
211 bool isInitialized() { return initialized_; }
212
214
216 void setDoCondEst(bool val){doCondEst_=val;}
217
219 Teuchos::ArrayView<MagnitudeType> getDiag() {
220 // NOTE (mfh 30 Jul 2015) See note on getOffDiag() below.
221 // getDiag() didn't actually throw for me in that case, but why
222 // not be cautious?
223 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
224 if (static_cast<size_type> (iter_) >= diag_.size ()) {
225 return diag_ ();
226 } else {
227 return diag_ (0, iter_);
228 }
229 }
230
232 Teuchos::ArrayView<MagnitudeType> getOffDiag() {
233 // NOTE (mfh 30 Jul 2015) The implementation as I found it
234 // returned "offdiag(0,iter_)". This breaks (Teuchos throws in
235 // debug mode) when the maximum number of iterations has been
236 // reached, because iter_ == offdiag_.size() in that case. The
237 // new logic fixes this.
238 typedef typename Teuchos::ArrayView<MagnitudeType>::size_type size_type;
239 if (static_cast<size_type> (iter_) >= offdiag_.size ()) {
240 return offdiag_ ();
241 } else {
242 return offdiag_ (0, iter_);
243 }
244 }
245
246 private:
247
248 //
249 // Classes inputed through constructor that define the linear problem to be solved.
250 //
251 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
252 const Teuchos::RCP<OutputManager<ScalarType> > om_;
253 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
254
255 //
256 // Algorithmic parameters
257 //
258 // numRHS_ is the current number of linear systems being solved.
260
261 //
262 // Current solver state
263 //
264 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
265 // is capable of running; _initialize is controlled by the initialize() member method
266 // For the implications of the state of initialized_, please see documentation for initialize()
268
269 // Current number of iterations performed.
270 int iter_;
271
272 // Assert that the matrix is positive definite
274
275 // Tridiagonal system for condition estimation (if needed)
276 Teuchos::ArrayRCP<MagnitudeType> diag_, offdiag_;
279
280 //
281 // State Storage
282 //
283 // Residual
284 Teuchos::RCP<MV> R_;
285 //
286 // Preconditioned residual
287 Teuchos::RCP<MV> Z_;
288 //
289 // Direction vector
290 Teuchos::RCP<MV> P_;
291 //
292 // Operator applied to direction vector
293 Teuchos::RCP<MV> AP_;
294
295 // Tolerance for ensemble breakdown
296 typename SVT::magnitudeType breakDownTol_;
297
298 };
299
301 // Constructor.
302 template<class StorageType, class MV, class OP>
304 PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
305 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
306 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
307 Teuchos::ParameterList &params ):
308 lp_(problem),
309 om_(printer),
310 stest_(tester),
311 numRHS_(0),
312 initialized_(false),
313 iter_(0),
314 assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) ),
315 numEntriesForCondEst_(params.get("Max Size For Condest",0) ),
316 breakDownTol_(params.get("Ensemble Breakdown Tolerance", 0.0))
317 {
318 }
319
320
322 // Initialize this iteration object
323 template <class StorageType, class MV, class OP>
325 initializeCG(CGIterationState<ScalarType,MV>& newstate)
326 {
327 // Check if there is any mltivector to clone from.
328 Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
329 Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
330 TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
331 "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
332
333 // Get the multivector that is not null.
334 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
335
336 // Get the number of right-hand sides we're solving for now.
337 int numRHS = MVT::GetNumberVecs(*tmp);
338 numRHS_ = numRHS;
339
340 // Initialize the state storage
341 // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
342 if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
343 R_ = MVT::Clone( *tmp, numRHS_ );
344 Z_ = MVT::Clone( *tmp, numRHS_ );
345 P_ = MVT::Clone( *tmp, numRHS_ );
346 AP_ = MVT::Clone( *tmp, numRHS_ );
347 }
348
349 // Tracking information for condition number estimation
350 if(numEntriesForCondEst_ > 0) {
351 diag_.resize(numEntriesForCondEst_);
352 offdiag_.resize(numEntriesForCondEst_-1);
353 }
354
355 // NOTE: In CGIter R_, the initial residual, is required!!!
356 //
357 std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
358
359 // Create convenience variables for zero and one.
360 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
361 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
362
363 if (!Teuchos::is_null(newstate.R)) {
364
365 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
366 std::invalid_argument, errstr );
367 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
368 std::invalid_argument, errstr );
369
370 // Copy basis vectors from newstate into V
371 if (newstate.R != R_) {
372 // copy over the initial residual (unpreconditioned).
373 MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
374 }
375
376 // Compute initial direction vectors
377 // Initially, they are set to the preconditioned residuals
378 //
379 if ( lp_->getLeftPrec() != Teuchos::null ) {
380 lp_->applyLeftPrec( *R_, *Z_ );
381 if ( lp_->getRightPrec() != Teuchos::null ) {
382 Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
383 lp_->applyRightPrec( *Z_, *tmp1 );
384 Z_ = tmp1;
385 }
386 }
387 else if ( lp_->getRightPrec() != Teuchos::null ) {
388 lp_->applyRightPrec( *R_, *Z_ );
389 }
390 else {
391 Z_ = R_;
392 }
393 MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
394 }
395 else {
396
397 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
398 "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
399 }
400
401 // The solver is initialized
402 initialized_ = true;
403 }
404
405
407 // Iterate until the status test informs us we should stop.
408 template <class StorageType, class MV, class OP>
410 iterate()
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_ );
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 // Scalars for condition estimation (if needed) - These will always use entry zero, for convenience
430 ScalarType pAp_old=one, beta_old=one ,rHz_old2=one;
431
432 // Get the current solution std::vector.
433 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
434
435 // Compute first <r,z> a.k.a. rHz
436 MVT::MvDot( *R_, *Z_, rHz );
437
438 if ( assertPositiveDefiniteness_ )
439 for (i=0; i<numRHS_; ++i)
440 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
441 CGIterateFailure,
442 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
443
445 // Iterate until the status test tells us to stop.
446 //
447 while (stest_->checkStatus(this) != Passed) {
448
449 // Increment the iteration
450 iter_++;
451
452 // Multiply the current direction std::vector by A and store in AP_
453 lp_->applyOp( *P_, *AP_ );
454
455 // Compute alpha := <R_,Z_> / <P_,AP_>
456 MVT::MvDot( *P_, *AP_, pAp );
457
458 for (i=0; i<numRHS_; ++i) {
459 if ( assertPositiveDefiniteness_ )
460 // Check that pAp[i] is a positive number!
461 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
462 CGIterateFailure,
463 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
464
465 const int ensemble_size = pAp[i].size();
466 for (int k=0; k<ensemble_size; ++k) {
467 if (SVT::magnitude(pAp[i].coeff(k)) <= breakDownTol_)
468 alpha(i,i).fastAccessCoeff(k) = 0.0;
469 else
470 alpha(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / pAp[i].coeff(k);
471 }
472 }
473
474 //
475 // Update the solution std::vector x := x + alpha * P_
476 //
477 MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
478 lp_->updateSolution();
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,
512 CGIterateFailure,
513 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
514 //
515 // Update the search directions.
516 for (i=0; i<numRHS_; ++i) {
517 const int ensemble_size = rHz[i].size();
518 for (int k=0; k<ensemble_size; ++k) {
519 if (SVT::magnitude(rHz[i].coeff(k)) <= breakDownTol_)
520 beta(i,i).fastAccessCoeff(k) = 0.0;
521 else
522 beta(i,i).fastAccessCoeff(k) = rHz[i].coeff(k) / rHz_old[i].coeff(k);
523 }
524 index[0] = i;
525 Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
526 Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
527 MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
528 }
529
530 // Condition estimate (if needed)
531 if(doCondEst_ > 0) {
532 if(iter_ > 1 ) {
533 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real((beta_old * beta_old * pAp_old + pAp[0]) / rHz_old[0]);
534 offdiag_[iter_-2] = -Teuchos::ScalarTraits<ScalarType>::real(beta_old * pAp_old / (sqrt( rHz_old[0] * rHz_old2)));
535 }
536 else {
537 diag_[iter_-1] = Teuchos::ScalarTraits<ScalarType>::real(pAp[0] / rHz_old[0]);
538 }
539 rHz_old2 = rHz_old[0];
540 beta_old = beta(0,0);
541 pAp_old = pAp[0];
542 }
543
544
545 //
546 } // end while (sTest_->checkStatus(this) != Passed)
547 }
548
549} // end Belos namespace
550
551#endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
expr val()
void iterate()
This method performs CG iterations on each linear system until the status test indicates the need to ...
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...
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
bool isInitialized()
States whether the solver has been initialized or not.
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.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)