Belos  Version of the Day
BelosCGIter.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_CG_ITER_HPP
43 #define BELOS_CG_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosCGIteration.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 
75 namespace Belos {
76 
77 template<class ScalarType, class MV, class OP>
78 class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
79 
80  public:
81 
82  //
83  // Convenience typedefs
84  //
89 
91 
92 
99  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
101  Teuchos::ParameterList &params );
102 
104  virtual ~CGIter() {};
106 
107 
109 
110 
123  void iterate();
124 
140 
144  void initialize()
145  {
147  initializeCG(empty);
148  }
149 
158  state.R = R_;
159  state.P = P_;
160  state.AP = AP_;
161  state.Z = Z_;
162  return state;
163  }
164 
166 
167 
169 
170 
172  int getNumIters() const { return iter_; }
173 
175  void resetNumIters( int iter = 0 ) { iter_ = iter; }
176 
179  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
180 
182 
184  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
185 
187 
189 
190 
192  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
193 
195  int getBlockSize() const { return 1; }
196 
198  void setBlockSize(int blockSize) {
199  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
200  "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
201  }
202 
204  bool isInitialized() { return initialized_; }
205 
207 
208  private:
209 
210  //
211  // Internal methods
212  //
214  void setStateSize();
215 
216  //
217  // Classes inputed through constructor that define the linear problem to be solved.
218  //
222 
223  //
224  // Current solver state
225  //
226  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
227  // is capable of running; _initialize is controlled by the initialize() member method
228  // For the implications of the state of initialized_, please see documentation for initialize()
229  bool initialized_;
230 
231  // stateStorageInitialized_ specifies that the state storage has been initialized.
232  // This initialization may be postponed if the linear problem was generated without
233  // the right-hand side or solution vectors.
234  bool stateStorageInitialized_;
235 
236  // Current number of iterations performed.
237  int iter_;
238 
239  //
240  // State Storage
241  //
242  // Residual
243  Teuchos::RCP<MV> R_;
244  //
245  // Preconditioned residual
246  Teuchos::RCP<MV> Z_;
247  //
248  // Direction vector
249  Teuchos::RCP<MV> P_;
250  //
251  // Operator applied to direction vector
252  Teuchos::RCP<MV> AP_;
253 
254 };
255 
257  // Constructor.
258  template<class ScalarType, class MV, class OP>
260  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
262  Teuchos::ParameterList &params ):
263  lp_(problem),
264  om_(printer),
265  stest_(tester),
266  initialized_(false),
267  stateStorageInitialized_(false),
268  iter_(0)
269  {
270  }
271 
273  // Setup the state storage.
274  template <class ScalarType, class MV, class OP>
276  {
277  if (!stateStorageInitialized_) {
278 
279  // Check if there is any multivector to clone from.
280  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
281  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
282  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
283  stateStorageInitialized_ = false;
284  return;
285  }
286  else {
287 
288  // Initialize the state storage
289  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
290  if (R_ == Teuchos::null) {
291  // Get the multivector that is not null.
292  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
293  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
294  "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
295  R_ = MVT::Clone( *tmp, 1 );
296  Z_ = MVT::Clone( *tmp, 1 );
297  P_ = MVT::Clone( *tmp, 1 );
298  AP_ = MVT::Clone( *tmp, 1 );
299  }
300 
301  // State storage has now been initialized.
302  stateStorageInitialized_ = true;
303  }
304  }
305  }
306 
307 
309  // Initialize this iteration object
310  template <class ScalarType, class MV, class OP>
312  {
313  // Initialize the state storage if it isn't already.
314  if (!stateStorageInitialized_)
315  setStateSize();
316 
317  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
318  "Belos::CGIter::initialize(): Cannot initialize state storage!");
319 
320  // NOTE: In CGIter R_, the initial residual, is required!!!
321  //
322  std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
323 
324  // Create convenience variables for zero and one.
325  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
326  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
327 
328  if (newstate.R != Teuchos::null) {
329 
330  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
331  std::invalid_argument, errstr );
332  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
333  std::invalid_argument, errstr );
334 
335  // Copy basis vectors from newstate into V
336  if (newstate.R != R_) {
337  // copy over the initial residual (unpreconditioned).
338  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
339  }
340 
341  // Compute initial direction vectors
342  // Initially, they are set to the preconditioned residuals
343  //
344  if ( lp_->getLeftPrec() != Teuchos::null ) {
345  lp_->applyLeftPrec( *R_, *Z_ );
346  if ( lp_->getRightPrec() != Teuchos::null ) {
347  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
348  lp_->applyRightPrec( *Z_, *tmp );
349  Z_ = tmp;
350  }
351  }
352  else if ( lp_->getRightPrec() != Teuchos::null ) {
353  lp_->applyRightPrec( *R_, *Z_ );
354  }
355  else {
356  Z_ = R_;
357  }
358  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
359  }
360  else {
361 
362  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
363  "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
364  }
365 
366  // The solver is initialized
367  initialized_ = true;
368  }
369 
370 
372  // Iterate until the status test informs us we should stop.
373  template <class ScalarType, class MV, class OP>
375  {
376  //
377  // Allocate/initialize data structures
378  //
379  if (initialized_ == false) {
380  initialize();
381  }
382 
383  // Allocate memory for scalars.
386  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
387 
388  // Create convenience variables for zero and one.
389  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
391 
392  // Get the current solution vector.
393  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
394 
395  // Check that the current solution vector only has one column.
396  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
397  "Belos::CGIter::iterate(): current linear system has more than one vector!" );
398 
399  // Compute first <r,z> a.k.a. rHz
400  MVT::MvTransMv( one, *R_, *Z_, rHz );
401 
403  // Iterate until the status test tells us to stop.
404  //
405  while (stest_->checkStatus(this) != Passed) {
406 
407  // Increment the iteration
408  iter_++;
409 
410  // Multiply the current direction vector by A and store in AP_
411  lp_->applyOp( *P_, *AP_ );
412 
413  // Compute alpha := <R_,Z_> / <P_,AP_>
414  MVT::MvTransMv( one, *P_, *AP_, pAp );
415  alpha(0,0) = rHz(0,0) / pAp(0,0);
416 
417  // Check that alpha is a positive number!
418  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
419  "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
420  //
421  // Update the solution vector x := x + alpha * P_
422  //
423  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
424  lp_->updateSolution();
425  //
426  // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
427  //
428  rHz_old(0,0) = rHz(0,0);
429  //
430  // Compute the new residual R_ := R_ - alpha * AP_
431  //
432  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
433  //
434  // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ],
435  // and the new direction vector p.
436  //
437  if ( lp_->getLeftPrec() != Teuchos::null ) {
438  lp_->applyLeftPrec( *R_, *Z_ );
439  if ( lp_->getRightPrec() != Teuchos::null ) {
440  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
441  lp_->applyRightPrec( *Z_, *tmp );
442  Z_ = tmp;
443  }
444  }
445  else if ( lp_->getRightPrec() != Teuchos::null ) {
446  lp_->applyRightPrec( *R_, *Z_ );
447  }
448  else {
449  Z_ = R_;
450  }
451  //
452  MVT::MvTransMv( one, *R_, *Z_, rHz );
453  //
454  beta(0,0) = rHz(0,0) / rHz_old(0,0);
455  //
456  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
457 
458  } // end while (sTest_->checkStatus(this) != Passed)
459  }
460 
461 } // end Belos namespace
462 
463 #endif /* BELOS_CG_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::ScalarTraits< ScalarType > SCT
Definition: BelosCGIter.hpp:87
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.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
Structure to contain pointers to CGIteration state variables.
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
Definition: BelosCGIter.hpp:85
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
CGIter(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)
CGIter constructor with linear problem, solver utilities, and parameter list of solver options...
CGIterateFailure is thrown when the CGIteration object is unable to compute the next iterate in the C...
Teuchos::RCP< const MV > AP
The matrix A applied to current decent direction vector.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
Definition: BelosCGIter.hpp:86
Traits class which defines basic operations on multivectors.
int getNumIters() const
Get the current iteration count.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
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::RCP< const MV > Z
The current preconditioned residual.
virtual ~CGIter()
Destructor.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
SCT::magnitudeType MagnitudeType
Definition: BelosCGIter.hpp:88
void resetNumIters(int iter=0)
Reset the iteration count.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > P
The current decent direction vector.
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.

Generated for Belos by doxygen 1.8.14