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 //
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_FIXEDPOINT_ITER_HPP
43 #define BELOS_FIXEDPOINT_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
52 
53 #include "BelosLinearProblem.hpp"
54 #include "BelosOutputManager.hpp"
55 #include "BelosStatusTest.hpp"
56 #include "BelosOperatorTraits.hpp"
57 #include "BelosMultiVecTraits.hpp"
58 
61 #include "Teuchos_ScalarTraits.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
74 namespace Belos {
75 
76 template<class ScalarType, class MV, class OP>
77 class FixedPointIter : virtual public FixedPointIteration<ScalarType,MV,OP> {
78 
79  public:
80 
81  //
82  // Convenience typedefs
83  //
88 
90 
91 
98  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
100  Teuchos::ParameterList &params );
101 
103  virtual ~FixedPointIter() {};
105 
106 
108 
109 
122  void iterate();
123 
139 
143  void initialize()
144  {
146  initializeFixedPoint(empty);
147  }
148 
157  state.R = R_;
158  state.Z = Z_;
159  return state;
160  }
161 
163 
164 
166 
167 
169  int getNumIters() const { return iter_; }
170 
172  void resetNumIters( int iter = 0 ) { iter_ = iter; }
173 
176  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
177 
179 
181  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
182 
184 
186 
187 
189  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
190 
192  int getBlockSize() const { return numRHS_; }
193 
195  void setBlockSize(int blockSize);
196 
198  bool isInitialized() { return initialized_; }
199 
201 
202  private:
203 
204  //
205  // Internal methods
206  //
208  void setStateSize();
209 
210  //
211  // Classes inputed through constructor that define the linear problem to be solved.
212  //
216 
217  // Algorithmic parameters
218  //
219  // blockSize_ is the solver block size.
220  int numRHS_;
221 
222  //
223  // Current solver state
224  //
225  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
226  // is capable of running; _initialize is controlled by the initialize() member method
227  // For the implications of the state of initialized_, please see documentation for initialize()
228  bool initialized_;
229 
230  // stateStorageInitialized_ specifies that the state storage has been initialized.
231  // This initialization may be postponed if the linear problem was generated without
232  // the right-hand side or solution vectors.
233  bool stateStorageInitialized_;
234 
235  // Current number of iterations performed.
236  int iter_;
237 
238  //
239  // State Storage
240  //
241  // Residual
242  Teuchos::RCP<MV> R_;
243  //
244  // Preconditioned residual
245  Teuchos::RCP<MV> Z_;
246  //
247 
248 };
249 
251  // Constructor.
252  template<class ScalarType, class MV, class OP>
254  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
256  Teuchos::ParameterList &params ):
257  lp_(problem),
258  om_(printer),
259  stest_(tester),
260  numRHS_(0),
261  initialized_(false),
262  stateStorageInitialized_(false),
263  iter_(0)
264  {
265  setBlockSize(params.get("Block Size",MVT::GetNumberVecs(*problem->getCurrRHSVec())));
266  }
267 
269  // Setup the state storage.
270  template <class ScalarType, class MV, class OP>
272  {
273  if (!stateStorageInitialized_) {
274  // Check if there is any multivector to clone from.
275  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
276  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
277  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
278  stateStorageInitialized_ = false;
279  return;
280  }
281  else {
282 
283  // Initialize the state storage
284  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
285  if (R_ == Teuchos::null) {
286  // Get the multivector that is not null.
287  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
288  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
289  "Belos::FixedPointIter::setStateSize(): linear problem does not specify multivectors to clone from.");
290  R_ = MVT::Clone( *tmp, numRHS_ );
291  Z_ = MVT::Clone( *tmp, numRHS_ );
292  }
293 
294  // State storage has now been initialized.
295  stateStorageInitialized_ = true;
296  }
297  }
298  }
299 
301  // Set the block size and make necessary adjustments.
302  template <class ScalarType, class MV, class OP>
304  {
305  // This routine only allocates space; it doesn't not perform any computation
306  // any change in size will invalidate the state of the solver.
307 
308  TEUCHOS_TEST_FOR_EXCEPTION(blockSize != MVT::GetNumberVecs(*lp_->getCurrRHSVec()), std::invalid_argument, "Belos::FixedPointIter::setBlockSize size must match # RHS.");
309 
310  TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument, "Belos::FixedPointIter::setBlockSize was passed a non-positive argument.");
311  if (blockSize == numRHS_) {
312  // do nothing
313  return;
314  }
315 
316  if (blockSize!=numRHS_)
317  stateStorageInitialized_ = false;
318 
319  numRHS_ = blockSize;
320 
321  initialized_ = false;
322 
323  // Use the current blockSize_ to initialize the state storage.
324  setStateSize();
325 
326  }
327 
329  // Initialize this iteration object
330  template <class ScalarType, class MV, class OP>
332  {
333  // Initialize the state storage if it isn't already.
334  if (!stateStorageInitialized_)
335  setStateSize();
336 
337  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
338  "Belos::FixedPointIter::initialize(): Cannot initialize state storage!");
339 
340  // NOTE: In FixedPointIter R_, the initial residual, is required!!!
341  //
342  std::string errstr("Belos::FixedPointIter::initialize(): Specified multivectors must have a consistent length and width.");
343 
344  // Create convenience variables for zero and one.
345  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
346  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
347 
348  if (newstate.R != Teuchos::null) {
349  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*R_) != MVT::GetNumberVecs(*newstate.R),
350  std::invalid_argument, errstr );
351 
352  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
353  std::invalid_argument, errstr );
354  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
355  std::invalid_argument, errstr );
356 
357  // Copy basis vectors from newstate into V
358  if (newstate.R != R_) {
359  // copy over the initial residual (unpreconditioned).
360  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
361  }
362 
363  }
364  else {
365  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
366  "Belos::FixedPointIter::initialize(): FixedPointIterationState does not have initial residual.");
367  }
368 
369  // The solver is initialized
370  initialized_ = true;
371  }
372 
373 
375  // Iterate until the status test informs us we should stop.
376  template <class ScalarType, class MV, class OP>
378  {
379  //
380  // Allocate/initialize data structures
381  //
382  if (initialized_ == false) {
383  initialize();
384  }
385 
386  // Create convenience variables for zero and one.
387  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
389 
390  // Get the current solution vector.
391  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
392 
393  // Temp vector
394  Teuchos::RCP<MV> tmp = MVT::Clone( *R_, numRHS_ );
395 
396  if (lp_->getRightPrec() != Teuchos::null) {
397  // Set rhs to initial residual
398  Teuchos::RCP<MV> rhs = MVT::CloneCopy( *R_ );
399 
400  // Zero initial guess
401  MVT::MvInit( *Z_, zero );
402 
404  // Iterate until the status test tells us to stop.
405  //
406  while (stest_->checkStatus(this) != Passed) {
407 
408  // Increment the iteration
409  iter_++;
410 
411  // Apply preconditioner
412  lp_->applyRightPrec( *R_, *tmp );
413 
414  // Update solution vector
415  MVT::MvAddMv( one, *cur_soln_vec, one, *tmp, *cur_soln_vec );
416  lp_->updateSolution();
417 
418  // Update solution vector
419  MVT::MvAddMv( one, *Z_, one, *tmp, *Z_ );
420 
421  // Compute new residual
422  lp_->applyOp (*Z_, *tmp );
423  MVT::MvAddMv( one, *rhs, -one, *tmp, *R_ );
424 
425  } // end while (sTest_->checkStatus(this) != Passed)
426 
427  } else {
428  Teuchos::RCP<const MV> rhs = lp_->getCurrRHSVec();
429 
431  // Iterate until the status test tells us to stop.
432  //
433  while (stest_->checkStatus(this) != Passed) {
434 
435  // Increment the iteration
436  iter_++;
437 
438  // Compute initial preconditioned residual
439  if ( lp_->getLeftPrec() != Teuchos::null ) {
440  lp_->applyLeftPrec( *R_, *Z_ );
441  }
442  else {
443  Z_ = R_;
444  }
445 
446  // Update solution vector
447  MVT::MvAddMv(one,*cur_soln_vec,one,*Z_,*cur_soln_vec);
448  lp_->updateSolution();
449 
450  // Compute new residual
451  lp_->applyOp(*cur_soln_vec,*tmp);
452  MVT::MvAddMv(one,*rhs,-one,*tmp,*R_);
453 
454  } // end while (sTest_->checkStatus(this) != Passed)
455  }
456  }
457 
458 } // end Belos namespace
459 
460 #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 Thu Mar 28 2024 09:27:40 for Belos by doxygen 1.8.5