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"
72 #include "Teuchos_SerialDenseMatrix.hpp"
73 #include "Teuchos_SerialDenseVector.hpp"
74 #include "Teuchos_ScalarTraits.hpp"
75 #include "Teuchos_ParameterList.hpp"
76 #include "Teuchos_TimeMonitor.hpp"
77 
89 namespace Belos {
90 
95  template <class ScalarType, class MV>
96  struct TFQMRIterState {
97 
99  Teuchos::RCP<const MV> R;
100  Teuchos::RCP<const MV> W;
101  Teuchos::RCP<const MV> U;
102  Teuchos::RCP<const MV> Rtilde;
103  Teuchos::RCP<const MV> D;
104  Teuchos::RCP<const MV> V;
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  //
151  typedef Teuchos::ScalarTraits<ScalarType> SCT;
152  typedef typename SCT::magnitudeType MagnitudeType;
153 
155 
156 
158  TFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
159  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
160  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
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  //
302  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
303  const Teuchos::RCP<OutputManager<ScalarType> > om_;
304  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
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<MagnitudeType> tau_, cs_, theta_;
313 
314  //
315  // Current solver state
316  //
317  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
318  // is capable of running; _initialize is controlled by the initialize() member method
319  // For the implications of the state of initialized_, please see documentation for initialize()
320  bool initialized_;
321 
322  // stateStorageInitialized_ specifies that the state storage has be initialized to the current
323  // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
324  // generated without the right-hand side or solution vectors.
325  bool stateStorageInitialized_;
326 
327  // Current subspace dimension, and number of iterations performed.
328  int iter_;
329 
330  //
331  // State Storage
332  //
333  Teuchos::RCP<MV> R_;
334  Teuchos::RCP<MV> W_;
335  Teuchos::RCP<MV> U_, AU_;
336  Teuchos::RCP<MV> Rtilde_;
337  Teuchos::RCP<MV> D_;
338  Teuchos::RCP<MV> V_;
339  Teuchos::RCP<MV> solnUpdate_;
340  };
341 
342 
343  //
344  // Implementation
345  //
346 
348  // Constructor.
349  template <class ScalarType, class MV, class OP>
351  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
352  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
353  Teuchos::ParameterList &params
354  ) :
355  lp_(problem),
356  om_(printer),
357  stest_(tester),
358  alpha_(1,1),
359  rho_(1,1),
360  rho_old_(1,1),
361  tau_(1),
362  cs_(1),
363  theta_(1),
364  initialized_(false),
365  stateStorageInitialized_(false),
366  iter_(0)
367  {
368  }
369 
371  // Compute native residual from TFQMR recurrence.
372  template <class ScalarType, class MV, class OP>
373  Teuchos::RCP<const MV>
374  TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
375  {
376  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
377  if (normvec)
378  (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( iter_ + one )*tau_[0];
379 
380  return Teuchos::null;
381  }
382 
383 
385  // Setup the state storage.
386  template <class ScalarType, class MV, class OP>
388  {
389  if (!stateStorageInitialized_) {
390 
391  // Check if there is any multivector to clone from.
392  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
393  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
394  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
395  stateStorageInitialized_ = false;
396  return;
397  }
398  else {
399 
400  // Initialize the state storage
401  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
402  if (R_ == Teuchos::null) {
403  // Get the multivector that is not null.
404  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
405  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
406  "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
407  R_ = MVT::Clone( *tmp, 1 );
408  AU_ = MVT::Clone( *tmp, 1 );
409  D_ = MVT::Clone( *tmp, 1 );
410  solnUpdate_ = MVT::Clone( *tmp, 1 );
411  V_ = 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();
439  const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::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::MvTransMv( one, *Rtilde_, *R_, 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();
499  const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::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  //
519  //--------------------------------------------------------
520  // Compute the new alpha if we need to
521  //--------------------------------------------------------
522  //
523  if (iter_%2 == 0) {
524  MVT::MvTransMv( STone, *Rtilde_, *V_, alpha_ ); // alpha = rho / (r_tilde, v)
525  alpha_(0,0) = rho_old_(0,0)/alpha_(0,0);
526  }
527  //
528  //--------------------------------------------------------
529  // Update w.
530  // w = w - alpha*Au
531  //--------------------------------------------------------
532  //
533  MVT::MvAddMv( STone, *W_, -alpha_(0,0), *AU_, *W_ );
534  //
535  //--------------------------------------------------------
536  // Update d.
537  // d = u + (theta^2/alpha)eta*d
538  //--------------------------------------------------------
539  //
540  MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_(0,0))*eta, *D_, *D_ );
541  //
542  //--------------------------------------------------------
543  // Update u if we need to.
544  // u = u - alpha*v
545  //
546  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
547  //--------------------------------------------------------
548  //
549  if (iter_%2 == 0) {
550  // Compute new U.
551  MVT::MvAddMv( STone, *U_, -alpha_(0,0), *V_, *U_ );
552 
553  // Update Au for the next iteration.
554  lp_->apply( *U_, *AU_ );
555  }
556  //
557  //--------------------------------------------------------
558  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
559  //--------------------------------------------------------
560  //
561  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
562  theta_[0] /= tau_[0];
563  // cs = 1.0 / sqrt(1.0 + theta^2)
564  cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
565  tau_[0] *= theta_[0]*cs_[0]; // tau = tau * theta * cs
566  eta = cs_[0]*cs_[0]*alpha_(0,0); // eta = cs^2 * alpha
567 
568  //
569  //--------------------------------------------------------
570  // Update the solution.
571  // Don't update the linear problem object, may incur additional preconditioner application.
572  //--------------------------------------------------------
573  //
574  MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
575  //
576  if (iter_%2) {
577  //
578  //--------------------------------------------------------
579  // Compute the new rho, beta if we need to.
580  //--------------------------------------------------------
581  //
582  MVT::MvTransMv( STone, *Rtilde_, *W_, rho_ ); // rho = (r_tilde, w)
583  beta = rho_(0,0)/rho_old_(0,0); // beta = rho / rho_old
584  rho_old_(0,0) = rho_(0,0); // rho_old = rho
585  //
586  //--------------------------------------------------------
587  // Update u, v, and Au if we need to.
588  // Note: We are updating v in two stages to be memory efficient
589  //--------------------------------------------------------
590  //
591  MVT::MvAddMv( STone, *W_, beta, *U_, *U_ ); // u = w + beta*u
592 
593  // First stage of v update.
594  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
595 
596  // Update Au.
597  lp_->apply( *U_, *AU_ ); // Au = A*u
598 
599  // Second stage of v update.
600  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
601  }
602 
603  // Increment the iteration
604  iter_++;
605 
606  } // end while (sTest_->checkStatus(this) != Passed)
607  }
608 
609 } // namespace Belos
610 //
611 #endif // BELOS_TFQMR_ITER_HPP
612 //
613 // End of file BelosTFQMRIter.hpp
614 
615 
Teuchos::RCP< const MV > V
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
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.
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.
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< const MV > U
void resetNumIters(int iter=0)
Reset the iteration count.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
Pure virtual base class which describes the basic interface to the linear solver iteration.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > Rtilde
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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)
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
int getNumIters() const
Get the current iteration count.
TFQMRIterInitFailure is thrown when the TFQMRIter object is unable to generate an initial iterate in ...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11