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 //
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_MINRES_ITER_HPP
43 #define BELOS_MINRES_ITER_HPP
44 
61 
62 #include "BelosConfigDefs.hpp"
63 #include "BelosTypes.hpp"
64 #include "BelosMinresIteration.hpp"
65 
66 #include "BelosLinearProblem.hpp"
67 #include "BelosOutputManager.hpp"
68 #include "BelosStatusTest.hpp"
69 #include "BelosOperatorTraits.hpp"
70 #include "BelosMultiVecTraits.hpp"
71 
74 #include "Teuchos_ScalarTraits.hpp"
76 #include "Teuchos_TimeMonitor.hpp"
77 #include "Teuchos_BLAS.hpp"
78 
79 namespace Belos {
80 
94 template<class ScalarType, class MV, class OP>
95 class MinresIter : virtual public MinresIteration<ScalarType,MV,OP> {
96 
97  public:
98 
99  //
100  // Convenience typedefs
101  //
107 
109 
110 
120  const Teuchos::RCP< OutputManager< ScalarType > > & printer,
122  const Teuchos::ParameterList& params);
123 
125  virtual ~MinresIter() {};
127 
128 
130 
131 
146  void iterate();
147 
163 
169  void initialize()
170  {
172  initializeMinres(empty);
173  }
174 
182  if (! isInitialized())
183  throw std::logic_error("getState() cannot be called unless "
184  "the state has been initialized");
186  state.Y = Y_;
187  state.R1 = R1_;
188  state.R2 = R2_;
189  state.W = W_;
190  state.W1 = W1_;
191  state.W2 = W2_;
192  return state;
193  }
194 
196 
197 
199 
200 
202  int getNumIters() const { return iter_; }
203 
205  void resetNumIters( int iter = 0 ) { iter_ = iter; }
206 
210  getNativeResiduals( std::vector<MagnitudeType> *norms ) const
211  {
212  if (norms != NULL)
213  {
214  std::vector<MagnitudeType>& theNorms = *norms;
215  if (theNorms.size() < 1)
216  theNorms.resize(1);
217  theNorms[0] = phibar_;
218  }
219  return Teuchos::null;
220  }
221 
223 
225  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
226 
228  void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
229 
231 
233 
234 
236  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
237 
239  int getBlockSize() const { return 1; }
240 
242  void setBlockSize(int blockSize) {
243  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
244  "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
245  }
246 
248  bool isInitialized() const { return initialized_; }
249  bool isInitialized() { return initialized_; }
250 
252 
253  private:
254 
255  //
256  // Internal methods
257  //
259  void setStateSize();
260 
261  //
262  // Classes inputed through constructor that define the linear problem to be solved.
263  //
267 
268 
276  bool initialized_;
277 
284  bool stateStorageInitialized_;
285 
287  int iter_;
288 
293  MagnitudeType phibar_;
294 
295  //
296  // State Storage
297  //
298 
302  Teuchos::RCP< MV > R1_;
304  Teuchos::RCP< MV > R2_;
308  Teuchos::RCP< MV > W1_;
310  Teuchos::RCP< MV > W2_;
311 
320 
321 };
322 
324  // Constructor.
325  template<class ScalarType, class MV, class OP>
327  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
329  const Teuchos::ParameterList &/* params */ ):
330  lp_(problem),
331  om_(printer),
332  stest_(tester),
333  initialized_(false),
334  stateStorageInitialized_(false),
335  iter_(0),
336  phibar_(0.0)
337  {
338  }
339 
341  // Setup the state storage.
342  template <class ScalarType, class MV, class OP>
344  {
345  if (!stateStorageInitialized_) {
346 
347  // Check if there is any multivector to clone from.
348  Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
349  Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
350  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
351  stateStorageInitialized_ = false;
352  return;
353  }
354  else {
355 
356  // Initialize the state storage
357  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
358  if (Y_ == Teuchos::null) {
359  // Get the multivector that is not null.
360  Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
361  TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
362  std::invalid_argument,
363  "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
364  Y_ = MVT::Clone( *tmp, 1 );
365  R1_ = MVT::Clone( *tmp, 1 );
366  R2_ = MVT::Clone( *tmp, 1 );
367  W_ = MVT::Clone( *tmp, 1 );
368  W1_ = MVT::Clone( *tmp, 1 );
369  W2_ = MVT::Clone( *tmp, 1 );
370  }
371  // State storage has now been initialized.
372  stateStorageInitialized_ = true;
373  }
374  }
375  }
376 
377 
379  // Initialize this iteration object
380  template <class ScalarType, class MV, class OP>
382  {
383  // Initialize the state storage if it isn't already.
384  if (!stateStorageInitialized_)
385  setStateSize();
386 
387  TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
388  std::invalid_argument,
389  "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
390 
391  TEUCHOS_TEST_FOR_EXCEPTION( newstate.Y == Teuchos::null,
392  std::invalid_argument,
393  "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
394 
395  std::string errstr("Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
396  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.Y) != MVT::GetGlobalLength(*Y_),
397  std::invalid_argument,
398  errstr );
399  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.Y) != 1,
400  std::invalid_argument,
401  errstr );
402 
403  // Create convenience variables for zero, one.
404  const ScalarType one = SCT::one();
405  const MagnitudeType zero = SMT::zero();
406 
407  // Set up y and v for the first Lanczos vector v_1.
408  // y = beta1_ P' v1, where P = C**(-1).
409  // v is really P' v1.
410  MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R2_ );
411  MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *R1_ );
412 
413  // Initialize the W's to 0.
414  MVT::MvInit ( *W_ );
415  MVT::MvInit ( *W2_ );
416 
417  if ( lp_->getLeftPrec() != Teuchos::null ) {
418  lp_->applyLeftPrec( *newstate.Y, *Y_ );
419  }
420  else {
421  if (newstate.Y != Y_) {
422  // copy over the initial residual (unpreconditioned).
423  MVT::MvAddMv( one, *newstate.Y, zero, *newstate.Y, *Y_ );
424  }
425  }
426 
427  // beta1_ = b'*y;
429  MVT::MvTransMv( one, *newstate.Y, *Y_, beta1_ );
430 
431  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < zero,
432  std::invalid_argument,
433  "The preconditioner is not positive definite." );
434 
435  if( SCT::magnitude(beta1_(0,0)) == zero )
436  {
437  // X = 0
438  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
439  MVT::MvInit( *cur_soln_vec );
440  }
441 
442  beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
443 
444  // The solver is initialized
445  initialized_ = true;
446  }
447 
448 
450  // Iterate until the status test informs us we should stop.
451  template <class ScalarType, class MV, class OP>
453  {
454  //
455  // Allocate/initialize data structures
456  //
457  if (initialized_ == false) {
458  initialize();
459  }
460 
462 
463  // Create convenience variables for zero and one.
464  const ScalarType one = SCT::one();
465  const MagnitudeType zero = SMT::zero();
466 
467  // Allocate memory for scalars.
470  phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
471  ScalarType shift = zero; // TODO Allow for proper shift.
472 
473  // Initialize a few variables.
474  ScalarType oldBeta = zero;
475  ScalarType epsln = zero;
476  ScalarType cs = -one;
477  ScalarType sn = zero;
478  ScalarType dbar = zero;
479 
480  // Declare a few others that will be initialized in the loop.
481  ScalarType oldeps;
482  ScalarType delta;
483  ScalarType gbar;
484  ScalarType phi;
485  ScalarType gamma;
486 
487  // Allocate workspace.
488  Teuchos::RCP<MV> V = MVT::Clone( *Y_, 1 );
489  Teuchos::RCP<MV> tmpY, tmpW; // Not allocated, just used to transfer ownership.
490 
491  // Get the current solution vector.
492  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
493 
494  // Check that the current solution vector only has one column.
495  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
497  "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
498 
500  // Iterate until the status test tells us to stop.
501  //
502  while (stest_->checkStatus(this) != Passed) {
503 
504  // Increment the iteration
505  iter_++;
506 
507  // Normalize previous vector.
508  // v = y / beta(0,0);
509  MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
510 
511  // Apply operator.
512  lp_->applyOp (*V, *Y_);
513 
514  // Apply shift
515  if (shift != zero)
516  MVT::MvAddMv (one, *Y_, -shift, *V, *Y_);
517 
518  if (iter_ > 1)
519  MVT::MvAddMv (one, *Y_, -beta(0,0)/oldBeta, *R1_, *Y_);
520 
521  // alpha := dot(V, Y_)
522  MVT::MvTransMv (one, *V, *Y_, alpha);
523 
524  // y := y - alpha/beta r2
525  MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
526 
527  // r1 = r2;
528  // r2 = y;
529  tmpY = R1_;
530  R1_ = R2_;
531  R2_ = Y_;
532  Y_ = tmpY;
533 
534  // apply left preconditioner
535  if ( lp_->getLeftPrec() != Teuchos::null ) {
536  lp_->applyLeftPrec( *R2_, *Y_ );
537  } // else "y = r2"
538  else {
539  MVT::MvAddMv( one, *R2_, zero, *R2_, *Y_ );
540  }
541 
542  // Get new beta.
543  oldBeta = beta(0,0);
544  MVT::MvTransMv( one, *R2_, *Y_, beta );
545 
546  // Intercept beta <= 0.
547  //
548  // Note: we don't try to test for nonzero imaginary component of
549  // beta, because (a) it could be small and nonzero due to
550  // rounding error in computing the inner product, and (b) it's
551  // hard to tell how big "not small" should be, without computing
552  // some error bounds (for example, by modifying the linear
553  // algebra library to compute a posteriori rounding error bounds
554  // for the inner product, and then changing
555  // Belos::MultiVecTraits to make this information available).
556  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) <= zero,
558  "Belos::MinresIter::iterate(): Encountered nonpositi"
559  "ve value " << beta(0,0) << " for r2^H*M*r2 at itera"
560  "tion " << iter_ << ": MINRES cannot continue." );
561  beta(0,0) = SCT::squareroot( beta(0,0) );
562 
563  // Apply previous rotation Q_{k-1} to get
564  //
565  // [delta_k epsln_{k+1}] = [cs sn][dbar_k 0 ]
566  // [gbar_k dbar_{k+1} ] [-sn cs][alpha_k beta_{k+1}].
567  //
568  oldeps = epsln;
569  delta = cs*dbar + sn*alpha(0,0);
570  gbar = sn*dbar - cs*alpha(0,0);
571  epsln = sn*beta(0,0);
572  dbar = - cs*beta(0,0);
573 
574  // Compute the next plane rotation Q_k.
575  this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
576 
577  phi = cs * phibar_; // phi_k
578  phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ ); // phibar_{k+1}
579 
580  // w1 = w2;
581  // w2 = w;
582  MVT::MvAddMv( one, *W_, zero, *W_, *W1_ );
583  tmpW = W1_;
584  W1_ = W2_;
585  W2_ = W_;
586  W_ = tmpW;
587 
588  // w = (v - oldeps*w1 - delta*w2) / gamma;
589  MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
590  MVT::MvAddMv( one, *W_, -delta, *W2_, *W_ );
591  MVT::MvScale( *W_, one / gamma );
592 
593  // Update x:
594  // x = x + phi*w;
595  MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
596  lp_->updateSolution();
597  } // end while (sTest_->checkStatus(this) != Passed)
598  }
599 
600 
602  // Compute the next plane rotation Qk.
603  // r = norm([a b]);
604  // c = a / r;
605  // s = b / r;
606  template <class ScalarType, class MV, class OP>
607  void MinresIter<ScalarType,MV,OP>::symOrtho( ScalarType a, ScalarType b,
608  ScalarType *c, ScalarType *s, ScalarType *r
609  )
610  {
611  const ScalarType one = SCT::one();
612  const ScalarType zero = SCT::zero();
613  const MagnitudeType m_zero = SMT::zero();
614  const MagnitudeType absA = SCT::magnitude( a );
615  const MagnitudeType absB = SCT::magnitude( b );
616  if ( absB == m_zero ) {
617  *s = zero;
618  *r = absA;
619  if ( absA == m_zero )
620  *c = one;
621  else
622  *c = a / absA;
623  } else if ( absA == m_zero ) {
624  *c = zero;
625  *s = b / absB;
626  *r = absB;
627  } else if ( absB >= absA ) { // && a!=0 && b!=0
628  ScalarType tau = a / b;
629  if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
630  *s = -one / SCT::squareroot( one+tau*tau );
631  else
632  *s = one / SCT::squareroot( one+tau*tau );
633  *c = *s * tau;
634  *r = b / *s;
635  } else { // (absA > absB) && a!=0 && b!=0
636  ScalarType tau = b / a;
637  if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
638  *c = -one / SCT::squareroot( one+tau*tau );
639  else
640  *c = one / SCT::squareroot( one+tau*tau );
641  *s = *c * tau;
642  *r = a / *c;
643  }
644  }
645 
646 } // end Belos namespace
647 
648 #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 Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5