Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosMinresIter.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_MINRES_ITER_HPP
11 #define BELOS_MINRES_ITER_HPP
12 
29 
30 #include "BelosConfigDefs.hpp"
31 #include "BelosTypes.hpp"
32 #include "BelosMinresIteration.hpp"
33 
34 #include "BelosLinearProblem.hpp"
35 #include "BelosOutputManager.hpp"
36 #include "BelosStatusTest.hpp"
37 #include "BelosOperatorTraits.hpp"
38 #include "BelosMultiVecTraits.hpp"
39 
42 #include "Teuchos_ScalarTraits.hpp"
44 #include "Teuchos_TimeMonitor.hpp"
45 
46 namespace Belos {
47 
61 template<class ScalarType, class MV, class OP>
62 class MinresIter : virtual public MinresIteration<ScalarType,MV,OP> {
63 
64  public:
65 
66  //
67  // Convenience typedefs
68  //
74 
76 
77 
87  const Teuchos::RCP< OutputManager< ScalarType > > & printer,
89  const Teuchos::ParameterList& params);
90 
92  virtual ~MinresIter() {};
94 
95 
97 
98 
113  void iterate();
114 
130 
136  void initialize()
137  {
139  initializeMinres(empty);
140  }
141 
149  if (! isInitialized())
150  throw std::logic_error("getState() cannot be called unless "
151  "the state has been initialized");
153  state.Y = Y_;
154  state.R1 = R1_;
155  state.R2 = R2_;
156  state.W = W_;
157  state.W1 = W1_;
158  state.W2 = W2_;
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 
177  getNativeResiduals( std::vector<MagnitudeType> *norms ) const
178  {
179  if (norms != NULL)
180  {
181  std::vector<MagnitudeType>& theNorms = *norms;
182  if (theNorms.size() < 1)
183  theNorms.resize(1);
184  theNorms[0] = phibar_;
185  }
186  return Teuchos::null;
187  }
188 
190 
192  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
193 
195  void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
196 
198 
200 
201 
203  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
204 
206  int getBlockSize() const { return 1; }
207 
209  void setBlockSize(int blockSize) {
210  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
211  "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
212  }
213 
215  bool isInitialized() const { return initialized_; }
216  bool isInitialized() { return initialized_; }
217 
219 
220  private:
221 
222  //
223  // Internal methods
224  //
226  void setStateSize();
227 
228  //
229  // Classes inputed through constructor that define the linear problem to be solved.
230  //
234 
235 
243  bool initialized_;
244 
251  bool stateStorageInitialized_;
252 
254  int iter_;
255 
260  MagnitudeType phibar_;
261 
262  //
263  // State Storage
264  //
265 
269  Teuchos::RCP< MV > R1_;
271  Teuchos::RCP< MV > R2_;
275  Teuchos::RCP< MV > W1_;
277  Teuchos::RCP< MV > W2_;
278 
287 
288 };
289 
291  // Constructor.
292  template<class ScalarType, class MV, class OP>
294  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
296  const Teuchos::ParameterList &/* params */ ):
297  lp_(problem),
298  om_(printer),
299  stest_(tester),
300  initialized_(false),
301  stateStorageInitialized_(false),
302  iter_(0),
303  phibar_(0.0)
304  {
305  }
306 
308  // Setup the state storage.
309  template <class ScalarType, class MV, class OP>
311  {
312  if (!stateStorageInitialized_) {
313 
314  // Check if there is any multivector to clone from.
315  Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
316  Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
317  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
318  stateStorageInitialized_ = false;
319  return;
320  }
321  else {
322 
323  // Initialize the state storage
324  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
325  if (Y_ == Teuchos::null) {
326  // Get the multivector that is not null.
327  Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
328  TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
329  std::invalid_argument,
330  "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
331  Y_ = MVT::Clone( *tmp, 1 );
332  R1_ = MVT::Clone( *tmp, 1 );
333  R2_ = MVT::Clone( *tmp, 1 );
334  W_ = MVT::Clone( *tmp, 1 );
335  W1_ = MVT::Clone( *tmp, 1 );
336  W2_ = MVT::Clone( *tmp, 1 );
337  }
338  // State storage has now been initialized.
339  stateStorageInitialized_ = true;
340  }
341  }
342  }
343 
344 
346  // Initialize this iteration object
347  template <class ScalarType, class MV, class OP>
349  {
350  // Initialize the state storage if it isn't already.
351  if (!stateStorageInitialized_)
352  setStateSize();
353 
354  TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
355  std::invalid_argument,
356  "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
357 
358  TEUCHOS_TEST_FOR_EXCEPTION( newstate.Y == Teuchos::null,
359  std::invalid_argument,
360  "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
361 
362  std::string errstr("Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
363  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.Y) != MVT::GetGlobalLength(*Y_),
364  std::invalid_argument,
365  errstr );
366  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.Y) != 1,
367  std::invalid_argument,
368  errstr );
369 
370  // Create convenience variables for zero, one.
371  const ScalarType one = SCT::one();
372  const MagnitudeType m_zero = SMT::zero();
373 
374  // Set up y and v for the first Lanczos vector v_1.
375  // y = beta1_ P' v1, where P = C**(-1).
376  // v is really P' v1.
377  MVT::Assign( *newstate.Y, *R2_ );
378  MVT::Assign( *newstate.Y, *R1_ );
379 
380  // Initialize the W's to 0.
381  MVT::MvInit ( *W_ );
382  MVT::MvInit ( *W2_ );
383 
384  if ( lp_->getLeftPrec() != Teuchos::null ) {
385  lp_->applyLeftPrec( *newstate.Y, *Y_ );
386  if ( lp_->getRightPrec() != Teuchos::null ) {
387  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
388  lp_->applyRightPrec( *tmp, *Y_ );
389  }
390  }
391  else if ( lp_->getRightPrec() != Teuchos::null ) {
392  lp_->applyRightPrec( *newstate.Y, *Y_ );
393  }
394  else {
395  if (newstate.Y != Y_) {
396  // copy over the initial residual (unpreconditioned).
397  MVT::Assign( *newstate.Y, *Y_ );
398  }
399  }
400 
401  // beta1_ = b'*y;
403  MVT::MvTransMv( one, *newstate.Y, *Y_, beta1_ );
404 
405  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < m_zero,
406  std::invalid_argument,
407  "The preconditioner is not positive definite." );
408 
409  if( SCT::magnitude(beta1_(0,0)) == m_zero )
410  {
411  // X = 0
412  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
413  MVT::MvInit( *cur_soln_vec );
414  }
415 
416  beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
417 
418  // The solver is initialized
419  initialized_ = true;
420  }
421 
422 
424  // Iterate until the status test informs us we should stop.
425  template <class ScalarType, class MV, class OP>
427  {
428  //
429  // Allocate/initialize data structures
430  //
431  if (initialized_ == false) {
432  initialize();
433  }
434 
435  // Create convenience variables for zero and one.
436  const ScalarType one = SCT::one();
437  const ScalarType zero = SCT::zero();
438  const MagnitudeType m_zero = SMT::zero();
439 
440  // Allocate memory for scalars.
443  phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
444 
445  // Initialize a few variables.
446  ScalarType oldBeta = zero;
447  ScalarType epsln = zero;
448  ScalarType cs = -one;
449  ScalarType sn = zero;
450  ScalarType dbar = zero;
451 
452  // Declare a few others that will be initialized in the loop.
453  ScalarType oldeps;
454  ScalarType delta;
455  ScalarType gbar;
456  ScalarType phi;
457  ScalarType gamma;
458 
459  // Allocate workspace.
460  Teuchos::RCP<MV> V = MVT::Clone( *Y_, 1 );
461  Teuchos::RCP<MV> tmpY, tmpW; // Not allocated, just used to transfer ownership.
462 
463  // Get the current solution vector.
464  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
465 
466  // Check that the current solution vector only has one column.
467  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
469  "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
470 
472  // Iterate until the status test tells us to stop.
473  //
474  while (stest_->checkStatus(this) != Passed) {
475 
476  // Increment the iteration
477  iter_++;
478 
479  // Normalize previous vector.
480  // v = y / beta(0,0);
481  MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
482 
483  // Apply operator.
484  lp_->applyOp (*V, *Y_);
485 
486  if (iter_ > 1)
487  MVT::MvAddMv (one, *Y_, -beta(0,0)/oldBeta, *R1_, *Y_);
488 
489  // alpha := dot(V, Y_)
490  MVT::MvTransMv (one, *V, *Y_, alpha);
491 
492  // y := y - alpha/beta r2
493  MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
494 
495  // r1 = r2;
496  // r2 = y;
497  tmpY = R1_;
498  R1_ = R2_;
499  R2_ = Y_;
500  Y_ = tmpY;
501 
502  // apply preconditioner
503  if ( lp_->getLeftPrec() != Teuchos::null ) {
504  lp_->applyLeftPrec( *R2_, *Y_ );
505  if ( lp_->getRightPrec() != Teuchos::null ) {
506  Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
507  lp_->applyRightPrec( *tmp, *Y_ );
508  }
509  }
510  else if ( lp_->getRightPrec() != Teuchos::null ) {
511  lp_->applyRightPrec( *R2_, *Y_ );
512  } // else "y = r2"
513  else {
514  MVT::Assign( *R2_, *Y_ );
515  }
516 
517  // Get new beta.
518  oldBeta = beta(0,0);
519  MVT::MvTransMv( one, *R2_, *Y_, beta );
520 
521  // Intercept beta <= 0.
522  //
523  // Note: we don't try to test for nonzero imaginary component of
524  // beta, because (a) it could be small and nonzero due to
525  // rounding error in computing the inner product, and (b) it's
526  // hard to tell how big "not small" should be, without computing
527  // some error bounds (for example, by modifying the linear
528  // algebra library to compute a posteriori rounding error bounds
529  // for the inner product, and then changing
530  // Belos::MultiVecTraits to make this information available).
531  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) < m_zero,
533  "Belos::MinresIter::iterate(): Encountered negative "
534  "value " << beta(0,0) << " for r2^H*M*r2 at itera"
535  "tion " << iter_ << ": MINRES cannot continue." );
536  beta(0,0) = SCT::squareroot( beta(0,0) );
537 
538  // Apply previous rotation Q_{k-1} to get
539  //
540  // [delta_k epsln_{k+1}] = [cs sn][dbar_k 0 ]
541  // [gbar_k dbar_{k+1} ] [-sn cs][alpha_k beta_{k+1}].
542  //
543  oldeps = epsln;
544  delta = cs*dbar + sn*alpha(0,0);
545  gbar = sn*dbar - cs*alpha(0,0);
546  epsln = sn*beta(0,0);
547  dbar = - cs*beta(0,0);
548 
549  // Compute the next plane rotation Q_k.
550  this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
551 
552  phi = cs * phibar_; // phi_k
553  phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ ); // phibar_{k+1}
554 
555  // w1 = w2;
556  // w2 = w;
557  MVT::Assign( *W_, *W1_ );
558  tmpW = W1_;
559  W1_ = W2_;
560  W2_ = W_;
561  W_ = tmpW;
562 
563  // w = (v - oldeps*w1 - delta*w2) / gamma;
564  MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
565  MVT::MvAddMv( one, *W_, -delta, *W2_, *W_ );
566  MVT::MvScale( *W_, one / gamma );
567 
568  // Update x:
569  // x = x + phi*w;
570  MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
571  lp_->updateSolution();
572  } // end while (sTest_->checkStatus(this) != Passed)
573  }
574 
575 
577  // Compute the next plane rotation Qk.
578  // r = norm([a b]);
579  // c = a / r;
580  // s = b / r;
581  template <class ScalarType, class MV, class OP>
582  void MinresIter<ScalarType,MV,OP>::symOrtho( ScalarType a, ScalarType b,
583  ScalarType *c, ScalarType *s, ScalarType *r
584  )
585  {
586  const ScalarType one = SCT::one();
587  const ScalarType zero = SCT::zero();
588  const MagnitudeType m_zero = SMT::zero();
589  const MagnitudeType absA = SCT::magnitude( a );
590  const MagnitudeType absB = SCT::magnitude( b );
591  if ( absB == m_zero ) {
592  *s = zero;
593  *r = absA;
594  if ( absA == m_zero )
595  *c = one;
596  else
597  *c = a / absA;
598  } else if ( absA == m_zero ) {
599  *c = zero;
600  *s = b / absB;
601  *r = absB;
602  } else if ( absB >= absA ) { // && a!=0 && b!=0
603  ScalarType tau = a / b;
604  if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
605  *s = -one / SCT::squareroot( one+tau*tau );
606  else
607  *s = one / SCT::squareroot( one+tau*tau );
608  *c = *s * tau;
609  *r = b / *s;
610  } else { // (absA > absB) && a!=0 && b!=0
611  ScalarType tau = b / a;
612  if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
613  *c = -one / SCT::squareroot( one+tau*tau );
614  else
615  *c = one / SCT::squareroot( one+tau*tau );
616  *s = *c * tau;
617  *r = a / *c;
618  }
619  }
620 
621 } // end Belos namespace
622 
623 #endif /* BELOS_MINRES_ITER_HPP */
Teuchos::RCP< const MV > R1
Previous residual.
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.
MINRES implementation.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void resetNumIters(int iter=0)
Reset the iteration count.
Declaration of basic traits for the multivector type.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
A pure virtual class for defining the status tests for the Belos iterative solvers.
SCT::magnitudeType MagnitudeType
Class which defines basic traits for the operator type.
virtual ~MinresIter()
Destructor.
MinresIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::ParameterList &params)
Constructor.
Traits class which defines basic operations on multivectors.
bool isInitialized()
States whether the solver has been initialized or not.
void iterate()
Perform MINRES iterations until convergence or error.
Teuchos::RCP< const MV > R2
Previous residual.
MinresIterateFailure is thrown when the MinresIteration object is unable to compute the next iterate ...
void symOrtho(ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Structure to contain pointers to MinresIteration state variables.
Teuchos::ScalarTraits< MagnitudeType > SMT
MultiVecTraits< ScalarType, MV > MVT
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > W2
Previous direction vector.
static magnitudeType magnitude(T a)
Pure virtual base class which augments the basic interface for a minimal residual linear solver itera...
OperatorTraits< ScalarType, MV, OP > OPT
void initialize()
Initialize the solver.
Teuchos::RCP< const MV > Y
The current residual.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initializeMinres(const MinresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
Class which defines basic traits for the operator type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
bool isInitialized() const
States whether the solver has been initialized or not.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > W1
Previous direction vector.
Teuchos::RCP< const MV > W
The current direction vector.

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