Belos Version of the Day
Loading...
Searching...
No Matches
BelosBiCGStabIter.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_BICGSTAB_ITER_HPP
43#define BELOS_BICGSTAB_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
80
81
86 template <class ScalarType, class MV>
88
90 Teuchos::RCP<const MV> R;
91
93 Teuchos::RCP<const MV> Rhat;
94
96 Teuchos::RCP<const MV> P;
97
99 Teuchos::RCP<const MV> V;
100
101 std::vector<ScalarType> rho_old, alpha, omega;
102
103 BiCGStabIterationState() : R(Teuchos::null), Rhat(Teuchos::null),
104 P(Teuchos::null), V(Teuchos::null)
105 {
106 rho_old.clear();
107 alpha.clear();
108 omega.clear();
109 }
110 };
111
112 template<class ScalarType, class MV, class OP>
113 class BiCGStabIter : virtual public Iteration<ScalarType,MV,OP> {
114
115 public:
116
117 //
118 // Convenience typedefs
119 //
122 typedef Teuchos::ScalarTraits<ScalarType> SCT;
123 typedef typename SCT::magnitudeType MagnitudeType;
124 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
125
127
128
134 BiCGStabIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
135 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
136 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
137 Teuchos::ParameterList &params );
138
140 virtual ~BiCGStabIter() {};
142
143
145
146
160 void iterate();
161
183
188 {
190 initializeBiCGStab(empty);
191 }
192
202 state.R = R_;
203 state.Rhat = Rhat_;
204 state.P = P_;
205 state.V = V_;
206 state.rho_old = rho_old_;
207 state.alpha = alpha_;
208 state.omega = omega_;
209 return state;
210 }
211
213
214
216
217
219 int getNumIters() const { return iter_; }
220
222 void resetNumIters( int iter = 0 ) { iter_ = iter; }
223
226 // amk TODO: are the residuals actually being set? What is a native residual?
227 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
228
230
232 // amk TODO: what is this supposed to be doing?
233 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
234
236 bool breakdownDetected() { return breakdown_; }
237
239
241
242
244 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
245
247 int getBlockSize() const { return 1; }
248
250 void setBlockSize(int blockSize) {
251 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
252 "Belos::BiCGStabIter::setBlockSize(): Cannot use a block size that is not one.");
253 }
254
256 bool isInitialized() { return initialized_; }
257
259
260 private:
261
262 void axpy(const ScalarType alpha, const MV & A,
263 const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus=false);
264
265 //
266 // Classes inputed through constructor that define the linear problem to be solved.
267 //
268 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
269 const Teuchos::RCP<OutputManager<ScalarType> > om_;
270 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
271
272 //
273 // Algorithmic parameters
274 //
275 // numRHS_ is the current number of linear systems being solved.
276 int numRHS_;
277
278 //
279 // Current solver state
280 //
281 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
282 // is capable of running; _initialize is controlled by the initialize() member method
283 // For the implications of the state of initialized_, please see documentation for initialize()
284 bool initialized_;
285
286 // Breakdown has been observed for at least one of the linear systems
287 bool breakdown_;
288
289 // Current number of iterations performed.
290 int iter_;
291
292 //
293 // State Storage
294 //
295 // Initial residual
296 Teuchos::RCP<MV> Rhat_;
297 //
298 // Residual
299 Teuchos::RCP<MV> R_;
300 //
301 // Direction vector 1
302 Teuchos::RCP<MV> P_;
303 //
304 // Operator applied to preconditioned direction vector 1
305 Teuchos::RCP<MV> V_;
306 //
307 std::vector<ScalarType> rho_old_, alpha_, omega_;
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 breakdown_(false),
323 iter_(0)
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::BiCGStabIter::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 Rhat_ = MVT::Clone( *tmp, numRHS_ );
351 P_ = MVT::Clone( *tmp, numRHS_ );
352 V_ = MVT::Clone( *tmp, numRHS_ );
353
354 rho_old_.resize(numRHS_);
355 alpha_.resize(numRHS_);
356 omega_.resize(numRHS_);
357 }
358
359 // Reset breakdown to false before initializing iteration
360 breakdown_ = false;
361
362 // NOTE: In BiCGStabIter R_, the initial residual, is required!!!
363 //
364 std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
365
366 // Create convenience variable for one.
367 const ScalarType one = SCT::one();
368
369 if (!Teuchos::is_null(newstate.R)) {
370
371 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
372 std::invalid_argument, errstr );
373 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
374 std::invalid_argument, errstr );
375
376 // Copy residual vectors from newstate into R
377 if (newstate.R != R_) {
378 // Assigned by the new state
379 MVT::Assign(*newstate.R, *R_);
380 }
381 else {
382 // Computed
383 lp_->computeCurrResVec(R_.get());
384 }
385
386 // Set Rhat
387 if (!Teuchos::is_null(newstate.Rhat) && newstate.Rhat != Rhat_) {
388 // Assigned by the new state
389 MVT::Assign(*newstate.Rhat, *Rhat_);
390 }
391 else {
392 // Set to be the initial residual
393 MVT::Assign(*R_, *Rhat_);
394 }
395
396 // Set V
397 if (!Teuchos::is_null(newstate.V) && newstate.V != V_) {
398 // Assigned by the new state
399 MVT::Assign(*newstate.V, *V_);
400 }
401 else {
402 // Initial V = 0
403 MVT::MvInit(*V_);
404 }
405
406 // Set P
407 if (!Teuchos::is_null(newstate.P) && newstate.P != P_) {
408 // Assigned by the new state
409 MVT::Assign(*newstate.P, *P_);
410 }
411 else {
412 // Initial P = 0
413 MVT::MvInit(*P_);
414 }
415
416 // Set rho_old
417 if (newstate.rho_old.size () == static_cast<size_t> (numRHS_)) {
418 // Assigned by the new state
419 rho_old_ = newstate.rho_old;
420 }
421 else {
422 // Initial rho = 1
423 rho_old_.assign(numRHS_,one);
424 }
425
426 // Set alpha
427 if (newstate.alpha.size() == static_cast<size_t> (numRHS_)) {
428 // Assigned by the new state
429 alpha_ = newstate.alpha;
430 }
431 else {
432 // Initial rho = 1
433 alpha_.assign(numRHS_,one);
434 }
435
436 // Set omega
437 if (newstate.omega.size() == static_cast<size_t> (numRHS_)) {
438 // Assigned by the new state
439 omega_ = newstate.omega;
440 }
441 else {
442 // Initial rho = 1
443 omega_.assign(numRHS_,one);
444 }
445
446 }
447 else {
448
449 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
450 "Belos::BiCGStabIter::initialize(): BiCGStabStateIterState does not have initial residual.");
451 }
452
453 // The solver is initialized
454 initialized_ = true;
455 }
456
457
459 // Iterate until the status test informs us we should stop.
460 template <class ScalarType, class MV, class OP>
462 {
463 using Teuchos::RCP;
464
465 //
466 // Allocate/initialize data structures
467 //
468 if (initialized_ == false) {
469 initialize();
470 }
471
472 // Allocate memory for scalars.
473 int i=0;
474 std::vector<ScalarType> rho_new( numRHS_ ), beta( numRHS_ );
475 std::vector<ScalarType> rhatV( numRHS_ ), tT( numRHS_ ), tS( numRHS_ );
476
477 // Create convenience variable for one.
478 const ScalarType one = SCT::one();
479
480 // TODO: We may currently be using more space than is required
481 RCP<MV> leftPrecVec, leftPrecVec2;
482
483 RCP<MV> Y, Z, S, T;
484 S = MVT::Clone( *R_, numRHS_ );
485 T = MVT::Clone( *R_, numRHS_ );
486 if (lp_->isLeftPrec() || lp_->isRightPrec()) {
487 Y = MVT::Clone( *R_, numRHS_ );
488 Z = MVT::Clone( *R_, numRHS_ );
489 }
490 else {
491 Y = P_;
492 Z = S;
493 }
494
495 // Get the current solution std::vector.
496 Teuchos::RCP<MV> X = lp_->getCurrLHSVec();
497
499 // Iterate until the status test tells us to stop.
500 //
501 while (stest_->checkStatus(this) != Passed && !breakdown_) {
502
503 // Increment the iteration
504 iter_++;
505
506 // rho_new = <R_, Rhat_>
507 MVT::MvDot(*R_,*Rhat_,rho_new);
508
509 // beta = ( rho_new / rho_old ) (alpha / omega )
510 // TODO: None of these loops are currently threaded
511 for(i=0; i<numRHS_; i++) {
512 // Catch breakdown in rho_old here, since
513 // it is just rho_new from the previous iteration.
514 if (SCT::magnitude(rho_new[i]) < MT::sfmin())
515 breakdown_ = true;
516
517 beta[i] = (rho_new[i] / rho_old_[i]) * (alpha_[i] / omega_[i]);
518 }
519
520 // p = r + beta (p - omega v)
521 // TODO: Is it safe to call MvAddMv with A or B = mv?
522 // TODO: Not all of these things have to be part of the state
523 axpy(one, *P_, omega_, *V_, *P_, true); // p = p - omega v
524 axpy(one, *R_, beta, *P_, *P_); // p = r + beta (p - omega v)
525
526 // y = K\p, unless K does not exist
527 // TODO: There may be a more efficient way to apply the preconditioners
528 if(lp_->isLeftPrec()) {
529 if(lp_->isRightPrec()) {
530 if(leftPrecVec == Teuchos::null) {
531 leftPrecVec = MVT::Clone( *R_, numRHS_ );
532 }
533 lp_->applyLeftPrec(*P_,*leftPrecVec);
534 lp_->applyRightPrec(*leftPrecVec,*Y);
535 }
536 else {
537 lp_->applyLeftPrec(*P_,*Y);
538 }
539 }
540 else if(lp_->isRightPrec()) {
541 lp_->applyRightPrec(*P_,*Y);
542 }
543
544 // v = Ay
545 lp_->applyOp(*Y,*V_);
546
547 // alpha = rho_new / <Rhat, V>
548 MVT::MvDot(*V_,*Rhat_,rhatV);
549 for(i=0; i<numRHS_; i++) {
550 if (SCT::magnitude(rhatV[i]) < MT::sfmin())
551 {
552 breakdown_ = true;
553 return;
554 }
555 else
556 alpha_[i] = rho_new[i] / rhatV[i];
557 }
558
559 // s = r - alpha v
560 axpy(one, *R_, alpha_, *V_, *S, true);
561
562 // z = K\s, unless K does not exist
563 if(lp_->isLeftPrec()) {
564 if(lp_->isRightPrec()) {
565 if(leftPrecVec == Teuchos::null) {
566 leftPrecVec = MVT::Clone( *R_, numRHS_ );
567 }
568 lp_->applyLeftPrec(*S,*leftPrecVec);
569 lp_->applyRightPrec(*leftPrecVec,*Z);
570 }
571 else {
572 lp_->applyLeftPrec(*S,*Z);
573 }
574 }
575 else if(lp_->isRightPrec()) {
576 lp_->applyRightPrec(*S,*Z);
577 }
578
579 // t = Az
580 lp_->applyOp(*Z,*T);
581
582 // omega = <K1\t,K1\s> / <K1\t,K1\t>
583 if(lp_->isLeftPrec()) {
584 if(leftPrecVec == Teuchos::null) {
585 leftPrecVec = MVT::Clone( *R_, numRHS_ );
586 }
587 if(leftPrecVec2 == Teuchos::null) {
588 leftPrecVec2 = MVT::Clone( *R_, numRHS_ );
589 }
590 lp_->applyLeftPrec(*T,*leftPrecVec2);
591 MVT::MvDot(*leftPrecVec2,*leftPrecVec2,tT);
592 MVT::MvDot(*leftPrecVec,*leftPrecVec2,tS);
593 }
594 else {
595 MVT::MvDot(*T,*T,tT);
596 MVT::MvDot(*T,*S,tS);
597 }
598 for(i=0; i<numRHS_; i++) {
599 if (SCT::magnitude(tT[i]) < MT::sfmin())
600 {
601 omega_[i] = SCT::zero();
602 breakdown_ = true;
603 }
604 else
605 omega_[i] = tS[i] / tT[i];
606 }
607
608 // x = x + alpha y + omega z
609 axpy(one, *X, alpha_, *Y, *X); // x = x + alpha y
610 axpy(one, *X, omega_, *Z, *X); // x = x + alpha y + omega z
611
612 // r = s - omega t
613 axpy(one, *S, omega_, *T, *R_, true);
614
615 // Update rho_old
616 rho_old_ = rho_new;
617 } // end while (sTest_->checkStatus(this) != Passed)
618 }
619
620
622 // Iterate until the status test informs us we should stop.
623 template <class ScalarType, class MV, class OP>
624 void BiCGStabIter<ScalarType,MV,OP>::axpy(const ScalarType alpha, const MV & A,
625 const std::vector<ScalarType> beta, const MV& B, MV& mv, bool minus)
626 {
627 Teuchos::RCP<const MV> A1, B1;
628 Teuchos::RCP<MV> mv1;
629 std::vector<int> index(1);
630
631 for(int i=0; i<numRHS_; i++) {
632 index[0] = i;
633 A1 = MVT::CloneView(A,index);
634 B1 = MVT::CloneView(B,index);
635 mv1 = MVT::CloneViewNonConst(mv,index);
636 if(minus) {
637 MVT::MvAddMv(alpha,*A1,-beta[i],*B1,*mv1);
638 }
639 else {
640 MVT::MvAddMv(alpha,*A1,beta[i],*B1,*mv1);
641 }
642 }
643 }
644
645} // end Belos namespace
646
647#endif /* BELOS_BICGSTAB_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.
This class implements the pseudo-block BiCGStab iteration, where the basic BiCGStab algorithm is perf...
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::ScalarTraits< ScalarType > SCT
bool breakdownDetected()
Has breakdown been detected in any linear system.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void resetNumIters(int iter=0)
Reset the iteration count.
BiCGStabIter(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)
BiCGStabIter constructor with linear problem, solver utilities, and parameter list of solver options.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
MultiVecTraits< ScalarType, MV > MVT
virtual ~BiCGStabIter()
Destructor.
void iterate()
This method performs BiCGStab iterations on each linear system until the status test indicates the ne...
int getNumIters() const
Get the current iteration count.
void initializeBiCGStab(BiCGStabIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
BiCGStabIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::ScalarTraits< MagnitudeType > MT
SCT::magnitudeType MagnitudeType
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
OperatorTraits< ScalarType, MV, OP > OPT
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...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to BiCGStabIteration state variables.
std::vector< ScalarType > omega
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Rhat
The initial residual.
std::vector< ScalarType > rho_old
Teuchos::RCP< const MV > P
The first decent direction vector.
std::vector< ScalarType > alpha
Teuchos::RCP< const MV > V
A * M * the first decent direction vector.

Generated for Belos by doxygen 1.10.0