Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
55#include "BelosStatusTest.hpp"
58
59#include "Teuchos_SerialDenseMatrix.hpp"
60#include "Teuchos_SerialDenseVector.hpp"
61#include "Teuchos_ScalarTraits.hpp"
62#include "Teuchos_ParameterList.hpp"
63#include "Teuchos_TimeMonitor.hpp"
64
74namespace Belos {
75
76template<class ScalarType, class MV, class OP>
77class FixedPointIter : virtual public FixedPointIteration<ScalarType,MV,OP> {
78
79 public:
80
81 //
82 // Convenience typedefs
83 //
86 typedef Teuchos::ScalarTraits<ScalarType> SCT;
87 typedef typename SCT::magnitudeType MagnitudeType;
88
90
91
97 FixedPointIter( 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 ~FixedPointIter() {};
105
106
108
109
122 void iterate();
123
139
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
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 //
213 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
214 const Teuchos::RCP<OutputManager<ScalarType> > om_;
215 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
216
217 // Algorithmic parameters
218 //
219 // blockSize_ is the solver block size.
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()
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.
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,
255 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
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();
388 const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero(); // unused
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 */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a fixed point linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
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 preconditioned fixed point iteration.
Teuchos::ScalarTraits< ScalarType > SCT
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
SCT::magnitudeType MagnitudeType
const Teuchos::RCP< OutputManager< ScalarType > > om_
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
virtual ~FixedPointIter()
Destructor.
void setStateSize()
Method for initalizing the state storage needed by FixedPoint.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
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.
int getNumIters() const
Get the current iteration count.
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...
void initializeFixedPoint(FixedPointIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
OperatorTraits< ScalarType, MV, OP > OPT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void iterate()
This method performs Fixed Point iterations until the status test indicates the need to stop or an er...
MultiVecTraits< ScalarType, MV > MVT
FixedPointIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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 FixedPointIteration state variables.
Teuchos::RCP< const MV > R
The current residual.
Teuchos::RCP< const MV > Z
The current preconditioned residual.