Belos  Version of the Day
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 "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_ScalarTraits.hpp"
73 #include "Teuchos_ParameterList.hpp"
74 #include "Teuchos_TimeMonitor.hpp"
75 
87 namespace Belos {
88 
93  template <class ScalarType, class MV>
95 
96  typedef Teuchos::ScalarTraits<ScalarType> SCT;
97  typedef typename SCT::magnitudeType MagnitudeType;
98 
100  Teuchos::RCP<const MV> W;
101  Teuchos::RCP<const MV> U;
102  Teuchos::RCP<const MV> AU;
103  Teuchos::RCP<const MV> Rtilde;
104  Teuchos::RCP<const MV> D;
105  Teuchos::RCP<const MV> V;
106  std::vector<ScalarType> alpha, eta, rho;
107  std::vector<MagnitudeType> tau, theta;
108 
109 
110  PseudoBlockTFQMRIterState() : W(Teuchos::null), U(Teuchos::null), AU(Teuchos::null),
111  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
112  {}
113  };
114 
115 
117 
118 
131  PseudoBlockTFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
132  {}};
133 
141  PseudoBlockTFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
142  {}};
143 
145 
146 
147  template <class ScalarType, class MV, class OP>
148  class PseudoBlockTFQMRIter : public Iteration<ScalarType,MV,OP> {
149  public:
150  //
151  // Convenience typedefs
152  //
155  typedef Teuchos::ScalarTraits<ScalarType> SCT;
156  typedef typename SCT::magnitudeType MagnitudeType;
157 
159 
160 
162  PseudoBlockTFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
163  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
164  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
165  Teuchos::ParameterList &params );
166 
168  virtual ~PseudoBlockTFQMRIter() {};
170 
171 
173 
174 
196  void iterate();
197 
220 
224  void initialize()
225  {
227  initializeTFQMR(empty);
228  }
229 
239 
240  // Copy over the vectors.
241  state.W = W_;
242  state.U = U_;
243  state.AU = AU_;
244  state.Rtilde = Rtilde_;
245  state.D = D_;
246  state.V = V_;
247 
248  // Copy over the scalars.
249  state.alpha = alpha_;
250  state.eta = eta_;
251  state.rho = rho_;
252  state.tau = tau_;
253  state.theta = theta_;
254 
255  return state;
256  }
257 
259 
260 
262 
263 
265  int getNumIters() const { return iter_; }
266 
268  void resetNumIters( int iter = 0 ) { iter_ = iter; }
269 
272  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
273 
275 
278  Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
279 
281 
282 
284 
285 
287  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
288 
290  int getBlockSize() const { return 1; }
291 
293  void setBlockSize(int blockSize) {
294  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
295  "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
296  }
297 
299  bool isInitialized() { return initialized_; }
300 
302 
303 
304  private:
305 
306  //
307  // Classes inputed through constructor that define the linear problem to be solved.
308  //
309  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
310  const Teuchos::RCP<OutputManager<ScalarType> > om_;
311  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
312 
313  //
314  // Algorithmic parameters
315  //
316 
317  // numRHS_ is the current number of linear systems being solved.
318  int numRHS_;
319 
320  // Storage for QR factorization of the least squares system.
321  std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
322  std::vector<MagnitudeType> tau_, theta_;
323 
324  //
325  // Current solver state
326  //
327  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
328  // is capable of running; _initialize is controlled by the initialize() member method
329  // For the implications of the state of initialized_, please see documentation for initialize()
330  bool initialized_;
331 
332  // Current subspace dimension, and number of iterations performed.
333  int iter_;
334 
335  //
336  // State Storage
337  //
338  Teuchos::RCP<MV> W_;
339  Teuchos::RCP<MV> U_, AU_;
340  Teuchos::RCP<MV> Rtilde_;
341  Teuchos::RCP<MV> D_;
342  Teuchos::RCP<MV> V_;
343  Teuchos::RCP<MV> solnUpdate_;
344 
345  };
346 
347 
348  //
349  // Implementation
350  //
351 
353  // Constructor.
354  template <class ScalarType, class MV, class OP>
356  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
357  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
358  Teuchos::ParameterList &params
359  ) :
360  lp_(problem),
361  om_(printer),
362  stest_(tester),
363  numRHS_(0),
364  initialized_(false),
365  iter_(0)
366  {
367  }
368 
370  // Compute native residual from TFQMR recurrence.
371  template <class ScalarType, class MV, class OP>
372  Teuchos::RCP<const MV>
373  PseudoBlockTFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
374  {
375  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
376  if (normvec) {
377  // Resize the vector passed in, if it is too small.
378  if ((int) normvec->size() < numRHS_)
379  normvec->resize( numRHS_ );
380 
381  // Compute the native residuals.
382  for (int i=0; i<numRHS_; i++) {
383  (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
384  }
385  }
386 
387  return Teuchos::null;
388  }
389 
391  // Initialize this iteration object
392  template <class ScalarType, class MV, class OP>
394  {
395  // Create convenience variables for zero and one.
396  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
397  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
398  const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
399 
400  // NOTE: In PseudoBlockTFQMRIter Rtilde_, the initial residual, is required!!!
401  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Rtilde == Teuchos::null,std::invalid_argument,
402  "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
403 
404  // Get the number of right-hand sides we're solving for now.
405  int numRHS = MVT::GetNumberVecs(*newstate.Rtilde);
406  numRHS_ = numRHS;
407 
408  // Initialize the state storage
409  // If the subspace has not be initialized before or we are reusing this solver object, generate it using Rtilde.
410  if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
411  {
412  // Create and/or initialize D_.
413  if ( Teuchos::is_null(D_) )
414  D_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
415  MVT::MvInit( *D_, STzero );
416 
417  // Create and/or initialize solnUpdate_;
418  if ( Teuchos::is_null(solnUpdate_) )
419  solnUpdate_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
420  MVT::MvInit( *solnUpdate_, STzero );
421 
422  // Create Rtilde_.
423  if (newstate.Rtilde != Rtilde_)
424  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
425  W_ = MVT::CloneCopy( *Rtilde_ );
426  U_ = MVT::CloneCopy( *Rtilde_ );
427  V_ = MVT::Clone( *Rtilde_, numRHS_ );
428 
429  // Multiply the current residual by Op and store in V_
430  // V_ = Op * R_
431  lp_->apply( *U_, *V_ );
432  AU_ = MVT::CloneCopy( *V_ );
433 
434  // Resize work vectors.
435  alpha_.resize( numRHS_, STone );
436  eta_.resize( numRHS_, STzero );
437  rho_.resize( numRHS_ );
438  rho_old_.resize( numRHS_ );
439  tau_.resize( numRHS_ );
440  theta_.resize( numRHS_, MTzero );
441 
442  MVT::MvNorm( *Rtilde_, tau_ ); // tau = ||r_0||
443  MVT::MvDot( *Rtilde_, *Rtilde_, rho_ ); // rho = (r_tilde, r0)
444  }
445  else
446  {
447  // If the subspace has changed sizes, clone it from the incoming state.
448  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
449  W_ = MVT::CloneCopy( *newstate.W );
450  U_ = MVT::CloneCopy( *newstate.U );
451  AU_ = MVT::CloneCopy( *newstate.AU );
452  D_ = MVT::CloneCopy( *newstate.D );
453  V_ = MVT::CloneCopy( *newstate.V );
454 
455  // The solution update is just set to zero, since the current update has already
456  // been added to the solution during deflation.
457  solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
458  MVT::MvInit( *solnUpdate_, STzero );
459 
460  // Copy work vectors.
461  alpha_ = newstate.alpha;
462  eta_ = newstate.eta;
463  rho_ = newstate.rho;
464  tau_ = newstate.tau;
465  theta_ = newstate.theta;
466  }
467 
468  // The solver is initialized
469  initialized_ = true;
470  }
471 
472 
474  // Iterate until the status test informs us we should stop.
475  template <class ScalarType, class MV, class OP>
477  {
478  //
479  // Allocate/initialize data structures
480  //
481  if (initialized_ == false) {
482  initialize();
483  }
484 
485  // Create convenience variables for zero and one.
486  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
487  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
488  const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
489  std::vector< ScalarType > beta(numRHS_,STzero);
490  std::vector<int> index(1);
491  //
492  // Start executable statements.
493  //
494 
496  // Iterate until the status test tells us to stop.
497  //
498  while (stest_->checkStatus(this) != Passed) {
499 
500  for (int iIter=0; iIter<2; iIter++)
501  {
502  //
503  //--------------------------------------------------------
504  // Compute the new alpha if we need to
505  //--------------------------------------------------------
506  //
507  if (iIter == 0) {
508  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
509  for (int i=0; i<numRHS_; i++) {
510  rho_old_[i] = rho_[i]; // rho_old = rho
511  alpha_[i] = rho_old_[i]/alpha_[i];
512  }
513  }
514  //
515  //--------------------------------------------------------
516  // Loop over all RHS and compute updates.
517  //--------------------------------------------------------
518  //
519  for (int i=0; i<numRHS_; ++i) {
520  index[0] = i;
521 
522  //
523  //--------------------------------------------------------
524  // Update w.
525  // w = w - alpha*Au
526  //--------------------------------------------------------
527  //
528  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
529  Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
530  MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
531  //
532  //--------------------------------------------------------
533  // Update d.
534  // d = u + (theta^2/alpha)eta*d
535  //--------------------------------------------------------
536  //
537  Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
538  Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
539  MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
540  //
541  //--------------------------------------------------------
542  // Update u if we need to.
543  // u = u - alpha*v
544  //
545  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
546  //--------------------------------------------------------
547  //
548  if (iIter == 0) {
549  // Compute new U.
550  Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
551  Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
552  MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
553  }
554  }
555  //
556  //--------------------------------------------------------
557  // Update Au for the next iteration.
558  //--------------------------------------------------------
559  //
560  if (iIter == 0) {
561  lp_->apply( *U_, *AU_ );
562  }
563  //
564  //--------------------------------------------------------
565  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
566  //--------------------------------------------------------
567  //
568  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
569 
570  for (int i=0; i<numRHS_; ++i) {
571  theta_[i] /= tau_[i];
572  // cs = 1.0 / sqrt(1.0 + theta^2)
573  MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
574  tau_[i] *= theta_[i]*cs; // tau = tau * theta * cs
575  eta_[i] = cs*cs*alpha_[i]; // eta = cs^2 * alpha
576  }
577  //
578  //--------------------------------------------------------
579  // Accumulate the update for the solution x := x + eta*D_
580  //--------------------------------------------------------
581  //
582  for (int i=0; i<numRHS_; ++i) {
583  index[0]=i;
584  Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
585  Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
586  MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
587  }
588  //
589  if (iIter == 1) {
590  //
591  //--------------------------------------------------------
592  // Compute the new rho, beta if we need to.
593  //--------------------------------------------------------
594  //
595  MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
596 
597  for (int i=0; i<numRHS_; ++i) {
598  beta[i] = rho_[i]/rho_old_[i]; // beta = rho / rho_old
599 
600  //
601  //--------------------------------------------------------
602  // Update u, v, and Au if we need to.
603  // Note: We are updating v in two stages to be memory efficient
604  //--------------------------------------------------------
605  //
606  index[0]=i;
607  Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
608  Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
609  MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i ); // u = w + beta*u
610 
611  // First stage of v update.
612  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
613  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
614  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
615  }
616 
617  // Update Au.
618  lp_->apply( *U_, *AU_ ); // Au = A*u
619 
620  // Second stage of v update.
621  for (int i=0; i<numRHS_; ++i) {
622  index[0]=i;
623  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
624  Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
625  MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
626  }
627  }
628  }
629 
630  // Increment the iteration
631  iter_++;
632 
633  } // end while (sTest_->checkStatus(this) != Passed)
634  }
635 
636 } // namespace Belos
637 //
638 #endif // BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
639 //
640 // End of file BelosPseudoBlockTFQMRIter.hpp
641 
642 
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.
int getNumIters() const
Get the current iteration count.
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.
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
PseudoBlockTFQMRIterInitFailure(const std::string &what_arg)
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::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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
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 block TFQMR iterations until the status test indicates the need to stop or an er...
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.

Generated for Belos by doxygen 1.8.14