Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosGCRODRIter.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_GCRODR_ITER_HPP
43#define BELOS_GCRODR_ITER_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51
55#include "BelosStatusTest.hpp"
58
59#include "Teuchos_BLAS.hpp"
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
65
78namespace Belos {
79
81
82
87 template <class ScalarType, class MV>
93 int curDim;
94
96 Teuchos::RCP<MV> V;
97
99 Teuchos::RCP<MV> U, C;
100
106 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H;
107
110 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B;
111
112 GCRODRIterState() : curDim(0), V(Teuchos::null),
113 U(Teuchos::null), C(Teuchos::null),
114 H(Teuchos::null), B(Teuchos::null)
115 {}
116 };
117
119
121
122
135 public:
136 GCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
137 };
138
146 public:
147 GCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
148 };
149
151
152
153 template<class ScalarType, class MV, class OP>
154 class GCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
155
156 public:
157
158 //
159 // Convenience typedefs
160 //
163 typedef Teuchos::ScalarTraits<ScalarType> SCT;
164 typedef typename SCT::magnitudeType MagnitudeType;
165
167
168
177 GCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
178 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
179 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
180 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
181 Teuchos::ParameterList &params );
182
184 virtual ~GCRODRIter() {};
186
187
189
190
212 void iterate();
213
236
240 void initialize() {
242 initialize(empty);
243 }
244
254 state.curDim = curDim_;
255 state.V = V_;
256 state.U = U_;
257 state.C = C_;
258 state.H = H_;
259 state.B = B_;
260 return state;
261 }
262
264
265
267
268
270 int getNumIters() const { return iter_; }
271
273 void resetNumIters( int iter = 0 ) { iter_ = iter; }
274
277 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
278
280
285 Teuchos::RCP<MV> getCurrentUpdate() const;
286
288
291 void updateLSQR( int dim = -1 );
292
294 int getCurSubspaceDim() const {
295 if (!initialized_) return 0;
296 return curDim_;
297 };
298
300 int getMaxSubspaceDim() const { return numBlocks_; }
301
303
304
306
307
310
312 int getNumBlocks() const { return numBlocks_; }
313
315 void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
316
318 int getRecycledBlocks() const { return recycledBlocks_; }
319
321 void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
322
324 int getBlockSize() const { return 1; }
325
327 void setBlockSize(int blockSize) {
328 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
329 }
330
332 void setSize( int recycledBlocks, int numBlocks ) {
333 // only call resize if size changed
334 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
335 recycledBlocks_ = recycledBlocks;
336 numBlocks_ = numBlocks;
337 cs_.sizeUninitialized( numBlocks_+1 );
338 sn_.sizeUninitialized( numBlocks_+1 );
339 z_.sizeUninitialized( numBlocks_+1 );
340 R_.shapeUninitialized( numBlocks_+1,numBlocks );
341 }
342 }
343
345 bool isInitialized() { return initialized_; }
346
348
349 private:
350
351 //
352 // Internal methods
353 //
354
355 //
356 // Classes inputed through constructor that define the linear problem to be solved.
357 //
358 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
359 const Teuchos::RCP<OutputManager<ScalarType> > om_;
360 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
361 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
362
363 //
364 // Algorithmic parameters
365 //
366 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
368
369 // recycledBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
371
372 // Storage for QR factorization of the least squares system.
373 Teuchos::SerialDenseVector<int,ScalarType> sn_;
374 Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
375
376 //
377 // Current solver state
378 //
379 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
380 // is capable of running; _initialize is controlled by the initialize() member method
381 // For the implications of the state of initialized_, please see documentation for initialize()
383
384 // Current subspace dimension, and number of iterations performed.
386
387 //
388 // State Storage
389 //
390 // Krylov vectors.
391 Teuchos::RCP<MV> V_;
392 //
393 // Recycled subspace vectors.
394 Teuchos::RCP<MV> U_, C_;
395 //
396 // Projected matrices
397 // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
398 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
399 //
400 // B_ : Projected matrix from the recycled subspace B = C^H*A*V
401 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
402 //
403 // QR decomposition of Projected matrices for solving the least squares system HY = B.
404 // R_: Upper triangular reduction of H
405 // z_: Q applied to right-hand side of the least squares system
406 Teuchos::SerialDenseMatrix<int,ScalarType> R_;
407 Teuchos::SerialDenseVector<int,ScalarType> z_;
408 };
409
411 // Constructor.
412 template<class ScalarType, class MV, class OP>
414 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
415 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
416 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
417 Teuchos::ParameterList &params ):
418 lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
419 numBlocks_ = 0;
420 recycledBlocks_ = 0;
421 initialized_ = false;
422 curDim_ = 0;
423 iter_ = 0;
424 V_ = Teuchos::null;
425 U_ = Teuchos::null;
426 C_ = Teuchos::null;
427 H_ = Teuchos::null;
428 B_ = Teuchos::null;
429
430 // Get the maximum number of blocks allowed for this Krylov subspace
431 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
432 int nb = Teuchos::getParameter<int>(params, "Num Blocks");
433
434 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
435 int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
436
437 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
438 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
439
440 numBlocks_ = nb;
441 recycledBlocks_ = rb;
442
443 cs_.sizeUninitialized( numBlocks_+1 );
444 sn_.sizeUninitialized( numBlocks_+1 );
445 z_.sizeUninitialized( numBlocks_+1 );
446 R_.shapeUninitialized( numBlocks_+1,numBlocks_ );
447
448 }
449
451 // Get the current update from this subspace.
452 template <class ScalarType, class MV, class OP>
454 //
455 // If this is the first iteration of the Arnoldi factorization,
456 // there is no update, so return Teuchos::null.
457 //
458 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
459 if (curDim_==0) {
460 return currentUpdate;
461 } else {
462 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
463 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464 Teuchos::BLAS<int,ScalarType> blas;
465 currentUpdate = MVT::Clone( *V_, 1 );
466 //
467 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
468 //
469 Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, z_, curDim_, 1 );
470 //
471 // Solve the least squares problem.
472 //
473 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
474 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
475 R_.values(), R_.stride(), y.values(), y.stride() );
476 //
477 // Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
478 //
479 std::vector<int> index(curDim_);
480 for ( int i=0; i<curDim_; i++ ) index[i] = i;
481 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
482 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
483 //
484 // Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
485 //
486 if (U_ != Teuchos::null) {
487 Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
488 Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
489 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
490 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
491 }
492 }
493 return currentUpdate;
494 }
495
496
498 // Get the native residuals stored in this iteration.
499 // Note: This method does not return a MultiVector of the residual vectors, only return Teuchos::null
500 template <class ScalarType, class MV, class OP>
501 Teuchos::RCP<const MV> GCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
502 //
503 // NOTE: Make sure the incoming std::vector is the correct size!
504 //
505 if ( norms && (int)norms->size()==0 )
506 norms->resize( 1 );
507
508 if (norms) {
509 Teuchos::BLAS<int,ScalarType> blas;
510 (*norms)[0] = blas.NRM2( 1, &z_(curDim_), 1);
511 }
512 return Teuchos::null;
513 }
514
515
516
518 // Initialize this iteration object
519 template <class ScalarType, class MV, class OP>
521
522 if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
523 curDim_ = newstate.curDim;
524 V_ = newstate.V;
525 U_ = newstate.U;
526 C_ = newstate.C;
527 H_ = newstate.H;
528 B_ = newstate.B;
529 }
530 else {
531 TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have V initialized.");
532 TEUCHOS_TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have H initialized.");
533 }
534
535 // the solver is initialized
536 initialized_ = true;
537
538 }
539
540
542 // Iterate until the status test informs us we should stop.
543 template <class ScalarType, class MV, class OP>
545
546 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, GCRODRIterInitFailure,"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
547
548 // Force call to setsize to ensure internal storage is correct dimension
549 setSize( recycledBlocks_, numBlocks_ );
550
551 Teuchos::RCP<MV> Vnext;
552 Teuchos::RCP<const MV> Vprev;
553 std::vector<int> curind(1);
554
555 // z_ must be zeroed out in order to compute Givens rotations correctly
556 z_.putScalar(0.0);
557
558 // Orthonormalize the new V_0
559 curind[0] = 0;
560 Vnext = MVT::CloneViewNonConst(*V_,curind);
561 // Orthonormalize first column
562 // First, get a monoelemental matrix to hold the orthonormalization coefficients
563 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z0 =
564 Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(1,1) );
565 int rank = ortho_->normalize( *Vnext, z0 );
566 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
567 // Copy the scalar coefficient back into the z_ vector
568 z_(0) = (*z0)(0,0);
569
570 std::vector<int> prevind(numBlocks_+1);
571
573 // iterate until the status test tells us to stop.
574 //
575 // also break if our basis is full
576 //
577 if (U_ == Teuchos::null) { // iterate without recycle space (because none is available yet)
578 while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
579 iter_++;
580 int lclDim = curDim_ + 1;
581
582 // Get next index in basis
583 curind[0] = lclDim;
584 Vnext = MVT::CloneViewNonConst(*V_,curind);
585
586 // Get previous index in basis
587 curind[0] = curDim_;
588 Vprev = MVT::CloneView(*V_,curind);
589
590 // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
591 lp_->apply(*Vprev,*Vnext);
592
593 // Remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.
594
595 // Get a view of all the previous vectors
596 prevind.resize(lclDim);
597 for (int i=0; i<lclDim; i++) { prevind[i] = i; }
598 Vprev = MVT::CloneView(*V_,prevind);
599 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
600
601 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
602 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
603 subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
604 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
605 AsubH.append( subH );
606
607 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
608 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
609 subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
610
611 // Project out the previous Krylov vectors and normalize the next vector.
612 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
613
614 // Copy over the coefficients to R just in case we run into an error.
615 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
616 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
617 subR2.assign(subH2);
618
619 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
620
621 // Update the QR factorization of the upper Hessenberg matrix
622 updateLSQR();
623
624 // Update basis dim
625 curDim_++;
626 } // end while (statusTest == false)
627 }
628 else { // iterate with recycle space
629 while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
630 iter_++;
631 int lclDim = curDim_ + 1;
632
633 // Get the current part of the basis.
634 curind[0] = lclDim;
635 Vnext = MVT::CloneViewNonConst(*V_,curind);
636
637 // Get a view of the previous vectors.
638 // This is used for orthogonalization and for computing V^H K H
639 curind[0] = curDim_;
640 Vprev = MVT::CloneView(*V_,curind);
641
642 // Compute the next std::vector in the Krylov basis: Vnext = Op*Vprev
643 lp_->apply(*Vprev,*Vnext);
644 Vprev = Teuchos::null;
645
646 // First, remove the recycled subspace (C) from Vnext and put coefficients in B.
647 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
648 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
649 subB = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
650 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
651 AsubB.append( subB );
652
653 // Project out the recycled subspace.
654 ortho_->project( *Vnext, AsubB, C );
655
656 // Now, remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.
657 // Get a view of all the previous vectors
658 prevind.resize(lclDim);
659 for (int i=0; i<lclDim; i++) { prevind[i] = i; }
660 Vprev = MVT::CloneView(*V_,prevind);
661 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
662
663 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
664 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
665 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
666 AsubH.append( subH );
667
668 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
669 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
670
671 // Project out the previous Krylov vectors and normalize the next vector.
672 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
673
674 // Copy over the coefficients to R just in case we run into an error.
675 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
676 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
677 subR2.assign(subH2);
678
679 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
680
681 // Update the QR factorization of the upper Hessenberg matrix
682 updateLSQR();
683
684 // Update basis dim
685 curDim_++;
686
687 } // end while (statusTest == false)
688 } // end if (U_ == Teuchos::null)
689
690 }
691
692
693 template<class ScalarType, class MV, class OP>
695
696 int i;
697 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
698
699 // Get correct dimension based on input "dim"
700 // Remember that ortho failures result in an exit before updateLSQR() is called.
701 // Therefore, it is possible that dim == curDim_.
702 int curDim = curDim_;
703 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
704 curDim = dim;
705
706 Teuchos::BLAS<int, ScalarType> blas;
707 //
708 // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
709 // system to upper-triangular form.
710 //
711 // QR factorization of Least-Squares system with Givens rotations
712 //
713 for (i=0; i<curDim; i++) {
714 //
715 // Apply previous Givens rotations to new column of Hessenberg matrix
716 //
717 blas.ROT( 1, &R_(i,curDim), 1, &R_(i+1, curDim), 1, &cs_[i], &sn_[i] );
718 }
719 //
720 // Calculate new Givens rotation
721 //
722 blas.ROTG( &R_(curDim,curDim), &R_(curDim+1,curDim), &cs_[curDim], &sn_[curDim] );
723 R_(curDim+1,curDim) = zero;
724 //
725 // Update RHS w/ new transformation
726 //
727 blas.ROT( 1, &z_(curDim), 1, &z_(curDim+1), 1, &cs_[curDim], &sn_[curDim] );
728
729 } // end updateLSQR()
730
731} // end Belos namespace
732
733#endif /* BELOS_GCRODR_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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
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.
GCRODRIterInitFailure is thrown when the GCRODRIter object is unable to generate an initial iterate i...
GCRODRIterInitFailure(const std::string &what_arg)
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
GCRODRIterOrthoFailure(const std::string &what_arg)
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
Teuchos::SerialDenseMatrix< int, ScalarType > R_
Teuchos::SerialDenseVector< int, ScalarType > sn_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
virtual ~GCRODRIter()
Destructor.
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::RCP< MV > C_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
Teuchos::SerialDenseVector< int, ScalarType > z_
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
OperatorTraits< ScalarType, MV, OP > OPT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void resetNumIters(int iter=0)
Reset the iteration count.
GCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
Teuchos::ScalarTraits< ScalarType > SCT
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
bool isInitialized()
States whether the solver has been initialized or not.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< MV > U_
Teuchos::RCP< MV > V_
int getNumIters() const
Get the current iteration count.
GCRODRIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
GCRODRIter constructor with linear problem, solver utilities, and parameter list of solver options.
void initialize()
Initialize the solver with empty data. Calling this method will result in error, as GCRODRIter must b...
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
SCT::magnitudeType MagnitudeType
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
Teuchos::SerialDenseVector< int, MagnitudeType > cs_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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 GCRODRIter state variables.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< MV > V
The current Krylov basis.
int curDim
The current dimension of the reduction.