Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "BelosOutputManager.hpp"
66 #include "BelosStatusTest.hpp"
67 #include "BelosOperatorTraits.hpp"
68 #include "BelosMultiVecTraits.hpp"
69 
70 #include "Teuchos_BLAS.hpp"
73 #include "Teuchos_ScalarTraits.hpp"
75 #include "Teuchos_TimeMonitor.hpp"
76 
88 namespace Belos {
89 
94  template <class ScalarType, class MV>
95  struct TFQMRIterState {
96 
104 
105  TFQMRIterState() : R(Teuchos::null), W(Teuchos::null), U(Teuchos::null),
106  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
107  {}
108  };
109 
110 
112 
113 
125  class TFQMRIterInitFailure : public BelosError {public:
126  TFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
127  {}};
128 
135  class TFQMRIterateFailure : public BelosError {public:
136  TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
137  {}};
138 
140 
141 
142  template <class ScalarType, class MV, class OP>
143  class TFQMRIter : public Iteration<ScalarType,MV,OP> {
144  public:
145  //
146  // Convenience typedefs
147  //
152 
154 
155 
158  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
160  Teuchos::ParameterList &params );
161 
163  virtual ~TFQMRIter() {};
165 
166 
168 
169 
180  void iterate();
181 
203  void initializeTFQMR(const TFQMRIterState<ScalarType,MV> & newstate);
204 
208  void initialize()
209  {
211  initializeTFQMR(empty);
212  }
213 
223  state.R = R_;
224  state.W = W_;
225  state.U = U_;
226  state.Rtilde = Rtilde_;
227  state.D = D_;
228  state.V = V_;
229  state.solnUpdate = solnUpdate_;
230  return state;
231  }
232 
234 
235 
237 
238 
240  int getNumIters() const { return iter_; }
241 
243  void resetNumIters( int iter = 0 ) { iter_ = iter; }
244 
247  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
248 
250 
253  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
254 
256 
257 
259 
260 
262  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
263 
265  int getBlockSize() const { return 1; }
266 
268  void setBlockSize(int blockSize) {
269  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
270  "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
271  }
272 
274  bool isInitialized() { return initialized_; }
275 
277 
278 
279  private:
280 
281  //
282  // Internal methods
283  //
285  void setStateSize();
286 
287  //
288  // Classes inputed through constructor that define the linear problem to be solved.
289  //
293 
294  //
295  // Algorithmic parameters
296  //
297 
298  // Storage for QR factorization of the least squares system.
299  // Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
300  std::vector<ScalarType> alpha_, rho_, rho_old_;
301  std::vector<MagnitudeType> tau_, cs_, theta_;
302 
303  //
304  // Current solver state
305  //
306  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
307  // is capable of running; _initialize is controlled by the initialize() member method
308  // For the implications of the state of initialized_, please see documentation for initialize()
309  bool initialized_;
310 
311  // stateStorageInitialized_ specifies that the state storage has be initialized to the current
312  // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
313  // generated without the right-hand side or solution vectors.
314  bool stateStorageInitialized_;
315 
316  // Current subspace dimension, and number of iterations performed.
317  int iter_;
318 
319  //
320  // State Storage
321  //
322  Teuchos::RCP<MV> R_;
323  Teuchos::RCP<MV> W_;
324  Teuchos::RCP<MV> U_, AU_;
325  Teuchos::RCP<MV> Rtilde_;
326  Teuchos::RCP<MV> D_;
327  Teuchos::RCP<MV> V_;
328  Teuchos::RCP<MV> solnUpdate_;
329  };
330 
331 
332  //
333  // Implementation
334  //
335 
337  // Constructor.
338  template <class ScalarType, class MV, class OP>
340  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
342  Teuchos::ParameterList &/* params */
343  ) :
344  lp_(problem),
345  om_(printer),
346  stest_(tester),
347  alpha_(1),
348  rho_(1),
349  rho_old_(1),
350  tau_(1),
351  cs_(1),
352  theta_(1),
353  initialized_(false),
354  stateStorageInitialized_(false),
355  iter_(0)
356  {
357  }
358 
360  // Compute native residual from TFQMR recurrence.
361  template <class ScalarType, class MV, class OP>
363  TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
364  {
366  if (normvec)
367  (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
368 
369  return Teuchos::null;
370  }
371 
372 
374  // Setup the state storage.
375  template <class ScalarType, class MV, class OP>
377  {
378  if (!stateStorageInitialized_) {
379 
380  // Check if there is any multivector to clone from.
381  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
382  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
383  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
384  stateStorageInitialized_ = false;
385  return;
386  }
387  else {
388 
389  // Initialize the state storage
390  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
391  if (R_ == Teuchos::null) {
392  // Get the multivector that is not null.
393  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
394  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
395  "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
396  R_ = MVT::Clone( *tmp, 1 );
397  D_ = MVT::Clone( *tmp, 1 );
398  V_ = MVT::Clone( *tmp, 1 );
399  solnUpdate_ = MVT::Clone( *tmp, 1 );
400  }
401 
402  // State storage has now been initialized.
403  stateStorageInitialized_ = true;
404  }
405  }
406  }
407 
409  // Initialize this iteration object
410  template <class ScalarType, class MV, class OP>
412  {
413  // Initialize the state storage if it isn't already.
414  if (!stateStorageInitialized_)
415  setStateSize();
416 
417  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
418  "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
419 
420  // NOTE: In TFQMRIter R_, the initial residual, is required!!!
421  //
422  std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
423 
424  // Create convenience variables for zero and one.
426 
427  if (newstate.R != Teuchos::null) {
428 
429  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
430  std::invalid_argument, errstr );
431  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
432  std::invalid_argument, errstr );
433 
434  // Copy basis vectors from newstate into V
435  if (newstate.R != R_) {
436  // copy over the initial residual (unpreconditioned).
437  MVT::Assign( *newstate.R, *R_ );
438  }
439 
440  // Compute initial vectors
441  // Initially, they are set to the preconditioned residuals
442  //
443  W_ = MVT::CloneCopy( *R_ );
444  U_ = MVT::CloneCopy( *R_ );
445  Rtilde_ = MVT::CloneCopy( *R_ );
446  MVT::MvInit( *D_ );
447  MVT::MvInit( *solnUpdate_ );
448  // Multiply the current residual by Op and store in V_
449  // V_ = Op * R_
450  //
451  lp_->apply( *U_, *V_ );
452  AU_ = MVT::CloneCopy( *V_ );
453  //
454  // Compute initial scalars: theta, eta, tau, rho_old
455  //
456  theta_[0] = MTzero;
457  MVT::MvNorm( *R_, tau_ ); // tau = ||r_0||
458  MVT::MvDot( *R_, *Rtilde_, rho_old_ ); // rho = (r_tilde, r0)
459  }
460  else {
461 
462  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
463  "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
464  }
465 
466  // The solver is initialized
467  initialized_ = true;
468  }
469 
470 
472  // Iterate until the status test informs us we should stop.
473  template <class ScalarType, class MV, class OP>
475  {
476  //
477  // Allocate/initialize data structures
478  //
479  if (initialized_ == false) {
480  initialize();
481  }
482 
483  // Create convenience variables for zero and one.
484  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
487  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
488  ScalarType eta = STzero, beta = STzero;
489  //
490  // Start executable statements.
491  //
492  // Get the current solution vector.
493  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
494 
495  // Check that the current solution vector only has one column.
496  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
497  "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
498 
499 
501  // Iterate until the status test tells us to stop.
502  //
503  while (stest_->checkStatus(this) != Passed) {
504 
505  for (int iIter=0; iIter<2; iIter++)
506  {
507  //
508  //--------------------------------------------------------
509  // Compute the new alpha if we need to
510  //--------------------------------------------------------
511  //
512  if (iIter == 0) {
513  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
514  alpha_[0] = rho_old_[0]/alpha_[0];
515  }
516  //
517  //--------------------------------------------------------
518  // Update w.
519  // w = w - alpha*Au
520  //--------------------------------------------------------
521  //
522  MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
523  //
524  //--------------------------------------------------------
525  // Update d.
526  // d = u + (theta^2/alpha)eta*d
527  //--------------------------------------------------------
528  //
529  MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
530  //
531  //--------------------------------------------------------
532  // Update u if we need to.
533  // u = u - alpha*v
534  //
535  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
536  //--------------------------------------------------------
537  //
538  if (iIter == 0) {
539  // Compute new U.
540  MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
541 
542  // Update Au for the next iteration.
543  lp_->apply( *U_, *AU_ );
544  }
545  //
546  //--------------------------------------------------------
547  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
548  //--------------------------------------------------------
549  //
550  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
551  theta_[0] /= tau_[0];
552  // cs = 1.0 / sqrt(1.0 + theta^2)
553  cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
554  tau_[0] *= theta_[0]*cs_[0]; // tau = tau * theta * cs
555  eta = cs_[0]*cs_[0]*alpha_[0]; // eta = cs^2 * alpha
556  //
557  //--------------------------------------------------------
558  // Update the solution.
559  // Don't update the linear problem object, may incur additional preconditioner application.
560  //--------------------------------------------------------
561  //
562  MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
563  //
564  //--------------------------------------------------------
565  // Check for breakdown before continuing.
566  //--------------------------------------------------------
567  if ( tau_[0] == MTzero ) {
568  break;
569  }
570  //
571  if (iIter == 1) {
572  //
573  //--------------------------------------------------------
574  // Compute the new rho, beta if we need to.
575  //--------------------------------------------------------
576  //
577  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
578  beta = rho_[0]/rho_old_[0]; // beta = rho / rho_old
579  rho_old_[0] = rho_[0]; // rho_old = rho
580  //
581  //--------------------------------------------------------
582  // Update u, v, and Au if we need to.
583  // Note: We are updating v in two stages to be memory efficient
584  //--------------------------------------------------------
585  //
586  MVT::MvAddMv( STone, *W_, beta, *U_, *U_ ); // u = w + beta*u
587 
588  // First stage of v update.
589  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
590 
591  // Update Au.
592  lp_->apply( *U_, *AU_ ); // Au = A*u
593 
594  // Second stage of v update.
595  MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
596  }
597 
598  }
599 
600  // Increment the iteration
601  iter_++;
602 
603  } // end while (sTest_->checkStatus(this) != Passed)
604  }
605 
606 } // namespace Belos
607 //
608 #endif // BELOS_TFQMR_ITER_HPP
609 //
610 // End of file BelosTFQMRIter.hpp
611 
612 
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.
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.
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.
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.
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
void iterate()
This method performs TFQMR iterations until the status test indicates the need to stop or an error oc...
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)
int getNumIters() const
Get the current iteration count.
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.

Generated on Fri Aug 14 2020 10:48:34 for Belos by doxygen 1.8.5