Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosFixedPointIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_FIXEDPOINT_ITER_HPP
11 #define BELOS_FIXEDPOINT_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosOutputManager.hpp"
23 #include "BelosStatusTest.hpp"
24 #include "BelosOperatorTraits.hpp"
25 #include "BelosMultiVecTraits.hpp"
26 
29 #include "Teuchos_ScalarTraits.hpp"
31 #include "Teuchos_TimeMonitor.hpp"
32 
42 namespace Belos {
43 
44 template<class ScalarType, class MV, class OP>
45 class FixedPointIter : virtual public FixedPointIteration<ScalarType,MV,OP> {
46 
47  public:
48 
49  //
50  // Convenience typedefs
51  //
56 
58 
59 
66  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
68  Teuchos::ParameterList &params );
69 
71  virtual ~FixedPointIter() {};
73 
74 
76 
77 
90  void iterate();
91 
107 
111  void initialize()
112  {
114  initializeFixedPoint(empty);
115  }
116 
125  state.R = R_;
126  state.Z = Z_;
127  return state;
128  }
129 
131 
132 
134 
135 
137  int getNumIters() const { return iter_; }
138 
140  void resetNumIters( int iter = 0 ) { iter_ = iter; }
141 
144  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
145 
147 
149  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
150 
152 
154 
155 
157  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
158 
160  int getBlockSize() const { return numRHS_; }
161 
163  void setBlockSize(int blockSize);
164 
166  bool isInitialized() { return initialized_; }
167 
169 
170  private:
171 
172  //
173  // Internal methods
174  //
176  void setStateSize();
177 
178  //
179  // Classes inputed through constructor that define the linear problem to be solved.
180  //
184 
185  // Algorithmic parameters
186  //
187  // blockSize_ is the solver block size.
188  int numRHS_;
189 
190  //
191  // Current solver state
192  //
193  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
194  // is capable of running; _initialize is controlled by the initialize() member method
195  // For the implications of the state of initialized_, please see documentation for initialize()
196  bool initialized_;
197 
198  // stateStorageInitialized_ specifies that the state storage has been initialized.
199  // This initialization may be postponed if the linear problem was generated without
200  // the right-hand side or solution vectors.
201  bool stateStorageInitialized_;
202 
203  // Current number of iterations performed.
204  int iter_;
205 
206  //
207  // State Storage
208  //
209  // Residual
210  Teuchos::RCP<MV> R_;
211  //
212  // Preconditioned residual
213  Teuchos::RCP<MV> Z_;
214  //
215 
216 };
217 
219  // Constructor.
220  template<class ScalarType, class MV, class OP>
222  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
224  Teuchos::ParameterList &params ):
225  lp_(problem),
226  om_(printer),
227  stest_(tester),
228  numRHS_(0),
229  initialized_(false),
230  stateStorageInitialized_(false),
231  iter_(0)
232  {
233  setBlockSize(params.get("Block Size",MVT::GetNumberVecs(*problem->getCurrRHSVec())));
234  }
235 
237  // Setup the state storage.
238  template <class ScalarType, class MV, class OP>
240  {
241  if (!stateStorageInitialized_) {
242  // Check if there is any multivector to clone from.
243  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
244  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
245  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
246  stateStorageInitialized_ = false;
247  return;
248  }
249  else {
250 
251  // Initialize the state storage
252  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
253  if (R_ == Teuchos::null) {
254  // Get the multivector that is not null.
255  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
256  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
257  "Belos::FixedPointIter::setStateSize(): linear problem does not specify multivectors to clone from.");
258  R_ = MVT::Clone( *tmp, numRHS_ );
259  Z_ = MVT::Clone( *tmp, numRHS_ );
260  }
261 
262  // State storage has now been initialized.
263  stateStorageInitialized_ = true;
264  }
265  }
266  }
267 
269  // Set the block size and make necessary adjustments.
270  template <class ScalarType, class MV, class OP>
272  {
273  // This routine only allocates space; it doesn't not perform any computation
274  // any change in size will invalidate the state of the solver.
275 
276  TEUCHOS_TEST_FOR_EXCEPTION(blockSize != MVT::GetNumberVecs(*lp_->getCurrRHSVec()), std::invalid_argument, "Belos::FixedPointIter::setBlockSize size must match # RHS.");
277 
278  TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument, "Belos::FixedPointIter::setBlockSize was passed a non-positive argument.");
279  if (blockSize == numRHS_) {
280  // do nothing
281  return;
282  }
283 
284  if (blockSize!=numRHS_)
285  stateStorageInitialized_ = false;
286 
287  numRHS_ = blockSize;
288 
289  initialized_ = false;
290 
291  // Use the current blockSize_ to initialize the state storage.
292  setStateSize();
293 
294  }
295 
297  // Initialize this iteration object
298  template <class ScalarType, class MV, class OP>
300  {
301  // Initialize the state storage if it isn't already.
302  if (!stateStorageInitialized_)
303  setStateSize();
304 
305  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
306  "Belos::FixedPointIter::initialize(): Cannot initialize state storage!");
307 
308  // NOTE: In FixedPointIter R_, the initial residual, is required!!!
309  //
310  std::string errstr("Belos::FixedPointIter::initialize(): Specified multivectors must have a consistent length and width.");
311 
312  // Create convenience variables for zero and one.
313  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
314  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
315 
316  if (newstate.R != Teuchos::null) {
317  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*R_) != MVT::GetNumberVecs(*newstate.R),
318  std::invalid_argument, errstr );
319 
320  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
321  std::invalid_argument, errstr );
322  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
323  std::invalid_argument, errstr );
324 
325  // Copy basis vectors from newstate into V
326  if (newstate.R != R_) {
327  // copy over the initial residual (unpreconditioned).
328  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
329  }
330 
331  }
332  else {
333  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
334  "Belos::FixedPointIter::initialize(): FixedPointIterationState does not have initial residual.");
335  }
336 
337  // The solver is initialized
338  initialized_ = true;
339  }
340 
341 
343  // Iterate until the status test informs us we should stop.
344  template <class ScalarType, class MV, class OP>
346  {
347  //
348  // Allocate/initialize data structures
349  //
350  if (initialized_ == false) {
351  initialize();
352  }
353 
354  // Create convenience variables for zero and one.
355  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
357 
358  // Get the current solution vector.
359  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
360 
361  // Temp vector
362  Teuchos::RCP<MV> tmp = MVT::Clone( *R_, numRHS_ );
363 
364  if (lp_->getRightPrec() != Teuchos::null) {
365  // Set rhs to initial residual
366  Teuchos::RCP<MV> rhs = MVT::CloneCopy( *R_ );
367 
368  // Zero initial guess
369  MVT::MvInit( *Z_, zero );
370 
372  // Iterate until the status test tells us to stop.
373  //
374  while (stest_->checkStatus(this) != Passed) {
375 
376  // Increment the iteration
377  iter_++;
378 
379  // Apply preconditioner
380  lp_->applyRightPrec( *R_, *tmp );
381 
382  // Update solution vector
383  MVT::MvAddMv( one, *cur_soln_vec, one, *tmp, *cur_soln_vec );
384  lp_->updateSolution();
385 
386  // Update solution vector
387  MVT::MvAddMv( one, *Z_, one, *tmp, *Z_ );
388 
389  // Compute new residual
390  lp_->applyOp (*Z_, *tmp );
391  MVT::MvAddMv( one, *rhs, -one, *tmp, *R_ );
392 
393  } // end while (sTest_->checkStatus(this) != Passed)
394 
395  } else {
396  Teuchos::RCP<const MV> rhs = lp_->getCurrRHSVec();
397 
399  // Iterate until the status test tells us to stop.
400  //
401  while (stest_->checkStatus(this) != Passed) {
402 
403  // Increment the iteration
404  iter_++;
405 
406  // Compute initial preconditioned residual
407  if ( lp_->getLeftPrec() != Teuchos::null ) {
408  lp_->applyLeftPrec( *R_, *Z_ );
409  }
410  else {
411  Z_ = R_;
412  }
413 
414  // Update solution vector
415  MVT::MvAddMv(one,*cur_soln_vec,one,*Z_,*cur_soln_vec);
416  lp_->updateSolution();
417 
418  // Compute new residual
419  lp_->applyOp(*cur_soln_vec,*tmp);
420  MVT::MvAddMv(one,*rhs,-one,*tmp,*R_);
421 
422  } // end while (sTest_->checkStatus(this) != Passed)
423  }
424  }
425 
426 } // end Belos namespace
427 
428 #endif /* BELOS_FIXEDPOINT_ITER_HPP */
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.
MultiVecTraits< ScalarType, MV > MVT
T & get(const std::string &name, T def_value)
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
Pure virtual base class for defining the status testing capabilities of Belos.
virtual ~FixedPointIter()
Destructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::ScalarTraits< ScalarType > SCT
void initializeFixedPoint(FixedPointIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Declaration of basic traits for the multivector type.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void iterate()
This method performs Fixed Point iterations until the status test indicates the need to stop or an er...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
int getNumIters() const
Get the current iteration count.
Traits class which defines basic operations on multivectors.
FixedPointIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void resetNumIters(int iter=0)
Reset the iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
FixedPointIter(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)
FixedPointIter constructor with linear problem, solver utilities, and parameter list of solver option...
Teuchos::RCP< const MV > R
The current residual.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SCT::magnitudeType MagnitudeType
Structure to contain pointers to FixedPointIteration state variables.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Class which defines basic traits for the operator type.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
Pure virtual base class which augments the basic interface for a fixed point linear solver iteration...
Belos header file which uses auto-configuration information to include necessary C++ headers...
This class implements the preconditioned fixed point iteration.
OperatorTraits< ScalarType, MV, OP > OPT

Generated on Mon Jul 15 2024 09:24:23 for Belos by doxygen 1.8.5