Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
66#include "BelosStatusTest.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
87namespace Belos {
88
93 template <class ScalarType, class MV>
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
119 class TFQMRIterateFailure : public BelosError {public:
120 TFQMRIterateFailure(const std::string& what_arg) : BelosError(what_arg)
121 {}};
122
124
125
126 template <class ScalarType, class MV, class OP>
127 class TFQMRIter : public Iteration<ScalarType,MV,OP> {
128 public:
129 //
130 // Convenience typedefs
131 //
134 typedef Teuchos::ScalarTraits<ScalarType> SCT;
135 typedef typename SCT::magnitudeType MagnitudeType;
136
138
139
141 TFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
142 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
143 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
144 Teuchos::ParameterList &params );
145
147 virtual ~TFQMRIter() {};
149
150
152
153
164 void iterate();
165
187 void initializeTFQMR(const TFQMRIterState<ScalarType,MV> & newstate);
188
193 {
195 initializeTFQMR(empty);
196 }
197
207 state.R = R_;
208 state.W = W_;
209 state.U = U_;
210 state.Rtilde = Rtilde_;
211 state.D = D_;
212 state.V = V_;
213 state.solnUpdate = solnUpdate_;
214 return state;
215 }
216
218
219
221
222
224 int getNumIters() const { return iter_; }
225
227 void resetNumIters( int iter = 0 ) { iter_ = iter; }
228
231 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
232
234
237 Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
238
240
241
243
244
247
249 int getBlockSize() const { return 1; }
250
252 void setBlockSize(int blockSize) {
253 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
254 "Belos::TFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
255 }
256
258 bool isInitialized() { return initialized_; }
259
261
262
263 private:
264
265 //
266 // Internal methods
267 //
269 void setStateSize();
270
271 //
272 // Classes inputed through constructor that define the linear problem to be solved.
273 //
274 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
275 const Teuchos::RCP<OutputManager<ScalarType> > om_;
276 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
277
278 //
279 // Algorithmic parameters
280 //
281
282 // Storage for QR factorization of the least squares system.
283 // Teuchos::SerialDenseMatrix<int,ScalarType> alpha_, rho_, rho_old_;
284 std::vector<ScalarType> alpha_, rho_, rho_old_;
285 std::vector<MagnitudeType> tau_, cs_, theta_;
286
287 //
288 // Current solver state
289 //
290 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
291 // is capable of running; _initialize is controlled by the initialize() member method
292 // For the implications of the state of initialized_, please see documentation for initialize()
294
295 // stateStorageInitialized_ specifies that the state storage has be initialized to the current
296 // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
297 // generated without the right-hand side or solution vectors.
299
300 // Current subspace dimension, and number of iterations performed.
301 int iter_;
302
303 //
304 // State Storage
305 //
306 Teuchos::RCP<MV> R_;
307 Teuchos::RCP<MV> W_;
308 Teuchos::RCP<MV> U_, AU_;
309 Teuchos::RCP<MV> Rtilde_;
310 Teuchos::RCP<MV> D_;
311 Teuchos::RCP<MV> V_;
312 Teuchos::RCP<MV> solnUpdate_;
313 };
314
315
316 //
317 // Implementation
318 //
319
321 // Constructor.
322 template <class ScalarType, class MV, class OP>
324 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
325 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
326 Teuchos::ParameterList &/* params */
327 ) :
328 lp_(problem),
329 om_(printer),
330 stest_(tester),
331 alpha_(1),
332 rho_(1),
333 rho_old_(1),
334 tau_(1),
335 cs_(1),
336 theta_(1),
337 initialized_(false),
338 stateStorageInitialized_(false),
339 iter_(0)
340 {
341 }
342
344 // Compute native residual from TFQMR recurrence.
345 template <class ScalarType, class MV, class OP>
346 Teuchos::RCP<const MV>
347 TFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
348 {
349 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
350 if (normvec)
351 (*normvec)[0] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[0];
352
353 return Teuchos::null;
354 }
355
356
358 // Setup the state storage.
359 template <class ScalarType, class MV, class OP>
361 {
362 if (!stateStorageInitialized_) {
363
364 // Check if there is any multivector to clone from.
365 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
366 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
367 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
368 stateStorageInitialized_ = false;
369 return;
370 }
371 else {
372
373 // Initialize the state storage
374 // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
375 if (R_ == Teuchos::null) {
376 // Get the multivector that is not null.
377 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
378 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
379 "Belos::TFQMRIter::setStateSize(): linear problem does not specify multivectors to clone from.");
380 R_ = MVT::Clone( *tmp, 1 );
381 D_ = MVT::Clone( *tmp, 1 );
382 V_ = MVT::Clone( *tmp, 1 );
383 solnUpdate_ = MVT::Clone( *tmp, 1 );
384 }
385
386 // State storage has now been initialized.
387 stateStorageInitialized_ = true;
388 }
389 }
390 }
391
393 // Initialize this iteration object
394 template <class ScalarType, class MV, class OP>
396 {
397 // Initialize the state storage if it isn't already.
398 if (!stateStorageInitialized_)
399 setStateSize();
400
401 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
402 "Belos::TFQMRIter::initialize(): Cannot initialize state storage!");
403
404 // NOTE: In TFQMRIter R_, the initial residual, is required!!!
405 //
406 std::string errstr("Belos::TFQMRIter::initialize(): Specified multivectors must have a consistent length and width.");
407
408 // Create convenience variables for zero and one.
409 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
410
411 if (newstate.R != Teuchos::null) {
412
413 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
414 std::invalid_argument, errstr );
415 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
416 std::invalid_argument, errstr );
417
418 // Copy basis vectors from newstate into V
419 if (newstate.R != R_) {
420 // copy over the initial residual (unpreconditioned).
421 MVT::Assign( *newstate.R, *R_ );
422 }
423
424 // Compute initial vectors
425 // Initially, they are set to the preconditioned residuals
426 //
427 W_ = MVT::CloneCopy( *R_ );
428 U_ = MVT::CloneCopy( *R_ );
429 Rtilde_ = MVT::CloneCopy( *R_ );
430 MVT::MvInit( *D_ );
431 MVT::MvInit( *solnUpdate_ );
432 // Multiply the current residual by Op and store in V_
433 // V_ = Op * R_
434 //
435 lp_->apply( *U_, *V_ );
436 AU_ = MVT::CloneCopy( *V_ );
437 //
438 // Compute initial scalars: theta, eta, tau, rho_old
439 //
440 theta_[0] = MTzero;
441 MVT::MvNorm( *R_, tau_ ); // tau = ||r_0||
442 MVT::MvDot( *R_, *Rtilde_, rho_old_ ); // rho = (r_tilde, r0)
443 }
444 else {
445
446 TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
447 "Belos::TFQMRIter::initialize(): TFQMRIterState does not have initial residual.");
448 }
449
450 // The solver is initialized
451 initialized_ = true;
452 }
453
454
456 // Iterate until the status test informs us we should stop.
457 template <class ScalarType, class MV, class OP>
459 {
460 //
461 // Allocate/initialize data structures
462 //
463 if (initialized_ == false) {
464 initialize();
465 }
466
467 // Create convenience variables for zero and one.
468 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
469 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
470 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
471 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
472 ScalarType eta = STzero, beta = STzero;
473 //
474 // Start executable statements.
475 //
476 // Get the current solution vector.
477 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
478
479 // Check that the current solution vector only has one column.
480 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, TFQMRIterateFailure,
481 "Belos::TFQMRIter::iterate(): current linear system has more than one vector!" );
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 alpha_[0] = rho_old_[0]/alpha_[0];
499 }
500 //
501 //--------------------------------------------------------
502 // Update w.
503 // w = w - alpha*Au
504 //--------------------------------------------------------
505 //
506 MVT::MvAddMv( STone, *W_, -alpha_[0], *AU_, *W_ );
507 //
508 //--------------------------------------------------------
509 // Update d.
510 // d = u + (theta^2/alpha)eta*d
511 //--------------------------------------------------------
512 //
513 MVT::MvAddMv( STone, *U_, (theta_[0]*theta_[0]/alpha_[0])*eta, *D_, *D_ );
514 //
515 //--------------------------------------------------------
516 // Update u if we need to.
517 // u = u - alpha*v
518 //
519 // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
520 //--------------------------------------------------------
521 //
522 if (iIter == 0) {
523 // Compute new U.
524 MVT::MvAddMv( STone, *U_, -alpha_[0], *V_, *U_ );
525
526 // Update Au for the next iteration.
527 lp_->apply( *U_, *AU_ );
528 }
529 //
530 //--------------------------------------------------------
531 // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
532 //--------------------------------------------------------
533 //
534 MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
535 theta_[0] /= tau_[0];
536 // cs = 1.0 / sqrt(1.0 + theta^2)
537 cs_[0] = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[0]*theta_[0]);
538 tau_[0] *= theta_[0]*cs_[0]; // tau = tau * theta * cs
539 eta = cs_[0]*cs_[0]*alpha_[0]; // eta = cs^2 * alpha
540 //
541 //--------------------------------------------------------
542 // Update the solution.
543 // Don't update the linear problem object, may incur additional preconditioner application.
544 //--------------------------------------------------------
545 //
546 MVT::MvAddMv( STone, *solnUpdate_, eta, *D_, *solnUpdate_ );
547 //
548 //--------------------------------------------------------
549 // Check for breakdown before continuing.
550 //--------------------------------------------------------
551 if ( tau_[0] == MTzero ) {
552 break;
553 }
554 //
555 if (iIter == 1) {
556 //
557 //--------------------------------------------------------
558 // Compute the new rho, beta if we need to.
559 //--------------------------------------------------------
560 //
561 MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
562 beta = rho_[0]/rho_old_[0]; // beta = rho / rho_old
563 rho_old_[0] = rho_[0]; // rho_old = rho
564 //
565 //--------------------------------------------------------
566 // Update u, v, and Au if we need to.
567 // Note: We are updating v in two stages to be memory efficient
568 //--------------------------------------------------------
569 //
570 MVT::MvAddMv( STone, *W_, beta, *U_, *U_ ); // u = w + beta*u
571
572 // First stage of v update.
573 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
574
575 // Update Au.
576 lp_->apply( *U_, *AU_ ); // Au = A*u
577
578 // Second stage of v update.
579 MVT::MvAddMv( STone, *AU_, beta, *V_, *V_ ); // v = Au + beta*v
580 }
581
582 }
583
584 // Increment the iteration
585 iter_++;
586
587 } // end while (sTest_->checkStatus(this) != Passed)
588 }
589
590} // namespace Belos
591//
592#endif // BELOS_TFQMR_ITER_HPP
593//
594// End of file BelosTFQMRIter.hpp
595
596
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which describes the basic interface to the linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
Teuchos::ScalarTraits< ScalarType > SCT
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< MV > D_
std::vector< MagnitudeType > theta_
Teuchos::RCP< MV > Rtilde_
virtual ~TFQMRIter()
Belos::TFQMRIter destructor.
void setBlockSize(int blockSize)
Set the blocksize.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
std::vector< ScalarType > rho_
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< MV > V_
void iterate()
This method performs TFQMR iterations until the status test indicates the need to stop or an error oc...
Teuchos::RCP< MV > R_
Teuchos::RCP< MV > W_
SCT::magnitudeType MagnitudeType
std::vector< ScalarType > rho_old_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< MV > U_
TFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > solnUpdate_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< MV > AU_
void setStateSize()
Method for initalizing the state storage needed by TFQMR.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MultiVecTraits< ScalarType, MV > MVT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
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.
int getNumIters() const
Get the current iteration count.
std::vector< MagnitudeType > cs_
std::vector< MagnitudeType > tau_
std::vector< ScalarType > alpha_
void initializeTFQMR(const TFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
TFQMRIterateFailure is thrown when the TFQMRIter object is unable to compute the next iterate in the ...
TFQMRIterateFailure(const std::string &what_arg)
Structure to contain pointers to TFQMRIter state variables.
Teuchos::RCP< const MV > W
Teuchos::RCP< const MV > V
Teuchos::RCP< const MV > Rtilde
Teuchos::RCP< const MV > R
The current residual basis.
Teuchos::RCP< const MV > D
Teuchos::RCP< const MV > U