Belos Package Browser (Single Doxygen Collection)  Development
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_ScalarTraits.hpp"
71 #include "Teuchos_ParameterList.hpp"
72 #include "Teuchos_TimeMonitor.hpp"
73 
85 namespace Belos {
86 
91  template <class ScalarType, class MV>
93 
94  typedef Teuchos::ScalarTraits<ScalarType> SCT;
95  typedef typename SCT::magnitudeType MagnitudeType;
96 
98  Teuchos::RCP<const MV> W;
99  Teuchos::RCP<const MV> U;
100  Teuchos::RCP<const MV> AU;
101  Teuchos::RCP<const MV> Rtilde;
102  Teuchos::RCP<const MV> D;
103  Teuchos::RCP<const MV> V;
104  std::vector<ScalarType> alpha, eta, rho;
105  std::vector<MagnitudeType> tau, theta;
106 
107 
108  PseudoBlockTFQMRIterState() : W(Teuchos::null), U(Teuchos::null), AU(Teuchos::null),
109  Rtilde(Teuchos::null), D(Teuchos::null), V(Teuchos::null)
110  {}
111  };
112 
113 
115 
116 
129  PseudoBlockTFQMRIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
130  {}};
131 
139  PseudoBlockTFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
140  {}};
141 
143 
144 
145  template <class ScalarType, class MV, class OP>
146  class PseudoBlockTFQMRIter : public Iteration<ScalarType,MV,OP> {
147  public:
148  //
149  // Convenience typedefs
150  //
153  typedef Teuchos::ScalarTraits<ScalarType> SCT;
154  typedef typename SCT::magnitudeType MagnitudeType;
155 
157 
158 
160  PseudoBlockTFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
161  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
162  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
163  Teuchos::ParameterList &params );
164 
166  virtual ~PseudoBlockTFQMRIter() {};
168 
169 
171 
172 
183  void iterate();
184 
207 
211  void initialize()
212  {
214  initializeTFQMR(empty);
215  }
216 
226 
227  // Copy over the vectors.
228  state.W = W_;
229  state.U = U_;
230  state.AU = AU_;
231  state.Rtilde = Rtilde_;
232  state.D = D_;
233  state.V = V_;
234 
235  // Copy over the scalars.
236  state.alpha = alpha_;
237  state.eta = eta_;
238  state.rho = rho_;
239  state.tau = tau_;
240  state.theta = theta_;
241 
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::PseudoBlockTFQMRIter::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  // Classes inputed through constructor that define the linear problem to be solved.
295  //
296  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
297  const Teuchos::RCP<OutputManager<ScalarType> > om_;
298  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
299 
300  //
301  // Algorithmic parameters
302  //
303 
304  // numRHS_ is the current number of linear systems being solved.
305  int numRHS_;
306 
307  // Storage for QR factorization of the least squares system.
308  std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
309  std::vector<MagnitudeType> tau_, theta_;
310 
311  //
312  // Current solver state
313  //
314  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
315  // is capable of running; _initialize is controlled by the initialize() member method
316  // For the implications of the state of initialized_, please see documentation for initialize()
318 
319  // Current subspace dimension, and number of iterations performed.
320  int iter_;
321 
322  //
323  // State Storage
324  //
325  Teuchos::RCP<MV> W_;
326  Teuchos::RCP<MV> U_, AU_;
327  Teuchos::RCP<MV> Rtilde_;
328  Teuchos::RCP<MV> D_;
329  Teuchos::RCP<MV> V_;
330  Teuchos::RCP<MV> solnUpdate_;
331 
332  };
333 
334 
335  //
336  // Implementation
337  //
338 
340  // Constructor.
341  template <class ScalarType, class MV, class OP>
343  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
344  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
345  Teuchos::ParameterList &/* params */
346  ) :
347  lp_(problem),
348  om_(printer),
349  stest_(tester),
350  numRHS_(0),
351  initialized_(false),
352  iter_(0)
353  {
354  }
355 
357  // Compute native residual from TFQMR recurrence.
358  template <class ScalarType, class MV, class OP>
359  Teuchos::RCP<const MV>
360  PseudoBlockTFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
361  {
362  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
363  if (normvec) {
364  // Resize the vector passed in, if it is too small.
365  if ((int) normvec->size() < numRHS_)
366  normvec->resize( numRHS_ );
367 
368  // Compute the native residuals.
369  for (int i=0; i<numRHS_; i++) {
370  (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
371  }
372  }
373 
374  return Teuchos::null;
375  }
376 
378  // Initialize this iteration object
379  template <class ScalarType, class MV, class OP>
381  {
382  // Create convenience variables for zero and one.
383  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
384  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
385  const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
386 
387  // NOTE: In PseudoBlockTFQMRIter Rtilde_, the initial residual, is required!!!
388  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Rtilde == Teuchos::null,std::invalid_argument,
389  "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
390 
391  // Get the number of right-hand sides we're solving for now.
392  int numRHS = MVT::GetNumberVecs(*newstate.Rtilde);
393  numRHS_ = numRHS;
394 
395  // Initialize the state storage
396  // If the subspace has not be initialized before or we are reusing this solver object, generate it using Rtilde.
397  if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
398  {
399  // Create and/or initialize D_.
400  if ( Teuchos::is_null(D_) )
401  D_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
402  MVT::MvInit( *D_, STzero );
403 
404  // Create and/or initialize solnUpdate_;
405  if ( Teuchos::is_null(solnUpdate_) )
406  solnUpdate_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
407  MVT::MvInit( *solnUpdate_, STzero );
408 
409  // Create Rtilde_.
410  if (newstate.Rtilde != Rtilde_)
411  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
412  W_ = MVT::CloneCopy( *Rtilde_ );
413  U_ = MVT::CloneCopy( *Rtilde_ );
414  V_ = MVT::Clone( *Rtilde_, numRHS_ );
415 
416  // Multiply the current residual by Op and store in V_
417  // V_ = Op * R_
418  lp_->apply( *U_, *V_ );
419  AU_ = MVT::CloneCopy( *V_ );
420 
421  // Resize work vectors.
422  alpha_.resize( numRHS_, STone );
423  eta_.resize( numRHS_, STzero );
424  rho_.resize( numRHS_ );
425  rho_old_.resize( numRHS_ );
426  tau_.resize( numRHS_ );
427  theta_.resize( numRHS_, MTzero );
428 
429  MVT::MvNorm( *Rtilde_, tau_ ); // tau = ||r_0||
430  MVT::MvDot( *Rtilde_, *Rtilde_, rho_ ); // rho = (r_tilde, r0)
431  }
432  else
433  {
434  // If the subspace has changed sizes, clone it from the incoming state.
435  Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
436  W_ = MVT::CloneCopy( *newstate.W );
437  U_ = MVT::CloneCopy( *newstate.U );
438  AU_ = MVT::CloneCopy( *newstate.AU );
439  D_ = MVT::CloneCopy( *newstate.D );
440  V_ = MVT::CloneCopy( *newstate.V );
441 
442  // The solution update is just set to zero, since the current update has already
443  // been added to the solution during deflation.
444  solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
445  MVT::MvInit( *solnUpdate_, STzero );
446 
447  // Copy work vectors.
448  alpha_ = newstate.alpha;
449  eta_ = newstate.eta;
450  rho_ = newstate.rho;
451  tau_ = newstate.tau;
452  theta_ = newstate.theta;
453  }
454 
455  // The solver is initialized
456  initialized_ = true;
457  }
458 
459 
461  // Iterate until the status test informs us we should stop.
462  template <class ScalarType, class MV, class OP>
464  {
465  //
466  // Allocate/initialize data structures
467  //
468  if (initialized_ == false) {
469  initialize();
470  }
471 
472  // Create convenience variables for zero and one.
473  const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
474  const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
475  const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
476  const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
477  std::vector< ScalarType > beta(numRHS_,STzero);
478  std::vector<int> index(1);
479  //
480  // Start executable statements.
481  //
482 
484  // Iterate until the status test tells us to stop.
485  //
486  while (stest_->checkStatus(this) != Passed) {
487 
488  for (int iIter=0; iIter<2; iIter++)
489  {
490  //
491  //--------------------------------------------------------
492  // Compute the new alpha if we need to
493  //--------------------------------------------------------
494  //
495  if (iIter == 0) {
496  MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
497  for (int i=0; i<numRHS_; i++) {
498  rho_old_[i] = rho_[i]; // rho_old = rho
499  alpha_[i] = rho_old_[i]/alpha_[i];
500  }
501  }
502  //
503  //--------------------------------------------------------
504  // Loop over all RHS and compute updates.
505  //--------------------------------------------------------
506  //
507  for (int i=0; i<numRHS_; ++i) {
508  index[0] = i;
509 
510  //
511  //--------------------------------------------------------
512  // Update w.
513  // w = w - alpha*Au
514  //--------------------------------------------------------
515  //
516  Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
517  Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
518  MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
519  //
520  //--------------------------------------------------------
521  // Update d.
522  // d = u + (theta^2/alpha)eta*d
523  //--------------------------------------------------------
524  //
525  Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
526  Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
527  MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
528  //
529  //--------------------------------------------------------
530  // Update u if we need to.
531  // u = u - alpha*v
532  //
533  // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
534  //--------------------------------------------------------
535  //
536  if (iIter == 0) {
537  // Compute new U.
538  Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
539  Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
540  MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
541  }
542  }
543  //
544  //--------------------------------------------------------
545  // Update Au for the next iteration.
546  //--------------------------------------------------------
547  //
548  if (iIter == 0) {
549  lp_->apply( *U_, *AU_ );
550  }
551  //
552  //--------------------------------------------------------
553  // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
554  //--------------------------------------------------------
555  //
556  MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
557 
558  bool breakdown=false;
559  for (int i=0; i<numRHS_; ++i) {
560  theta_[i] /= tau_[i];
561  // cs = 1.0 / sqrt(1.0 + theta^2)
562  MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
563  tau_[i] *= theta_[i]*cs; // tau = tau * theta * cs
564  if ( tau_[i] == MTzero ) {
565  breakdown = true;
566  }
567  eta_[i] = cs*cs*alpha_[i]; // eta = cs^2 * alpha
568  }
569  //
570  //--------------------------------------------------------
571  // Accumulate the update for the solution x := x + eta*D_
572  //--------------------------------------------------------
573  //
574  for (int i=0; i<numRHS_; ++i) {
575  index[0]=i;
576  Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
577  Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
578  MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
579  }
580  //
581  //--------------------------------------------------------
582  // Breakdown was detected above, return to status test to
583  // remove converged solutions.
584  //--------------------------------------------------------
585  if ( breakdown ) {
586  break;
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
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
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.
std::vector< MagnitudeType > theta_
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.
std::vector< MagnitudeType > tau_
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockTFQMRIterateFailure(const std::string &what_arg)
const Teuchos::RCP< OutputManager< ScalarType > > om_
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.