Belos  Version of the Day
BelosCGSingleRedIter.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_SINGLE_RED_ITER_HPP
43 #define BELOS_CG_SINGLE_RED_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 
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 
75 namespace Belos {
76 
77 template<class ScalarType, class MV, class OP>
78 class CGSingleRedIter : virtual public CGIteration<ScalarType,MV,OP> {
79 
80  public:
81 
82  //
83  // Convenience typedefs
84  //
87  typedef Teuchos::ScalarTraits<ScalarType> SCT;
88  typedef typename SCT::magnitudeType MagnitudeType;
89 
91 
92 
98  CGSingleRedIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
99  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
100  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
101  Teuchos::ParameterList &params );
102 
104  virtual ~CGSingleRedIter() {};
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::CGSingleRedIter::setBlockSize(): Cannot use a block size that is not one.");
201  }
202 
204  bool isInitialized() { return initialized_; }
205 
207  void setDoCondEst(bool val){/*ignored*/}
208 
210  Teuchos::ArrayView<MagnitudeType> getDiag() {
211  Teuchos::ArrayView<MagnitudeType> temp;
212  return temp;
213  }
214 
216  Teuchos::ArrayView<MagnitudeType> getOffDiag() {
217  Teuchos::ArrayView<MagnitudeType> temp;
218  return temp;
219  }
220 
222 
223  private:
224 
225  //
226  // Internal methods
227  //
229  void setStateSize();
230 
231  //
232  // Classes inputed through constructor that define the linear problem to be solved.
233  //
234  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
235  const Teuchos::RCP<OutputManager<ScalarType> > om_;
236  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
237 
238  //
239  // Current solver state
240  //
241  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
242  // is capable of running; _initialize is controlled by the initialize() member method
243  // For the implications of the state of initialized_, please see documentation for initialize()
244  bool initialized_;
245 
246  // stateStorageInitialized_ specifies that the state storage has been initialized.
247  // This initialization may be postponed if the linear problem was generated without
248  // the right-hand side or solution vectors.
249  bool stateStorageInitialized_;
250 
251  // Current number of iterations performed.
252  int iter_;
253 
254  //
255  // State Storage
256  //
257  // Residual
258  Teuchos::RCP<MV> R_;
259  //
260  // Preconditioned residual
261  Teuchos::RCP<MV> Z_;
262  //
263  // Operator applied to preconditioned residual
264  Teuchos::RCP<MV> AZ_;
265  //
266  // This is the additional storage needed for single-reduction CG (Chronopoulos/Gear)
267  // R_ views the first vector and AZ_ views the second vector
268  Teuchos::RCP<MV> S_;
269  //
270  // Direction vector
271  Teuchos::RCP<MV> P_;
272  //
273  // Operator applied to direction vector (S)
274  Teuchos::RCP<MV> AP_;
275 
276 };
277 
279  // Constructor.
280  template<class ScalarType, class MV, class OP>
282  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
283  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
284  Teuchos::ParameterList &params ):
285  lp_(problem),
286  om_(printer),
287  stest_(tester),
288  initialized_(false),
289  stateStorageInitialized_(false),
290  iter_(0)
291  {
292  }
293 
295  // Setup the state storage.
296  template <class ScalarType, class MV, class OP>
298  {
299  if (!stateStorageInitialized_) {
300 
301  // Check if there is any multivector to clone from.
302  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
303  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
304  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
305  stateStorageInitialized_ = false;
306  return;
307  }
308  else {
309 
310  // Initialize the state storage
311  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
312  if (R_ == Teuchos::null) {
313  // Get the multivector that is not null.
314  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
315  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
316  "Belos::CGSingleRedIter::setStateSize(): linear problem does not specify multivectors to clone from.");
317  S_ = MVT::Clone( *tmp, 2 );
318  Z_ = MVT::Clone( *tmp, 1 );
319  P_ = MVT::Clone( *tmp, 1 );
320  AP_ = MVT::Clone( *tmp, 1 );
321 
322  // R_ will view the first column of S_, AZ_ will view the second.
323  std::vector<int> index(1,0);
324  R_ = MVT::CloneViewNonConst( *S_, index );
325  index[0] = 1;
326  AZ_ = MVT::CloneViewNonConst( *S_, index );
327  }
328 
329  // State storage has now been initialized.
330  stateStorageInitialized_ = true;
331  }
332  }
333  }
334 
335 
337  // Initialize this iteration object
338  template <class ScalarType, class MV, class OP>
340  {
341  // Initialize the state storage if it isn't already.
342  if (!stateStorageInitialized_)
343  setStateSize();
344 
345  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
346  "Belos::CGSingleRedIter::initialize(): Cannot initialize state storage!");
347 
348  // NOTE: In CGSingleRedIter R_, the initial residual, is required!!!
349  //
350  std::string errstr("Belos::CGSingleRedIter::initialize(): Specified multivectors must have a consistent length and width.");
351 
352  if (newstate.R != Teuchos::null) {
353 
354  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
355  std::invalid_argument, errstr );
356  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
357  std::invalid_argument, errstr );
358 
359  // Copy basis vectors from newstate into V
360  if (newstate.R != R_) {
361  // copy over the initial residual (unpreconditioned).
362  MVT::Assign( *newstate.R, *R_ );
363  }
364 
365  // Compute initial direction vectors
366  // Initially, they are set to the preconditioned residuals
367  //
368  if ( lp_->getLeftPrec() != Teuchos::null ) {
369  lp_->applyLeftPrec( *R_, *Z_ );
370  if ( lp_->getRightPrec() != Teuchos::null ) {
371  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
372  lp_->applyRightPrec( *Z_, *tmp );
373  Z_ = tmp;
374  }
375  }
376  else if ( lp_->getRightPrec() != Teuchos::null ) {
377  lp_->applyRightPrec( *R_, *Z_ );
378  }
379  else {
380  Z_ = R_;
381  }
382  MVT::Assign( *Z_, *P_ );
383 
384  // Multiply the current preconditioned residual vector by A and store in AZ_
385  lp_->applyOp( *Z_, *AZ_ );
386 
387  // Logically, AP_ = AZ_
388  MVT::Assign( *AZ_, *AP_ );
389  }
390  else {
391 
392  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
393  "Belos::CGSingleRedIter::initialize(): CGIterationState does not have initial residual.");
394  }
395 
396  // The solver is initialized
397  initialized_ = true;
398  }
399 
400 
402  // Iterate until the status test informs us we should stop.
403  template <class ScalarType, class MV, class OP>
405  {
406  //
407  // Allocate/initialize data structures
408  //
409  if (initialized_ == false) {
410  initialize();
411  }
412 
413  // Allocate memory for scalars.
414  Teuchos::SerialDenseMatrix<int,ScalarType> sHz( 2, 1 );
415  ScalarType rHz, rHz_old, alpha, beta, delta;
416 
417  // Create convenience variables for zero and one.
418  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
419  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
420 
421  // Get the current solution vector.
422  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
423 
424  // Check that the current solution vector only has one column.
425  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
426  "Belos::CGSingleRedIter::iterate(): current linear system has more than one vector!" );
427 
428  // Compute first <s,z> a.k.a. <r,z> and <Az,z> combined
429  MVT::MvTransMv( one, *S_, *Z_, sHz );
430  rHz = sHz(0,0);
431  delta = sHz(1,0);
432  alpha = rHz / delta;
433 
434  // Check that alpha is a positive number!
435  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGIterateFailure,
436  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
437 
439  // Iterate until the status test tells us to stop.
440  //
441  while (1) {
442 
443  // Update the solution vector x := x + alpha * P_
444  //
445  MVT::MvAddMv( one, *cur_soln_vec, alpha, *P_, *cur_soln_vec );
446  lp_->updateSolution();
447  //
448  // Compute the new residual R_ := R_ - alpha * AP_
449  //
450  MVT::MvAddMv( one, *R_, -alpha, *AP_, *R_ );
451  //
452  // Check the status test, now that the solution and residual have been updated
453  //
454  if (stest_->checkStatus(this) == Passed) {
455  break;
456  }
457  // Increment the iteration
458  iter_++;
459  //
460  // Apply preconditioner to new residual to update Z_
461  //
462  if ( lp_->getLeftPrec() != Teuchos::null ) {
463  lp_->applyLeftPrec( *R_, *Z_ );
464  if ( lp_->getRightPrec() != Teuchos::null ) {
465  Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
466  lp_->applyRightPrec( *Z_, *tmp );
467  Z_ = tmp;
468  }
469  }
470  else if ( lp_->getRightPrec() != Teuchos::null ) {
471  lp_->applyRightPrec( *R_, *Z_ );
472  }
473  else {
474  Z_ = R_;
475  }
476  //
477  // Multiply the current preconditioned residual vector by A and store in AZ_
478  lp_->applyOp( *Z_, *AZ_ );
479  //
480  // Compute <S_,Z_> a.k.a. <R_,Z_> and <AZ_,Z_> combined
481  MVT::MvTransMv( one, *S_, *Z_, sHz );
482  //
483  // Update scalars.
484  rHz_old = rHz;
485  rHz = sHz(0,0);
486  delta = sHz(1,0);
487  //
488  beta = rHz / rHz_old;
489  alpha = rHz / (delta - (beta*rHz / alpha));
490  //
491  // Check that alpha is a positive number!
492  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha) <= zero, CGIterateFailure,
493  "Belos::CGSingleRedIter::iterate(): non-positive value for p^H*A*p encountered!" );
494  //
495  // Update the direction vector P_ := Z_ + beta * P_
496  //
497  MVT::MvAddMv( one, *Z_, beta, *P_, *P_ );
498  //
499  // Update AP_ through recurrence relation AP_ := AZ_ + beta * AP_
500  // NOTE: This increases the number of vector updates by 1, in exchange for
501  // reducing the collectives from 2 to 1.
502  //
503  MVT::MvAddMv( one, *AZ_, beta, *AP_, *AP_ );
504  //
505  } // end while (sTest_->checkStatus(this) != Passed)
506  }
507 
508 } // end Belos namespace
509 
510 #endif /* BELOS_CG_SINGLE_RED_ITER_HPP */
Teuchos::RCP< const MV > R
The current residual.
int getNumIters() const
Get the current iteration count.
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.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Declaration of basic traits for the multivector type.
void iterate()
This method performs CG iterations until the status test indicates the need to stop or an error occur...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
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.
MultiVecTraits< ScalarType, MV > MVT
Traits class which defines basic operations on multivectors.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::ArrayView< MagnitudeType > getOffDiag()
Gets the off-diagonal for condition estimation (NOT_IMPLEMENTED)
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
bool isInitialized()
States whether the solver has been initialized or not.
CGIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
virtual ~CGSingleRedIter()
Destructor.
SCT::magnitudeType MagnitudeType
void initializeCG(CGIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Z
The current preconditioned residual.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
CGSingleRedIter(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)
CGSingleRedIter constructor with linear problem, solver utilities, and parameter list of solver optio...
void setDoCondEst(bool val)
Sets whether or not to store the diagonal for condition estimation.
Class which defines basic traits for the operator type.
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Teuchos::ScalarTraits< ScalarType > SCT
Belos header file which uses auto-configuration information to include necessary C++ headers...
void resetNumIters(int iter=0)
Reset the iteration count.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > P
The current decent direction vector.
Teuchos::ArrayView< MagnitudeType > getDiag()
Gets the diagonal for condition estimation (NOT_IMPLEMENTED)

Generated for Belos by doxygen 1.8.14