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

Generated for Belos by doxygen 1.8.14