Belos  Version of the Day
BelosTFQMRIter.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 // This file contains an implementation of the TFQMR iteration
43 // for solving non-Hermitian linear systems of equations Ax = b,
44 // where b is a single-vector and x is the corresponding solution.
45 //
46 // The implementation is a slight modification on the TFQMR iteration
47 // found in Saad's "Iterative Methods for Sparse Linear Systems".
48 //
49 
50 #ifndef BELOS_TFQMR_ITER_HPP
51 #define BELOS_TFQMR_ITER_HPP
52 
60 #include "BelosConfigDefs.hpp"
61 #include "BelosIteration.hpp"
62 #include "BelosTypes.hpp"
63 
64 #include "BelosLinearProblem.hpp"
65 #include "BelosMatOrthoManager.hpp"
66 #include "BelosOutputManager.hpp"
67 #include "BelosStatusTest.hpp"
68 #include "BelosOperatorTraits.hpp"
69 #include "BelosMultiVecTraits.hpp"
70 
71 #include "Teuchos_BLAS.hpp"
74 #include "Teuchos_ScalarTraits.hpp"
76 #include "Teuchos_TimeMonitor.hpp"
77 
89 namespace Belos {
90 
95  template <class ScalarType, class MV>
96  struct TFQMRIterState {
97 
105 
106  TFQMRIterState() : R(Teuchos::null), W(Teuchos::null), U(Teuchos::null),
107  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
108  {}
109  };
110 
111 
113 
114 
126  class TFQMRIterInitFailure : public BelosError {public:
127  TFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
128  {}};
129 
136  class TFQMRIterateFailure : public BelosError {public:
137  TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
138  {}};
139 
141 
142 
143  template <class ScalarType, class MV, class OP>
144  class TFQMRIter : public Iteration<ScalarType,MV,OP> {
145  public:
146  //
147  // Convenience typedefs
148  //
153 
155 
156 
159  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
161  Teuchos::ParameterList &params );
162 
164  virtual ~TFQMRIter() {};
166 
167 
169 
170 
192  void iterate();
193 
215  void initializeTFQMR(const TFQMRIterState<ScalarType,MV> & newstate);
216 
220  void initialize()
221  {
223  initializeTFQMR(empty);
224  }
225 
235  state.R = R_;
236  state.W = W_;
237  state.U = U_;
238  state.Rtilde = Rtilde_;
239  state.D = D_;
240  state.V = V_;
241  state.solnUpdate = solnUpdate_;
242  return state;
243  }
244 
246 
247 
249 
250 
252  int getNumIters() const { return iter_; }
253 
255  void resetNumIters( int iter = 0 ) { iter_ = iter; }
256 
259  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
260 
262 
265  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
266 
268 
269 
271 
272 
274  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
275 
277  int getBlockSize() const { return 1; }
278 
280  void setBlockSize(int blockSize) {
281  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
282  "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
283  }
284 
286  bool isInitialized() { return initialized_; }
287 
289 
290 
291  private:
292 
293  //
294  // Internal methods
295  //
297  void setStateSize();
298 
299  //
300  // Classes inputed through constructor that define the linear problem to be solved.
301  //
305 
306  //
307  // Algorithmic parameters
308  //
309 
310  // Storage for QR factorization of the least squares system.
311  // Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
312  std::vector<ScalarType> alpha_, rho_, rho_old_;
313  std::vector<MagnitudeType> tau_, cs_, theta_;
314 
315  //
316  // Current solver state
317  //
318  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
319  // is capable of running; _initialize is controlled by the initialize() member method
320  // For the implications of the state of initialized_, please see documentation for initialize()
321  bool initialized_;
322 
323  // stateStorageInitialized_ specifies that the state storage has be initialized to the current
324  // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
325  // generated without the right-hand side or solution vectors.
326  bool stateStorageInitialized_;
327 
328  // Current subspace dimension, and number of iterations performed.
329  int iter_;
330 
331  //
332  // State Storage
333  //
334  Teuchos::RCP<MV> R_;
335  Teuchos::RCP<MV> W_;
336  Teuchos::RCP<MV> U_, AU_;
337  Teuchos::RCP<MV> Rtilde_;
338  Teuchos::RCP<MV> D_;
339  Teuchos::RCP<MV> V_;
340  Teuchos::RCP<MV> solnUpdate_;
341  };
342 
343 
344  //
345  // Implementation
346  //
347 
349  // Constructor.
350  template <class ScalarType, class MV, class OP>
352  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
354  Teuchos::ParameterList &params
355  ) :
356  lp_(problem),
357  om_(printer),
358  stest_(tester),
359  alpha_(1),
360  rho_(1),
361  rho_old_(1),
362  tau_(1),
363  cs_(1),
364  theta_(1),
365  initialized_(false),
366  stateStorageInitialized_(false),
367  iter_(0)
368  {
369  }
370 
372  // Compute native residual from TFQMR recurrence.
373  template <class ScalarType, class MV, class OP>
375  TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
376  {
378  if (normvec)
379  (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
380 
381  return Teuchos::null;
382  }
383 
384 
386  // Setup the state storage.
387  template <class ScalarType, class MV, class OP>
389  {
390  if (!stateStorageInitialized_) {
391 
392  // Check if there is any multivector to clone from.
393  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
394  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
395  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
396  stateStorageInitialized_ = false;
397  return;
398  }
399  else {
400 
401  // Initialize the state storage
402  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
403  if (R_ == Teuchos::null) {
404  // Get the multivector that is not null.
405  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
406  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
407  "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
408  R_ = MVT::Clone( *tmp, 1 );
409  D_ = MVT::Clone( *tmp, 1 );
410  V_ = MVT::Clone( *tmp, 1 );
411  solnUpdate_ = MVT::Clone( *tmp, 1 );
412  }
413 
414  // State storage has now been initialized.
415  stateStorageInitialized_ = true;
416  }
417  }
418  }
419 
421  // Initialize this iteration object
422  template <class ScalarType, class MV, class OP>
424  {
425  // Initialize the state storage if it isn't already.
426  if (!stateStorageInitialized_)
427  setStateSize();
428 
429  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
430  "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
431 
432  // NOTE: In TFQMRIter R_, the initial residual, is required!!!
433  //
434  std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
435 
436  // Create convenience variables for zero and one.
437  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
438  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
440 
441  if (newstate.R != Teuchos::null) {
442 
443  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
444  std::invalid_argument, errstr );
445  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
446  std::invalid_argument, errstr );
447 
448  // Copy basis vectors from newstate into V
449  if (newstate.R != R_) {
450  // copy over the initial residual (unpreconditioned).
451  MVT::MvAddMv( one, *newstate.R, STzero, *newstate.R, *R_ );
452  }
453 
454  // Compute initial vectors
455  // Initially, they are set to the preconditioned residuals
456  //
457  W_ = MVT::CloneCopy( *R_ );
458  U_ = MVT::CloneCopy( *R_ );
459  Rtilde_ = MVT::CloneCopy( *R_ );
460  MVT::MvInit( *D_ );
461  MVT::MvInit( *solnUpdate_ );
462  // Multiply the current residual by Op and store in V_
463  // V_ = Op * R_
464  //
465  lp_->apply( *U_, *V_ );
466  AU_ = MVT::CloneCopy( *V_ );
467  //
468  // Compute initial scalars: theta, eta, tau, rho_old
469  //
470  theta_[0] = MTzero;
471  MVT::MvNorm( *R_, tau_ ); // tau = ||r_0||
472  MVT::MvDot( *R_, *Rtilde_, rho_old_ ); // rho = (r_tilde, r0)
473  }
474  else {
475 
476  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
477  "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
478  }
479 
480  // The solver is initialized
481  initialized_ = true;
482  }
483 
484 
486  // Iterate until the status test informs us we should stop.
487  template <class ScalarType, class MV, class OP>
489  {
490  //
491  // Allocate/initialize data structures
492  //
493  if (initialized_ == false) {
494  initialize();
495  }
496 
497  // Create convenience variables for zero and one.
498  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
500  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
501  ScalarType eta = STzero, beta = STzero;
502  //
503  // Start executable statements.
504  //
505  // Get the current solution vector.
506  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
507 
508  // Check that the current solution vector only has one column.
509  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
510  "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
511 
512 
514  // Iterate until the status test tells us to stop.
515  //
516  while (stest_->checkStatus(this) != Passed) {
517 
518  for (int iIter=0; iIter<2; iIter++)
519  {
520  //
521  //--------------------------------------------------------
522  // Compute the new alpha if we need to
523  //--------------------------------------------------------
524  //
525  if (iIter == 0) {
526  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
527  alpha_[0] = rho_old_[0]/alpha_[0];
528  }
529  //
530  //--------------------------------------------------------
531  // Update w.
532  // w = w - alpha*Au
533  //--------------------------------------------------------
534  //
535  MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
536  //
537  //--------------------------------------------------------
538  // Update d.
539  // d = u + (theta^2/alpha)eta*d
540  //--------------------------------------------------------
541  //
542  MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
543  //
544  //--------------------------------------------------------
545  // Update u if we need to.
546  // u = u - alpha*v
547  //
548  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
549  //--------------------------------------------------------
550  //
551  if (iIter == 0) {
552  // Compute new U.
553  MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
554 
555  // Update Au for the next iteration.
556  lp_->apply( *U_, *AU_ );
557  }
558  //
559  //--------------------------------------------------------
560  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
561  //--------------------------------------------------------
562  //
563  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
564  theta_[0] /= tau_[0];
565  // cs = 1.0 / sqrt(1.0 + theta^2)
566  cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
567  tau_[0] *= theta_[0]*cs_[0]; // tau = tau * theta * cs
568  eta = cs_[0]*cs_[0]*alpha_[0]; // eta = cs^2 * alpha
569 
570  //
571  //--------------------------------------------------------
572  // Update the solution.
573  // Don't update the linear problem object, may incur additional preconditioner application.
574  //--------------------------------------------------------
575  //
576  MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
577  //
578  if (iIter == 1) {
579  //
580  //--------------------------------------------------------
581  // Compute the new rho, beta if we need to.
582  //--------------------------------------------------------
583  //
584  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
585  beta = rho_[0]/rho_old_[0]; // beta = rho / rho_old
586  rho_old_[0] = rho_[0]; // rho_old = rho
587  //
588  //--------------------------------------------------------
589  // Update u, v, and Au if we need to.
590  // Note: We are updating v in two stages to be memory efficient
591  //--------------------------------------------------------
592  //
593  MVT::MvAddMv( STone, *W_, beta, *U_, *U_ ); // u = w + beta*u
594 
595  // First stage of v update.
596  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
597 
598  // Update Au.
599  lp_->apply( *U_, *AU_ ); // Au = A*u
600 
601  // Second stage of v update.
602  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
603  }
604 
605  }
606 
607  // Increment the iteration
608  iter_++;
609 
610  } // end while (sTest_->checkStatus(this) != Passed)
611  }
612 
613 } // namespace Belos
614 //
615 #endif // BELOS_TFQMR_ITER_HPP
616 //
617 // End of file BelosTFQMRIter.hpp
618 
619 
Teuchos::RCP< const MV > V
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Collection of types and exceptions used within the Belos solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
TFQMRIterateFailure(const std::string &what_arg)
void setBlockSize(int blockSize)
Set the blocksize.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
TFQMRIter(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)
Belos::TFQMRIter constructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
int getNumIters() const
Get the current iteration count.
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > U
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Rtilde
void iterate()
This method performs block TFQMR iterations until the status test indicates the need to stop or an er...
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< const MV > W
A linear system to solve, and its associated information.
Structure to contain pointers to TFQMRIter state variables.
Class which describes the linear problem to be solved by the iterative solver.
TFQMRIterInitFailure(const std::string &what_arg)
TFQMRIterInitFailure is thrown when the TFQMRIter object is unable to generate an initial iterate in ...
MultiVecTraits< ScalarType, MV > MVT
TFQMRIterateFailure is thrown when the TFQMRIter object is unable to compute the next iterate in the ...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Teuchos::RCP< const MV > D
Belos header file which uses auto-configuration information to include necessary C++ headers...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated for Belos by doxygen 1.8.14