Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockTFQMRIter.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_PSEUDO_BLOCK_TFQMR_ITER_HPP
51 #define BELOS_PSEUDO_BLOCK_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"
71 #include "Teuchos_ScalarTraits.hpp"
73 #include "Teuchos_TimeMonitor.hpp"
74 
86 namespace Belos {
87 
92  template <class ScalarType, class MV>
94 
97 
105  std::vector<ScalarType> alpha, eta, rho;
106  std::vector<MagnitudeType> tau, theta;
107 
108 
109  PseudoBlockTFQMRIterState() : W(Teuchos::null), U(Teuchos::null), AU(Teuchos::null),
110  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
111  {}
112  };
113 
114 
116 
117 
130  PseudoBlockTFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
131  {}};
132 
140  PseudoBlockTFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
141  {}};
142 
144 
145 
146  template <class ScalarType, class MV, class OP>
147  class PseudoBlockTFQMRIter : public Iteration<ScalarType,MV,OP> {
148  public:
149  //
150  // Convenience typedefs
151  //
156 
158 
159 
162  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
164  Teuchos::ParameterList &params );
165 
167  virtual ~PseudoBlockTFQMRIter() {};
169 
170 
172 
173 
184  void iterate();
185 
208 
212  void initialize()
213  {
215  initializeTFQMR(empty);
216  }
217 
227 
228  // Copy over the vectors.
229  state.W = W_;
230  state.U = U_;
231  state.AU = AU_;
232  state.Rtilde = Rtilde_;
233  state.D = D_;
234  state.V = V_;
235 
236  // Copy over the scalars.
237  state.alpha = alpha_;
238  state.eta = eta_;
239  state.rho = rho_;
240  state.tau = tau_;
241  state.theta = theta_;
242 
243  return state;
244  }
245 
247 
248 
250 
251 
253  int getNumIters() const { return iter_; }
254 
256  void resetNumIters( int iter = 0 ) { iter_ = iter; }
257 
260  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
261 
263 
266  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
267 
269 
270 
272 
273 
275  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
276 
278  int getBlockSize() const { return 1; }
279 
281  void setBlockSize(int blockSize) {
282  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
283  "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
284  }
285 
287  bool isInitialized() { return initialized_; }
288 
290 
291 
292  private:
293 
294  //
295  // Classes inputed through constructor that define the linear problem to be solved.
296  //
300 
301  //
302  // Algorithmic parameters
303  //
304 
305  // numRHS_ is the current number of linear systems being solved.
306  int numRHS_;
307 
308  // Storage for QR factorization of the least squares system.
309  std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
310  std::vector<MagnitudeType> tau_, theta_;
311 
312  //
313  // Current solver state
314  //
315  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
316  // is capable of running; _initialize is controlled by the initialize() member method
317  // For the implications of the state of initialized_, please see documentation for initialize()
318  bool initialized_;
319 
320  // Current subspace dimension, and number of iterations performed.
321  int iter_;
322 
323  //
324  // State Storage
325  //
326  Teuchos::RCP<MV> W_;
327  Teuchos::RCP<MV> U_, AU_;
328  Teuchos::RCP<MV> Rtilde_;
329  Teuchos::RCP<MV> D_;
330  Teuchos::RCP<MV> V_;
331  Teuchos::RCP<MV> solnUpdate_;
332 
333  };
334 
335 
336  //
337  // Implementation
338  //
339 
341  // Constructor.
342  template <class ScalarType, class MV, class OP>
344  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
346  Teuchos::ParameterList &/* params */
347  ) :
348  lp_(problem),
349  om_(printer),
350  stest_(tester),
351  numRHS_(0),
352  initialized_(false),
353  iter_(0)
354  {
355  }
356 
358  // Compute native residual from TFQMR recurrence.
359  template <class ScalarType, class MV, class OP>
361  PseudoBlockTFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
362  {
364  if (normvec) {
365  // Resize the vector passed in, if it is too small.
366  if ((int) normvec->size() < numRHS_)
367  normvec->resize( numRHS_ );
368 
369  // Compute the native residuals.
370  for (int i=0; i<numRHS_; i++) {
371  (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
372  }
373  }
374 
375  return Teuchos::null;
376  }
377 
379  // Initialize this iteration object
380  template <class ScalarType, class MV, class OP>
382  {
383  // Create convenience variables for zero and one.
384  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
385  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
387 
388  // NOTE: In PseudoBlockTFQMRIter Rtilde_, the initial residual, is required!!!
389  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Rtilde == Teuchos::null,std::invalid_argument,
390  "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
391 
392  // Get the number of right-hand sides we're solving for now.
393  int numRHS = MVT::GetNumberVecs(*newstate.Rtilde);
394  numRHS_ = numRHS;
395 
396  // Initialize the state storage
397  // If the subspace has not be initialized before or we are reusing this solver object, generate it using Rtilde.
398  if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
399  {
400  // Create and/or initialize D_.
401  if ( Teuchos::is_null(D_) )
402  D_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
403  MVT::MvInit( *D_, STzero );
404 
405  // Create and/or initialize solnUpdate_;
406  if ( Teuchos::is_null(solnUpdate_) )
407  solnUpdate_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
408  MVT::MvInit( *solnUpdate_, STzero );
409 
410  // Create Rtilde_.
411  if (newstate.Rtilde != Rtilde_)
412  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
413  W_ = MVT::CloneCopy( *Rtilde_ );
414  U_ = MVT::CloneCopy( *Rtilde_ );
415  V_ = MVT::Clone( *Rtilde_, numRHS_ );
416 
417  // Multiply the current residual by Op and store in V_
418  // V_ = Op * R_
419  lp_->apply( *U_, *V_ );
420  AU_ = MVT::CloneCopy( *V_ );
421 
422  // Resize work vectors.
423  alpha_.resize( numRHS_, STone );
424  eta_.resize( numRHS_, STzero );
425  rho_.resize( numRHS_ );
426  rho_old_.resize( numRHS_ );
427  tau_.resize( numRHS_ );
428  theta_.resize( numRHS_, MTzero );
429 
430  MVT::MvNorm( *Rtilde_, tau_ ); // tau = ||r_0||
431  MVT::MvDot( *Rtilde_, *Rtilde_, rho_ ); // rho = (r_tilde, r0)
432  }
433  else
434  {
435  // If the subspace has changed sizes, clone it from the incoming state.
436  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
437  W_ = MVT::CloneCopy( *newstate.W );
438  U_ = MVT::CloneCopy( *newstate.U );
439  AU_ = MVT::CloneCopy( *newstate.AU );
440  D_ = MVT::CloneCopy( *newstate.D );
441  V_ = MVT::CloneCopy( *newstate.V );
442 
443  // The solution update is just set to zero, since the current update has already
444  // been added to the solution during deflation.
445  solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
446  MVT::MvInit( *solnUpdate_, STzero );
447 
448  // Copy work vectors.
449  alpha_ = newstate.alpha;
450  eta_ = newstate.eta;
451  rho_ = newstate.rho;
452  tau_ = newstate.tau;
453  theta_ = newstate.theta;
454  }
455 
456  // The solver is initialized
457  initialized_ = true;
458  }
459 
460 
462  // Iterate until the status test informs us we should stop.
463  template <class ScalarType, class MV, class OP>
465  {
466  //
467  // Allocate/initialize data structures
468  //
469  if (initialized_ == false) {
470  initialize();
471  }
472 
473  // Create convenience variables for zero and one.
474  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
475  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
478  std::vector< ScalarType > beta(numRHS_,STzero);
479  std::vector<int> index(1);
480  //
481  // Start executable statements.
482  //
483 
485  // Iterate until the status test tells us to stop.
486  //
487  while (stest_->checkStatus(this) != Passed) {
488 
489  for (int iIter=0; iIter<2; iIter++)
490  {
491  //
492  //--------------------------------------------------------
493  // Compute the new alpha if we need to
494  //--------------------------------------------------------
495  //
496  if (iIter == 0) {
497  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
498  for (int i=0; i<numRHS_; i++) {
499  rho_old_[i] = rho_[i]; // rho_old = rho
500  alpha_[i] = rho_old_[i]/alpha_[i];
501  }
502  }
503  //
504  //--------------------------------------------------------
505  // Loop over all RHS and compute updates.
506  //--------------------------------------------------------
507  //
508  for (int i=0; i<numRHS_; ++i) {
509  index[0] = i;
510 
511  //
512  //--------------------------------------------------------
513  // Update w.
514  // w = w - alpha*Au
515  //--------------------------------------------------------
516  //
517  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
518  Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
519  MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
520  //
521  //--------------------------------------------------------
522  // Update d.
523  // d = u + (theta^2/alpha)eta*d
524  //--------------------------------------------------------
525  //
526  Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
527  Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
528  MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
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  Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
540  Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
541  MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
542  }
543  }
544  //
545  //--------------------------------------------------------
546  // Update Au for the next iteration.
547  //--------------------------------------------------------
548  //
549  if (iIter == 0) {
550  lp_->apply( *U_, *AU_ );
551  }
552  //
553  //--------------------------------------------------------
554  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
555  //--------------------------------------------------------
556  //
557  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
558 
559  bool breakdown=false;
560  for (int i=0; i<numRHS_; ++i) {
561  theta_[i] /= tau_[i];
562  // cs = 1.0 / sqrt(1.0 + theta^2)
563  MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
564  tau_[i] *= theta_[i]*cs; // tau = tau * theta * cs
565  if ( tau_[i] == MTzero ) {
566  breakdown = true;
567  }
568  eta_[i] = cs*cs*alpha_[i]; // eta = cs^2 * alpha
569  }
570  //
571  //--------------------------------------------------------
572  // Accumulate the update for the solution x := x + eta*D_
573  //--------------------------------------------------------
574  //
575  for (int i=0; i<numRHS_; ++i) {
576  index[0]=i;
577  Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
578  Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
579  MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
580  }
581  //
582  //--------------------------------------------------------
583  // Breakdown was detected above, return to status test to
584  // remove converged solutions.
585  //--------------------------------------------------------
586  if ( breakdown ) {
587  break;
588  }
589  //
590  if (iIter == 1) {
591  //
592  //--------------------------------------------------------
593  // Compute the new rho, beta if we need to.
594  //--------------------------------------------------------
595  //
596  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
597 
598  for (int i=0; i<numRHS_; ++i) {
599  beta[i] = rho_[i]/rho_old_[i]; // beta = rho / rho_old
600 
601  //
602  //--------------------------------------------------------
603  // Update u, v, and Au if we need to.
604  // Note: We are updating v in two stages to be memory efficient
605  //--------------------------------------------------------
606  //
607  index[0]=i;
608  Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
609  Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
610  MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i ); // u = w + beta*u
611 
612  // First stage of v update.
613  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
614  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
615  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
616  }
617 
618  // Update Au.
619  lp_->apply( *U_, *AU_ ); // Au = A*u
620 
621  // Second stage of v update.
622  for (int i=0; i<numRHS_; ++i) {
623  index[0]=i;
624  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
625  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
626  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
627  }
628  }
629  }
630 
631  // Increment the iteration
632  iter_++;
633 
634  } // end while (sTest_->checkStatus(this) != Passed)
635  }
636 
637 } // namespace Belos
638 //
639 #endif // BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
640 //
641 // End of file BelosPseudoBlockTFQMRIter.hpp
642 
643 
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...
OperatorTraits< ScalarType, MV, OP > OPT
Class which manages the output and verbosity of the Belos solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
PseudoBlockTFQMRIterateFailure is thrown when the PseudoBlockTFQMRIter object is unable to compute th...
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
bool is_null(const std::shared_ptr< T > &p)
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getNumIters() const
Get the current iteration count.
Pure virtual base class which describes the basic interface to the linear solver iteration.
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.
Class which defines basic traits for the operator type.
PseudoBlockTFQMRIterInitFailure(const std::string &what_arg)
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Traits class which defines basic operations on multivectors.
PseudoBlockTFQMRIter(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::PseudoBlockTFQMRIter constructor.
PseudoBlockTFQMRIterInitFailure is thrown when the PseudoBlockTFQMRIter object is unable to generate ...
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::ScalarTraits< ScalarType > SCT
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockTFQMRIterateFailure(const std::string &what_arg)
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
Teuchos::RCP< const MV > W
The current residual basis.
Class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType > SCT
void iterate()
This method performs pseudo-block TFQMR iterations until the status test indicates the need to stop o...
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
void resetNumIters(int iter=0)
Reset the iteration count.
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool isInitialized()
States whether the solver has been initialized or not.
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.

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