Belos Version of the Day
Loading...
Searching...
No Matches
BelosMinresIter.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#ifndef BELOS_MINRES_ITER_HPP
43#define BELOS_MINRES_ITER_HPP
44
61
62#include "BelosConfigDefs.hpp"
63#include "BelosTypes.hpp"
65
68#include "BelosStatusTest.hpp"
71
72#include "Teuchos_SerialDenseMatrix.hpp"
73#include "Teuchos_SerialDenseVector.hpp"
74#include "Teuchos_ScalarTraits.hpp"
75#include "Teuchos_ParameterList.hpp"
76#include "Teuchos_TimeMonitor.hpp"
77
78namespace Belos {
79
93template<class ScalarType, class MV, class OP>
94class MinresIter : virtual public MinresIteration<ScalarType,MV,OP> {
95
96 public:
97
98 //
99 // Convenience typedefs
100 //
103 typedef Teuchos::ScalarTraits< ScalarType > SCT;
104 typedef typename SCT::magnitudeType MagnitudeType;
105 typedef Teuchos::ScalarTraits< MagnitudeType > SMT;
106
108
109
118 MinresIter (const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > >& problem,
119 const Teuchos::RCP< OutputManager< ScalarType > > & printer,
120 const Teuchos::RCP< StatusTest< ScalarType, MV, OP > >& tester,
121 const Teuchos::ParameterList& params);
122
124 virtual ~MinresIter() {};
126
127
129
130
145 void iterate();
146
162
169 {
171 initializeMinres(empty);
172 }
173
181 if (! isInitialized())
182 throw std::logic_error("getState() cannot be called unless "
183 "the state has been initialized");
185 state.Y = Y_;
186 state.R1 = R1_;
187 state.R2 = R2_;
188 state.W = W_;
189 state.W1 = W1_;
190 state.W2 = W2_;
191 return state;
192 }
193
195
196
198
199
201 int getNumIters() const { return iter_; }
202
204 void resetNumIters( int iter = 0 ) { iter_ = iter; }
205
208 Teuchos::RCP<const MV>
209 getNativeResiduals( std::vector<MagnitudeType> *norms ) const
210 {
211 if (norms != NULL)
212 {
213 std::vector<MagnitudeType>& theNorms = *norms;
214 if (theNorms.size() < 1)
215 theNorms.resize(1);
216 theNorms[0] = phibar_;
217 }
218 return Teuchos::null;
219 }
220
222
224 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
225
227 void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
228
230
232
233
235 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
236
238 int getBlockSize() const { return 1; }
239
241 void setBlockSize(int blockSize) {
242 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
243 "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
244 }
245
247 bool isInitialized() const { return initialized_; }
248 bool isInitialized() { return initialized_; }
249
251
252 private:
253
254 //
255 // Internal methods
256 //
258 void setStateSize();
259
260 //
261 // Classes inputed through constructor that define the linear problem to be solved.
262 //
263 const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_;
264 const Teuchos::RCP< OutputManager< ScalarType > > om_;
265 const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_;
266
267
275 bool initialized_;
276
283 bool stateStorageInitialized_;
284
286 int iter_;
287
292 MagnitudeType phibar_;
293
294 //
295 // State Storage
296 //
297
299 Teuchos::RCP< MV > Y_;
301 Teuchos::RCP< MV > R1_;
303 Teuchos::RCP< MV > R2_;
305 Teuchos::RCP< MV > W_;
307 Teuchos::RCP< MV > W1_;
309 Teuchos::RCP< MV > W2_;
310
318 Teuchos::SerialDenseMatrix<int,ScalarType> beta1_;
319
320};
321
323 // Constructor.
324 template<class ScalarType, class MV, class OP>
326 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
327 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
328 const Teuchos::ParameterList &/* params */ ):
329 lp_(problem),
330 om_(printer),
331 stest_(tester),
332 initialized_(false),
333 stateStorageInitialized_(false),
334 iter_(0),
335 phibar_(0.0)
336 {
337 }
338
340 // Setup the state storage.
341 template <class ScalarType, class MV, class OP>
343 {
344 if (!stateStorageInitialized_) {
345
346 // Check if there is any multivector to clone from.
347 Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
348 Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ = false;
351 return;
352 }
353 else {
354
355 // Initialize the state storage
356 // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
357 if (Y_ == Teuchos::null) {
358 // Get the multivector that is not null.
359 Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
360 TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
361 std::invalid_argument,
362 "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
363 Y_ = MVT::Clone( *tmp, 1 );
364 R1_ = MVT::Clone( *tmp, 1 );
365 R2_ = MVT::Clone( *tmp, 1 );
366 W_ = MVT::Clone( *tmp, 1 );
367 W1_ = MVT::Clone( *tmp, 1 );
368 W2_ = MVT::Clone( *tmp, 1 );
369 }
370 // State storage has now been initialized.
371 stateStorageInitialized_ = true;
372 }
373 }
374 }
375
376
378 // Initialize this iteration object
379 template <class ScalarType, class MV, class OP>
381 {
382 // Initialize the state storage if it isn't already.
383 if (!stateStorageInitialized_)
384 setStateSize();
385
386 TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
387 std::invalid_argument,
388 "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
389
390 TEUCHOS_TEST_FOR_EXCEPTION( newstate.Y == Teuchos::null,
391 std::invalid_argument,
392 "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
393
394 std::string errstr("Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
395 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.Y) != MVT::GetGlobalLength(*Y_),
396 std::invalid_argument,
397 errstr );
398 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.Y) != 1,
399 std::invalid_argument,
400 errstr );
401
402 // Create convenience variables for zero, one.
403 const ScalarType one = SCT::one();
404 const MagnitudeType m_zero = SMT::zero();
405
406 // Set up y and v for the first Lanczos vector v_1.
407 // y = beta1_ P' v1, where P = C**(-1).
408 // v is really P' v1.
409 MVT::Assign( *newstate.Y, *R2_ );
410 MVT::Assign( *newstate.Y, *R1_ );
411
412 // Initialize the W's to 0.
413 MVT::MvInit ( *W_ );
414 MVT::MvInit ( *W2_ );
415
416 if ( lp_->getLeftPrec() != Teuchos::null ) {
417 lp_->applyLeftPrec( *newstate.Y, *Y_ );
418 if ( lp_->getRightPrec() != Teuchos::null ) {
419 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
420 lp_->applyRightPrec( *tmp, *Y_ );
421 }
422 }
423 else if ( lp_->getRightPrec() != Teuchos::null ) {
424 lp_->applyRightPrec( *newstate.Y, *Y_ );
425 }
426 else {
427 if (newstate.Y != Y_) {
428 // copy over the initial residual (unpreconditioned).
429 MVT::Assign( *newstate.Y, *Y_ );
430 }
431 }
432
433 // beta1_ = b'*y;
434 beta1_ = Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 );
435 MVT::MvTransMv( one, *newstate.Y, *Y_, beta1_ );
436
437 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < m_zero,
438 std::invalid_argument,
439 "The preconditioner is not positive definite." );
440
441 if( SCT::magnitude(beta1_(0,0)) == m_zero )
442 {
443 // X = 0
444 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
445 MVT::MvInit( *cur_soln_vec );
446 }
447
448 beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
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 one = SCT::one();
469 const ScalarType zero = SCT::zero();
470 const MagnitudeType m_zero = SMT::zero();
471
472 // Allocate memory for scalars.
473 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
474 Teuchos::SerialDenseMatrix<int,ScalarType> beta( beta1_ );
475 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
476
477 // Initialize a few variables.
478 ScalarType oldBeta = zero;
479 ScalarType epsln = zero;
480 ScalarType cs = -one;
481 ScalarType sn = zero;
482 ScalarType dbar = zero;
483
484 // Declare a few others that will be initialized in the loop.
485 ScalarType oldeps;
486 ScalarType delta;
487 ScalarType gbar;
488 ScalarType phi;
489 ScalarType gamma;
490
491 // Allocate workspace.
492 Teuchos::RCP<MV> V = MVT::Clone( *Y_, 1 );
493 Teuchos::RCP<MV> tmpY, tmpW; // Not allocated, just used to transfer ownership.
494
495 // Get the current solution vector.
496 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
497
498 // Check that the current solution vector only has one column.
499 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
501 "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
502
504 // Iterate until the status test tells us to stop.
505 //
506 while (stest_->checkStatus(this) != Passed) {
507
508 // Increment the iteration
509 iter_++;
510
511 // Normalize previous vector.
512 // v = y / beta(0,0);
513 MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
514
515 // Apply operator.
516 lp_->applyOp (*V, *Y_);
517
518 if (iter_ > 1)
519 MVT::MvAddMv (one, *Y_, -beta(0,0)/oldBeta, *R1_, *Y_);
520
521 // alpha := dot(V, Y_)
522 MVT::MvTransMv (one, *V, *Y_, alpha);
523
524 // y := y - alpha/beta r2
525 MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
526
527 // r1 = r2;
528 // r2 = y;
529 tmpY = R1_;
530 R1_ = R2_;
531 R2_ = Y_;
532 Y_ = tmpY;
533
534 // apply preconditioner
535 if ( lp_->getLeftPrec() != Teuchos::null ) {
536 lp_->applyLeftPrec( *R2_, *Y_ );
537 if ( lp_->getRightPrec() != Teuchos::null ) {
538 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
539 lp_->applyRightPrec( *tmp, *Y_ );
540 }
541 }
542 else if ( lp_->getRightPrec() != Teuchos::null ) {
543 lp_->applyRightPrec( *R2_, *Y_ );
544 } // else "y = r2"
545 else {
546 MVT::Assign( *R2_, *Y_ );
547 }
548
549 // Get new beta.
550 oldBeta = beta(0,0);
551 MVT::MvTransMv( one, *R2_, *Y_, beta );
552
553 // Intercept beta <= 0.
554 //
555 // Note: we don't try to test for nonzero imaginary component of
556 // beta, because (a) it could be small and nonzero due to
557 // rounding error in computing the inner product, and (b) it's
558 // hard to tell how big "not small" should be, without computing
559 // some error bounds (for example, by modifying the linear
560 // algebra library to compute a posteriori rounding error bounds
561 // for the inner product, and then changing
562 // Belos::MultiVecTraits to make this information available).
563 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) < m_zero,
565 "Belos::MinresIter::iterate(): Encountered negative "
566 "value " << beta(0,0) << " for r2^H*M*r2 at itera"
567 "tion " << iter_ << ": MINRES cannot continue." );
568 beta(0,0) = SCT::squareroot( beta(0,0) );
569
570 // Apply previous rotation Q_{k-1} to get
571 //
572 // [delta_k epsln_{k+1}] = [cs sn][dbar_k 0 ]
573 // [gbar_k dbar_{k+1} ] [-sn cs][alpha_k beta_{k+1}].
574 //
575 oldeps = epsln;
576 delta = cs*dbar + sn*alpha(0,0);
577 gbar = sn*dbar - cs*alpha(0,0);
578 epsln = sn*beta(0,0);
579 dbar = - cs*beta(0,0);
580
581 // Compute the next plane rotation Q_k.
582 this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
583
584 phi = cs * phibar_; // phi_k
585 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ ); // phibar_{k+1}
586
587 // w1 = w2;
588 // w2 = w;
589 MVT::Assign( *W_, *W1_ );
590 tmpW = W1_;
591 W1_ = W2_;
592 W2_ = W_;
593 W_ = tmpW;
594
595 // w = (v - oldeps*w1 - delta*w2) / gamma;
596 MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
597 MVT::MvAddMv( one, *W_, -delta, *W2_, *W_ );
598 MVT::MvScale( *W_, one / gamma );
599
600 // Update x:
601 // x = x + phi*w;
602 MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
603 lp_->updateSolution();
604 } // end while (sTest_->checkStatus(this) != Passed)
605 }
606
607
609 // Compute the next plane rotation Qk.
610 // r = norm([a b]);
611 // c = a / r;
612 // s = b / r;
613 template <class ScalarType, class MV, class OP>
614 void MinresIter<ScalarType,MV,OP>::symOrtho( ScalarType a, ScalarType b,
615 ScalarType *c, ScalarType *s, ScalarType *r
616 )
617 {
618 const ScalarType one = SCT::one();
619 const ScalarType zero = SCT::zero();
620 const MagnitudeType m_zero = SMT::zero();
621 const MagnitudeType absA = SCT::magnitude( a );
622 const MagnitudeType absB = SCT::magnitude( b );
623 if ( absB == m_zero ) {
624 *s = zero;
625 *r = absA;
626 if ( absA == m_zero )
627 *c = one;
628 else
629 *c = a / absA;
630 } else if ( absA == m_zero ) {
631 *c = zero;
632 *s = b / absB;
633 *r = absB;
634 } else if ( absB >= absA ) { // && a!=0 && b!=0
635 ScalarType tau = a / b;
636 if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
637 *s = -one / SCT::squareroot( one+tau*tau );
638 else
639 *s = one / SCT::squareroot( one+tau*tau );
640 *c = *s * tau;
641 *r = b / *s;
642 } else { // (absA > absB) && a!=0 && b!=0
643 ScalarType tau = b / a;
644 if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
645 *c = -one / SCT::squareroot( one+tau*tau );
646 else
647 *c = one / SCT::squareroot( one+tau*tau );
648 *s = *c * tau;
649 *r = a / *c;
650 }
651 }
652
653} // end Belos namespace
654
655#endif /* BELOS_MINRES_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Pure virtual base class which augments the basic interface for a minimal residual linear solver itera...
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.
MINRES implementation.
void initializeMinres(const MinresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
MinresIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::ParameterList &params)
Constructor.
bool isInitialized()
States whether the solver has been initialized or not.
int getNumIters() const
Get the current iteration count.
bool isInitialized() const
States whether the solver has been initialized or not.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::ScalarTraits< MagnitudeType > SMT
void initialize()
Initialize the solver.
void symOrtho(ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r)
void iterate()
Perform MINRES iterations until convergence or error.
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
OperatorTraits< ScalarType, MV, OP > OPT
void resetNumIters(int iter=0)
Reset the iteration count.
virtual ~MinresIter()
Destructor.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MinresIterateFailure is thrown when the MinresIteration object is unable to compute the next iterate ...
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.
Structure to contain pointers to MinresIteration state variables.
Teuchos::RCP< const MV > W
The current direction vector.
Teuchos::RCP< const MV > R2
Previous residual.
Teuchos::RCP< const MV > W1
Previous direction vector.
Teuchos::RCP< const MV > Y
The current residual.
Teuchos::RCP< const MV > W2
Previous direction vector.
Teuchos::RCP< const MV > R1
Previous residual.

Generated for Belos by doxygen 1.10.0