Belos Version of the Day
Loading...
Searching...
No Matches
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
66#include "BelosStatusTest.hpp"
69
70#include "Teuchos_ScalarTraits.hpp"
71#include "Teuchos_ParameterList.hpp"
72#include "Teuchos_TimeMonitor.hpp"
73
85namespace 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 template <class ScalarType, class MV, class OP>
114 class PseudoBlockTFQMRIter : public Iteration<ScalarType,MV,OP> {
115 public:
116 //
117 // Convenience typedefs
118 //
121 typedef Teuchos::ScalarTraits<ScalarType> SCT;
122 typedef typename SCT::magnitudeType MagnitudeType;
123
125
126
128 PseudoBlockTFQMRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
129 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
130 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
131 Teuchos::ParameterList &params );
132
136
137
139
140
151 void iterate();
152
175
180 {
182 initializeTFQMR(empty);
183 }
184
194
195 // Copy over the vectors.
196 state.W = W_;
197 state.U = U_;
198 state.AU = AU_;
199 state.Rtilde = Rtilde_;
200 state.D = D_;
201 state.V = V_;
202
203 // Copy over the scalars.
204 state.alpha = alpha_;
205 state.eta = eta_;
206 state.rho = rho_;
207 state.tau = tau_;
208 state.theta = theta_;
209
210 return state;
211 }
212
214
215
217
218
220 int getNumIters() const { return iter_; }
221
223 void resetNumIters( int iter = 0 ) { iter_ = iter; }
224
227 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
228
230
233 Teuchos::RCP<MV> getCurrentUpdate() const { return solnUpdate_; }
234
236
237
239
240
242 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
243
245 int getBlockSize() const { return 1; }
246
248 void setBlockSize(int blockSize) {
249 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
250 "Belos::PseudoBlockTFQMRIter::setBlockSize(): Cannot use a block size that is not one.");
251 }
252
254 bool isInitialized() { return initialized_; }
255
257
258
259 private:
260
261 //
262 // Classes inputed through constructor that define the linear problem to be solved.
263 //
264 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
265 const Teuchos::RCP<OutputManager<ScalarType> > om_;
266 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
267
268 //
269 // Algorithmic parameters
270 //
271
272 // numRHS_ is the current number of linear systems being solved.
273 int numRHS_;
274
275 // Storage for QR factorization of the least squares system.
276 std::vector<ScalarType> alpha_, eta_, rho_, rho_old_;
277 std::vector<MagnitudeType> tau_, theta_;
278
279 //
280 // Current solver state
281 //
282 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
283 // is capable of running; _initialize is controlled by the initialize() member method
284 // For the implications of the state of initialized_, please see documentation for initialize()
285 bool initialized_;
286
287 // Current subspace dimension, and number of iterations performed.
288 int iter_;
289
290 //
291 // State Storage
292 //
293 Teuchos::RCP<MV> W_;
294 Teuchos::RCP<MV> U_, AU_;
295 Teuchos::RCP<MV> Rtilde_;
296 Teuchos::RCP<MV> D_;
297 Teuchos::RCP<MV> V_;
298 Teuchos::RCP<MV> solnUpdate_;
299
300 };
301
302
303 //
304 // Implementation
305 //
306
308 // Constructor.
309 template <class ScalarType, class MV, class OP>
311 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
312 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
313 Teuchos::ParameterList &/* params */
314 ) :
315 lp_(problem),
316 om_(printer),
317 stest_(tester),
318 numRHS_(0),
319 initialized_(false),
320 iter_(0)
321 {
322 }
323
325 // Compute native residual from TFQMR recurrence.
326 template <class ScalarType, class MV, class OP>
327 Teuchos::RCP<const MV>
328 PseudoBlockTFQMRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *normvec ) const
329 {
330 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
331 if (normvec) {
332 // Resize the vector passed in, if it is too small.
333 if ((int) normvec->size() < numRHS_)
334 normvec->resize( numRHS_ );
335
336 // Compute the native residuals.
337 for (int i=0; i<numRHS_; i++) {
338 (*normvec)[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( 2*iter_ + one )*tau_[i];
339 }
340 }
341
342 return Teuchos::null;
343 }
344
346 // Initialize this iteration object
347 template <class ScalarType, class MV, class OP>
349 {
350 // Create convenience variables for zero and one.
351 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
352 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
353 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
354
355 // NOTE: In PseudoBlockTFQMRIter Rtilde_, the initial residual, is required!!!
356 TEUCHOS_TEST_FOR_EXCEPTION(newstate.Rtilde == Teuchos::null,std::invalid_argument,
357 "Belos::PseudoBlockTFQMRIter::initialize(): PseudoBlockTFQMRIterState does not have initial residual.");
358
359 // Get the number of right-hand sides we're solving for now.
360 int numRHS = MVT::GetNumberVecs(*newstate.Rtilde);
361 numRHS_ = numRHS;
362
363 // Initialize the state storage
364 // If the subspace has not be initialized before or we are reusing this solver object, generate it using Rtilde.
365 if ( Teuchos::is_null(Rtilde_) || (MVT::GetNumberVecs(*Rtilde_) == numRHS_) )
366 {
367 // Create and/or initialize D_.
368 if ( Teuchos::is_null(D_) )
369 D_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
370 MVT::MvInit( *D_, STzero );
371
372 // Create and/or initialize solnUpdate_;
373 if ( Teuchos::is_null(solnUpdate_) )
374 solnUpdate_ = MVT::Clone( *newstate.Rtilde, numRHS_ );
375 MVT::MvInit( *solnUpdate_, STzero );
376
377 // Create Rtilde_.
378 if (newstate.Rtilde != Rtilde_)
379 Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
380 W_ = MVT::CloneCopy( *Rtilde_ );
381 U_ = MVT::CloneCopy( *Rtilde_ );
382 V_ = MVT::Clone( *Rtilde_, numRHS_ );
383
384 // Multiply the current residual by Op and store in V_
385 // V_ = Op * R_
386 lp_->apply( *U_, *V_ );
387 AU_ = MVT::CloneCopy( *V_ );
388
389 // Resize work vectors.
390 alpha_.resize( numRHS_, STone );
391 eta_.resize( numRHS_, STzero );
392 rho_.resize( numRHS_ );
393 rho_old_.resize( numRHS_ );
394 tau_.resize( numRHS_ );
395 theta_.resize( numRHS_, MTzero );
396
397 MVT::MvNorm( *Rtilde_, tau_ ); // tau = ||r_0||
398 MVT::MvDot( *Rtilde_, *Rtilde_, rho_ ); // rho = (r_tilde, r0)
399 }
400 else
401 {
402 // If the subspace has changed sizes, clone it from the incoming state.
403 Rtilde_ = MVT::CloneCopy( *newstate.Rtilde );
404 W_ = MVT::CloneCopy( *newstate.W );
405 U_ = MVT::CloneCopy( *newstate.U );
406 AU_ = MVT::CloneCopy( *newstate.AU );
407 D_ = MVT::CloneCopy( *newstate.D );
408 V_ = MVT::CloneCopy( *newstate.V );
409
410 // The solution update is just set to zero, since the current update has already
411 // been added to the solution during deflation.
412 solnUpdate_ = MVT::Clone( *Rtilde_, numRHS_ );
413 MVT::MvInit( *solnUpdate_, STzero );
414
415 // Copy work vectors.
416 alpha_ = newstate.alpha;
417 eta_ = newstate.eta;
418 rho_ = newstate.rho;
419 tau_ = newstate.tau;
420 theta_ = newstate.theta;
421 }
422
423 // The solver is initialized
424 initialized_ = true;
425 }
426
427
429 // Iterate until the status test informs us we should stop.
430 template <class ScalarType, class MV, class OP>
432 {
433 //
434 // Allocate/initialize data structures
435 //
436 if (initialized_ == false) {
437 initialize();
438 }
439
440 // Create convenience variables for zero and one.
441 const ScalarType STone = Teuchos::ScalarTraits<ScalarType>::one();
442 const ScalarType STzero = Teuchos::ScalarTraits<ScalarType>::zero();
443 const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
444 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
445 std::vector< ScalarType > beta(numRHS_,STzero);
446 std::vector<int> index(1);
447 //
448 // Start executable statements.
449 //
450
452 // Iterate until the status test tells us to stop.
453 //
454 while (stest_->checkStatus(this) != Passed) {
455
456 for (int iIter=0; iIter<2; iIter++)
457 {
458 //
459 //--------------------------------------------------------
460 // Compute the new alpha if we need to
461 //--------------------------------------------------------
462 //
463 if (iIter == 0) {
464 MVT::MvDot( *V_, *Rtilde_, alpha_ ); // alpha = rho / (r_tilde, v)
465 for (int i=0; i<numRHS_; i++) {
466 rho_old_[i] = rho_[i]; // rho_old = rho
467 alpha_[i] = rho_old_[i]/alpha_[i];
468 }
469 }
470 //
471 //--------------------------------------------------------
472 // Loop over all RHS and compute updates.
473 //--------------------------------------------------------
474 //
475 for (int i=0; i<numRHS_; ++i) {
476 index[0] = i;
477
478 //
479 //--------------------------------------------------------
480 // Update w.
481 // w = w - alpha*Au
482 //--------------------------------------------------------
483 //
484 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
485 Teuchos::RCP<MV> W_i = MVT::CloneViewNonConst( *W_, index );
486 MVT::MvAddMv( STone, *W_i, -alpha_[i], *AU_i, *W_i );
487 //
488 //--------------------------------------------------------
489 // Update d.
490 // d = u + (theta^2/alpha)eta*d
491 //--------------------------------------------------------
492 //
493 Teuchos::RCP<const MV> U_i = MVT::CloneView( *U_, index );
494 Teuchos::RCP<MV> D_i = MVT::CloneViewNonConst( *D_, index );
495 MVT::MvAddMv( STone, *U_i, (theta_[i]*theta_[i]/alpha_[i])*eta_[i], *D_i, *D_i );
496 //
497 //--------------------------------------------------------
498 // Update u if we need to.
499 // u = u - alpha*v
500 //
501 // Note: This is usually computed with alpha (above), but we're trying be memory efficient.
502 //--------------------------------------------------------
503 //
504 if (iIter == 0) {
505 // Compute new U.
506 Teuchos::RCP<const MV> V_i = MVT::CloneView( *V_, index );
507 Teuchos::RCP<MV> U2_i = MVT::CloneViewNonConst( *U_, index );
508 MVT::MvAddMv( STone, *U2_i, -alpha_[i], *V_i, *U2_i );
509 }
510 }
511 //
512 //--------------------------------------------------------
513 // Update Au for the next iteration.
514 //--------------------------------------------------------
515 //
516 if (iIter == 0) {
517 lp_->apply( *U_, *AU_ );
518 }
519 //
520 //--------------------------------------------------------
521 // Compute the new theta, c, eta, tau; i.e. the update to the least squares solution.
522 //--------------------------------------------------------
523 //
524 MVT::MvNorm( *W_, theta_ ); // theta = ||w|| / tau
525
526 bool breakdown=false;
527 for (int i=0; i<numRHS_; ++i) {
528 theta_[i] /= tau_[i];
529 // cs = 1.0 / sqrt(1.0 + theta^2)
530 MagnitudeType cs = MTone / Teuchos::ScalarTraits<MagnitudeType>::squareroot(MTone + theta_[i]*theta_[i]);
531 tau_[i] *= theta_[i]*cs; // tau = tau * theta * cs
532 if ( tau_[i] == MTzero ) {
533 breakdown = true;
534 }
535 eta_[i] = cs*cs*alpha_[i]; // eta = cs^2 * alpha
536 }
537 //
538 //--------------------------------------------------------
539 // Accumulate the update for the solution x := x + eta*D_
540 //--------------------------------------------------------
541 //
542 for (int i=0; i<numRHS_; ++i) {
543 index[0]=i;
544 Teuchos::RCP<const MV> D_i = MVT::CloneView( *D_, index );
545 Teuchos::RCP<MV> update_i = MVT::CloneViewNonConst( *solnUpdate_, index );
546 MVT::MvAddMv( STone, *update_i, eta_[i], *D_i, *update_i );
547 }
548 //
549 //--------------------------------------------------------
550 // Breakdown was detected above, return to status test to
551 // remove converged solutions.
552 //--------------------------------------------------------
553 if ( breakdown ) {
554 break;
555 }
556 //
557 if (iIter == 1) {
558 //
559 //--------------------------------------------------------
560 // Compute the new rho, beta if we need to.
561 //--------------------------------------------------------
562 //
563 MVT::MvDot( *W_, *Rtilde_, rho_ ); // rho = (r_tilde, w)
564
565 for (int i=0; i<numRHS_; ++i) {
566 beta[i] = rho_[i]/rho_old_[i]; // beta = rho / rho_old
567
568 //
569 //--------------------------------------------------------
570 // Update u, v, and Au if we need to.
571 // Note: We are updating v in two stages to be memory efficient
572 //--------------------------------------------------------
573 //
574 index[0]=i;
575 Teuchos::RCP<const MV> W_i = MVT::CloneView( *W_, index );
576 Teuchos::RCP<MV> U_i = MVT::CloneViewNonConst( *U_, index );
577 MVT::MvAddMv( STone, *W_i, beta[i], *U_i, *U_i ); // u = w + beta*u
578
579 // First stage of v update.
580 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
581 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
582 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
583 }
584
585 // Update Au.
586 lp_->apply( *U_, *AU_ ); // Au = A*u
587
588 // Second stage of v update.
589 for (int i=0; i<numRHS_; ++i) {
590 index[0]=i;
591 Teuchos::RCP<const MV> AU_i = MVT::CloneView( *AU_, index );
592 Teuchos::RCP<MV> V_i = MVT::CloneViewNonConst( *V_, index );
593 MVT::MvAddMv( STone, *AU_i, beta[i], *V_i, *V_i ); // v = Au + beta*v
594 }
595 }
596 }
597
598 // Increment the iteration
599 iter_++;
600
601 } // end while (sTest_->checkStatus(this) != Passed)
602 }
603
604} // namespace Belos
605//
606#endif // BELOS_PSEUDO_BLOCK_TFQMR_ITER_HPP
607//
608// End of file BelosPseudoBlockTFQMRIter.hpp
609
610
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.
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...
This class implements the preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
int getNumIters() const
Get the current iteration count.
void iterate()
This method performs pseudo-block TFQMR iterations until the status test indicates the need to stop o...
void resetNumIters(int iter=0)
Reset the iteration count.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
virtual ~PseudoBlockTFQMRIter()
Belos::PseudoBlockTFQMRIter destructor.
PseudoBlockTFQMRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
MultiVecTraits< ScalarType, MV > MVT
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.
Teuchos::ScalarTraits< ScalarType > SCT
void initializeTFQMR(const PseudoBlockTFQMRIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void setBlockSize(int blockSize)
Set the blocksize.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
OperatorTraits< ScalarType, MV, OP > OPT
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
Teuchos::RCP< const MV > W
The current residual basis.
Teuchos::ScalarTraits< ScalarType > SCT

Generated for Belos by doxygen 1.9.6