Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosBlockGmresIter.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_BLOCK_GMRES_ITER_HPP
43#define BELOS_BLOCK_GMRES_ITER_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
52
56#include "BelosStatusTest.hpp"
59
60#include "Teuchos_BLAS.hpp"
61#include "Teuchos_LAPACK.hpp"
62#include "Teuchos_SerialDenseMatrix.hpp"
63#include "Teuchos_SerialDenseVector.hpp"
64#include "Teuchos_ScalarTraits.hpp"
65#include "Teuchos_ParameterList.hpp"
66#include "Teuchos_TimeMonitor.hpp"
67
81namespace Belos {
82
83template<class ScalarType, class MV, class OP>
84class BlockGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
85
86 public:
87
88 //
89 // Convenience typedefs
90 //
93 typedef Teuchos::ScalarTraits<ScalarType> SCT;
94 typedef typename SCT::magnitudeType MagnitudeType;
95
97
98
108 BlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
109 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
110 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
111 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
112 Teuchos::ParameterList &params );
113
115 virtual ~BlockGmresIter() {};
117
118
120
121
143 void iterate();
144
167
172 {
174 initializeGmres(empty);
175 }
176
186 state.curDim = curDim_;
187 state.V = V_;
188 state.H = H_;
189 state.R = R_;
190 state.z = z_;
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> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
209
211
216 Teuchos::RCP<MV> getCurrentUpdate() const;
217
219
222 void updateLSQR( int dim = -1 );
223
225 int getCurSubspaceDim() const {
226 if (!initialized_) return 0;
227 return curDim_;
228 };
229
232
234
235
237
238
241
243 int getBlockSize() const { return blockSize_; }
244
246 void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }
247
249 int getNumBlocks() const { return numBlocks_; }
250
252 void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }
253
260 void setSize(int blockSize, int numBlocks);
261
263 bool isInitialized() { return initialized_; }
264
266
267 private:
268
269 //
270 // Internal structs
271 //
272 struct CheckList {
273 bool checkV;
275 CheckList() : checkV(false), checkArn(false) {};
276 };
277 //
278 // Internal methods
279 //
281 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
282
284 void setStateSize();
285
286 //
287 // Classes inputed through constructor that define the linear problem to be solved.
288 //
289 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
290 const Teuchos::RCP<OutputManager<ScalarType> > om_;
291 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
292 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
293
294 //
295 // Algorithmic parameters
296 //
297 // blockSize_ is the solver block size.
298 // It controls the number of vectors added to the basis on each iteration.
300 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
302
303 // Storage for QR factorization of the least squares system.
304 Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
305 Teuchos::SerialDenseVector<int,MagnitudeType> cs;
306
307 //
308 // Current solver state
309 //
310 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
311 // is capable of running; _initialize is controlled by the initialize() member method
312 // For the implications of the state of initialized_, please see documentation for initialize()
314
315 // stateStorageInitialized_ specifies that the state storage has be initialized to the current
316 // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
317 // generated without the right-hand side or solution vectors.
319
320 // keepHessenberg_ specifies that the iteration must keep the Hessenberg matrix formed via the
321 // Arnoldi factorization and the upper triangular matrix that is the Hessenberg matrix reduced via
322 // QR factorization separate.
324
325 // initHessenberg_ specifies that the iteration should reinitialize the Hessenberg matrix by zeroing
326 // out all entries before an iteration is started.
328
329 // Current subspace dimension, and number of iterations performed.
331
332 //
333 // State Storage
334 //
335 Teuchos::RCP<MV> V_;
336 //
337 // Projected matrices
338 // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
339 //
340 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
341 //
342 // QR decomposition of Projected matrices for solving the least squares system HY = B.
343 // R_: Upper triangular reduction of H
344 // z_: Q applied to right-hand side of the least squares system
345 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
346 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;
347};
348
350 // Constructor.
351 template<class ScalarType, class MV, class OP>
353 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
354 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
355 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
356 Teuchos::ParameterList &params ):
357 lp_(problem),
358 om_(printer),
359 stest_(tester),
360 ortho_(ortho),
361 blockSize_(0),
362 numBlocks_(0),
363 initialized_(false),
364 stateStorageInitialized_(false),
365 keepHessenberg_(false),
366 initHessenberg_(false),
367 curDim_(0),
368 iter_(0)
369 {
370 // Find out whether we are saving the Hessenberg matrix.
371 if ( om_->isVerbosity( Debug ) )
372 keepHessenberg_ = true;
373 else
374 keepHessenberg_ = params.get("Keep Hessenberg", false);
375
376 // Find out whether we are initializing the Hessenberg matrix.
377 initHessenberg_ = params.get("Initialize Hessenberg", false);
378
379 // Get the maximum number of blocks allowed for this Krylov subspace
380 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
381 "Belos::BlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
382 int nb = Teuchos::getParameter<int>(params, "Num Blocks");
383
384 // Set the block size and allocate data
385 int bs = params.get("Block Size", 1);
386 setSize( bs, nb );
387 }
388
390 // Set the block size and make necessary adjustments.
391 template <class ScalarType, class MV, class OP>
392 void BlockGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
393 {
394 // This routine only allocates space; it doesn't not perform any computation
395 // any change in size will invalidate the state of the solver.
396
397 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setSize was passed a non-positive argument.");
398 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
399 // do nothing
400 return;
401 }
402
403 if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
404 stateStorageInitialized_ = false;
405
406 blockSize_ = blockSize;
407 numBlocks_ = numBlocks;
408
409 initialized_ = false;
410 curDim_ = 0;
411
412 // Use the current blockSize_ and numBlocks_ to initialize the state storage.
413 setStateSize();
414
415 }
416
418 // Setup the state storage.
419 template <class ScalarType, class MV, class OP>
421 {
422 if (!stateStorageInitialized_) {
423
424 // Check if there is any multivector to clone from.
425 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
426 Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
427 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
428 stateStorageInitialized_ = false;
429 return;
430 }
431 else {
432
434 // blockSize*numBlocks dependent
435 //
436 int newsd = blockSize_*(numBlocks_+1);
437
438 if (blockSize_==1) {
439 cs.resize( newsd );
440 sn.resize( newsd );
441 }
442 else {
443 beta.resize( newsd );
444 }
445
446 // Initialize the state storage
447 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_*static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
448 "Belos::BlockGmresIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
449
450 // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
451 if (V_ == Teuchos::null) {
452 // Get the multivector that is not null.
453 Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
454 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
455 "Belos::BlockGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
456 V_ = MVT::Clone( *tmp, newsd );
457 }
458 else {
459 // Generate V_ by cloning itself ONLY if more space is needed.
460 if (MVT::GetNumberVecs(*V_) < newsd) {
461 Teuchos::RCP<const MV> tmp = V_;
462 V_ = MVT::Clone( *tmp, newsd );
463 }
464 }
465
466 // Generate R_ only if it doesn't exist, otherwise resize it.
467 if (R_ == Teuchos::null) {
468 R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
469 }
470 if (initHessenberg_) {
471 R_->shape( newsd, newsd-blockSize_ );
472 }
473 else {
474 if (R_->numRows() < newsd || R_->numCols() < newsd-blockSize_) {
475 R_->shapeUninitialized( newsd, newsd-blockSize_ );
476 }
477 }
478
479 // Generate H_ only if it doesn't exist, and we are keeping the upper Hessenberg matrix.
480 if (keepHessenberg_) {
481 if (H_ == Teuchos::null) {
482 H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
483 }
484 if (initHessenberg_) {
485 H_->shape( newsd, newsd-blockSize_ );
486 }
487 else {
488 if (H_->numRows() < newsd || H_->numCols() < newsd-blockSize_) {
489 H_->shapeUninitialized( newsd, newsd-blockSize_ );
490 }
491 }
492 }
493 else {
494 // Point H_ and R_ at the same object.
495 H_ = R_;
496 }
497
498 // Generate z_ only if it doesn't exist, otherwise resize it.
499 if (z_ == Teuchos::null) {
500 z_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
501 }
502 if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
503 z_->shapeUninitialized( newsd, blockSize_ );
504 }
505
506 // State storage has now been initialized.
507 stateStorageInitialized_ = true;
508 }
509 }
510 }
511
513 // Get the current update from this subspace.
514 template <class ScalarType, class MV, class OP>
516 {
517 //
518 // If this is the first iteration of the Arnoldi factorization,
519 // there is no update, so return Teuchos::null.
520 //
521 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
522 if (curDim_==0) {
523 return currentUpdate;
524 } else {
525 const ScalarType one = SCT::one();
526 const ScalarType zero = SCT::zero();
527 Teuchos::BLAS<int,ScalarType> blas;
528 currentUpdate = MVT::Clone( *V_, blockSize_ );
529 //
530 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
531 //
532 Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
533 //
534 // Solve the least squares problem.
535 //
536 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
537 Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
538 R_->values(), R_->stride(), y.values(), y.stride() );
539 //
540 // Compute the current update.
541 //
542 std::vector<int> index(curDim_);
543 for ( int i=0; i<curDim_; i++ ) {
544 index[i] = i;
545 }
546 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
547 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
548 }
549 return currentUpdate;
550 }
551
552
554 // Get the native residuals stored in this iteration.
555 // Note: No residual std::vector will be returned by Gmres.
556 template <class ScalarType, class MV, class OP>
557 Teuchos::RCP<const MV> BlockGmresIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const
558 {
559 //
560 // NOTE: Make sure the incoming std::vector is the correct size!
561 //
562 if ( norms && (int)norms->size() < blockSize_ )
563 norms->resize( blockSize_ );
564
565 if (norms) {
566 Teuchos::BLAS<int,ScalarType> blas;
567 for (int j=0; j<blockSize_; j++) {
568 (*norms)[j] = blas.NRM2( blockSize_, &(*z_)(curDim_, j), 1);
569 }
570 }
571 return Teuchos::null;
572 }
573
574
575
577 // Initialize this iteration object
578 template <class ScalarType, class MV, class OP>
580 {
581 // Initialize the state storage if it isn't already.
582 if (!stateStorageInitialized_)
583 setStateSize();
584
585 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
586 "Belos::BlockGmresIter::initialize(): Cannot initialize state storage!");
587
588 const ScalarType zero = SCT::zero(), one = SCT::one();
589
590 // NOTE: In BlockGmresIter, V and Z are required!!!
591 // inconsitent multivectors widths and lengths will not be tolerated, and
592 // will be treated with exceptions.
593 //
594 std::string errstr("Belos::BlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
595
596 if (newstate.V != Teuchos::null && newstate.z != Teuchos::null) {
597
598 // initialize V_,z_, and curDim_
599
600 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
601 std::invalid_argument, errstr );
602 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
603 std::invalid_argument, errstr );
604 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*(numBlocks_+1),
605 std::invalid_argument, errstr );
606
607 curDim_ = newstate.curDim;
608 int lclDim = MVT::GetNumberVecs(*newstate.V);
609
610 // check size of Z
611 TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_, std::invalid_argument, errstr);
612
613
614 // copy basis vectors from newstate into V
615 if (newstate.V != V_) {
616 // only copy over the first block and print a warning.
617 if (curDim_ == 0 && lclDim > blockSize_) {
618 om_->stream(Warnings) << "Belos::BlockGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
619 << "The block size however is only " << blockSize_ << std::endl
620 << "The last " << lclDim - blockSize_ << " vectors will be discarded." << std::endl;
621 }
622 std::vector<int> nevind(curDim_+blockSize_);
623 for (int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
624 Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
625 Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_, nevind );
626 MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
627
628 // done with local pointers
629 lclV = Teuchos::null;
630 }
631
632 // put data into z_, make sure old information is not still hanging around.
633 if (newstate.z != z_) {
634 z_->putScalar();
635 Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+blockSize_,blockSize_);
636 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclZ;
637 lclZ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
638 lclZ->assign(newZ);
639
640 // done with local pointers
641 lclZ = Teuchos::null;
642 }
643
644 }
645 else {
646
647 TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
648 "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
649
650 TEUCHOS_TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
651 "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial norms z_0.");
652 }
653
654 // the solver is initialized
655 initialized_ = true;
656
657 if (om_->isVerbosity( Debug ) ) {
658 // Check almost everything here
659 CheckList chk;
660 chk.checkV = true;
661 chk.checkArn = true;
662 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
663 }
664
665 }
666
667
669 // Iterate until the status test informs us we should stop.
670 template <class ScalarType, class MV, class OP>
672 {
673 //
674 // Allocate/initialize data structures
675 //
676 if (initialized_ == false) {
677 initialize();
678 }
679
680 // Compute the current search dimension.
681 int searchDim = blockSize_*numBlocks_;
682
684 // iterate until the status test tells us to stop.
685 //
686 // also break if our basis is full
687 //
688 while (stest_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
689
690 iter_++;
691
692 // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
693 int lclDim = curDim_ + blockSize_;
694
695 // Get the current part of the basis.
696 std::vector<int> curind(blockSize_);
697 for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
698 Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
699
700 // Get a view of the previous vectors.
701 // This is used for orthogonalization and for computing V^H K H
702 for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
703 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
704
705 // Compute the next std::vector in the Krylov basis: Vnext = Op*Vprev
706 lp_->apply(*Vprev,*Vnext);
707 Vprev = Teuchos::null;
708
709 // Remove all previous Krylov basis vectors from Vnext
710 // Get a view of all the previous vectors
711 std::vector<int> prevind(lclDim);
712 for (int i=0; i<lclDim; i++) { prevind[i] = i; }
713 Vprev = MVT::CloneView(*V_,prevind);
714 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
715
716 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
717 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
718 subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
719 ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
720 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
721 AsubH.append( subH );
722
723 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
724 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
725 subH2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
726 ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
727 subH2->putScalar(); // Initialize subdiagonal to zero
728 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
729
730 // Copy over the coefficients if we are saving the upper Hessenberg matrix,
731 // just in case we run into an error.
732 if (keepHessenberg_) {
733 // Copy over the orthogonalization coefficients.
734 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
735 subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
736 ( Teuchos::View,*R_,lclDim,blockSize_,0,curDim_ ) );
737 subR->assign(*subH);
738
739 // Copy over the lower diagonal block of the Hessenberg matrix.
740 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
741 subR2 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
742 ( Teuchos::View,*R_,blockSize_,blockSize_,lclDim,curDim_ ) );
743 subR2->assign(*subH2);
744 }
745
746 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,GmresIterationOrthoFailure,
747 "Belos::BlockGmresIter::iterate(): couldn't generate basis of full rank.");
748 //
749 // V has been extended, and H has been extended.
750 //
751 // Update the QR factorization of the upper Hessenberg matrix
752 //
753 updateLSQR();
754 //
755 // Update basis dim and release all pointers.
756 //
757 Vnext = Teuchos::null;
758 curDim_ += blockSize_;
759 //
760 // When required, monitor some orthogonalities
761 if (om_->isVerbosity( Debug ) ) {
762 // Check almost everything here
763 CheckList chk;
764 chk.checkV = true;
765 chk.checkArn = true;
766 om_->print( Debug, accuracyCheck(chk, ": after local update") );
767 }
768 else if (om_->isVerbosity( OrthoDetails ) ) {
769 CheckList chk;
770 chk.checkV = true;
771 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
772 }
773
774 } // end while (statusTest == false)
775
776 }
777
778
779 template<class ScalarType, class MV, class OP>
781 {
782 int i, j, maxidx;
783 ScalarType sigma, mu, vscale, maxelem;
784 const ScalarType zero = SCT::zero();
785
786 // Get correct dimension based on input "dim"
787 // Remember that ortho failures result in an exit before updateLSQR() is called.
788 // Therefore, it is possible that dim == curDim_.
789 int curDim = curDim_;
790 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
791 curDim = dim;
792 }
793
794 Teuchos::BLAS<int, ScalarType> blas;
795 //
796 // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
797 // system to upper-triangular form.
798 //
799 if (blockSize_ == 1) {
800 //
801 // QR factorization of Least-Squares system with Givens rotations
802 //
803 for (i=0; i<curDim; i++) {
804 //
805 // Apply previous Givens rotations to new column of Hessenberg matrix
806 //
807 blas.ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
808 }
809 //
810 // Calculate new Givens rotation
811 //
812 blas.ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
813 (*R_)(curDim+1,curDim) = zero;
814 //
815 // Update RHS w/ new transformation
816 //
817 blas.ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
818 }
819 else {
820 //
821 // QR factorization of Least-Squares system with Householder reflectors
822 //
823 for (j=0; j<blockSize_; j++) {
824 //
825 // Apply previous Householder reflectors to new block of Hessenberg matrix
826 //
827 for (i=0; i<curDim+j; i++) {
828 sigma = blas.DOT( blockSize_, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
829 sigma += (*R_)(i,curDim+j);
830 sigma *= SCT::conjugate(beta[i]);
831 blas.AXPY(blockSize_, ScalarType(-sigma), &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
832 (*R_)(i,curDim+j) -= sigma;
833 }
834 //
835 // Compute new Householder reflector
836 //
837 maxidx = blas.IAMAX( blockSize_+1, &(*R_)(curDim+j,curDim+j), 1 );
838 maxelem = SCT::magnitude((*R_)(curDim+j+maxidx-1,curDim+j));
839 for (i=0; i<blockSize_+1; i++)
840 (*R_)(curDim+j+i,curDim+j) /= maxelem;
841 sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j), 1,
842 &(*R_)(curDim+j+1,curDim+j), 1 );
843 MagnitudeType sign_Rjj = -SCT::real((*R_)(curDim+j,curDim+j)) /
844 SCT::magnitude(SCT::real(((*R_)(curDim+j,curDim+j))));
845 if (sigma == zero) {
846 beta[curDim + j] = zero;
847 } else {
848 mu = SCT::squareroot(SCT::conjugate((*R_)(curDim+j,curDim+j))*(*R_)(curDim+j,curDim+j)+sigma);
849 vscale = (*R_)(curDim+j,curDim+j) - Teuchos::as<ScalarType>(sign_Rjj)*mu;
850 beta[curDim+j] = -Teuchos::as<ScalarType>(sign_Rjj) * vscale / mu;
851 (*R_)(curDim+j,curDim+j) = Teuchos::as<ScalarType>(sign_Rjj)*maxelem*mu;
852 for (i=0; i<blockSize_; i++)
853 (*R_)(curDim+j+1+i,curDim+j) /= vscale;
854 }
855 //
856 // Apply new Householder reflector to rhs
857 //
858 for (i=0; i<blockSize_; i++) {
859 sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j),
860 1, &(*z_)(curDim+j+1,i), 1);
861 sigma += (*z_)(curDim+j,i);
862 sigma *= SCT::conjugate(beta[curDim+j]);
863 blas.AXPY(blockSize_, ScalarType(-sigma), &(*R_)(curDim+j+1,curDim+j),
864 1, &(*z_)(curDim+j+1,i), 1);
865 (*z_)(curDim+j,i) -= sigma;
866 }
867 }
868 } // end if (blockSize_ == 1)
869
870 // If the least-squares problem is updated wrt "dim" then update the curDim_.
871 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
872 curDim_ = dim + blockSize_;
873 }
874
875 } // end updateLSQR()
876
878 // Check accuracy, orthogonality, and other debugging stuff
879 //
880 // bools specify which tests we want to run (instead of running more than we actually care about)
881 //
882 // checkV : V orthonormal
883 //
884 // checkArn: check the Arnoldi factorization
885 //
886 // NOTE: This method needs to check the current dimension of the subspace, since it is possible to
887 // call this method when curDim_ = 0 (after initialization).
888 template <class ScalarType, class MV, class OP>
889 std::string BlockGmresIter<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
890 {
891 std::stringstream os;
892 os.precision(2);
893 os.setf(std::ios::scientific, std::ios::floatfield);
894 MagnitudeType tmp;
895
896 os << " Debugging checks: iteration " << iter_ << where << std::endl;
897
898 // index vectors for V and F
899 std::vector<int> lclind(curDim_);
900 for (int i=0; i<curDim_; i++) lclind[i] = i;
901 std::vector<int> bsind(blockSize_);
902 for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
903
904 Teuchos::RCP<const MV> lclV,lclF;
905 Teuchos::RCP<MV> lclAV;
906 if (curDim_)
907 lclV = MVT::CloneView(*V_,lclind);
908 lclF = MVT::CloneView(*V_,bsind);
909
910 if (chk.checkV) {
911 if (curDim_) {
912 tmp = ortho_->orthonormError(*lclV);
913 os << " >> Error in V^H M V == I : " << tmp << std::endl;
914 }
915 tmp = ortho_->orthonormError(*lclF);
916 os << " >> Error in F^H M F == I : " << tmp << std::endl;
917 if (curDim_) {
918 tmp = ortho_->orthogError(*lclV,*lclF);
919 os << " >> Error in V^H M F == 0 : " << tmp << std::endl;
920 }
921 }
922
923 if (chk.checkArn) {
924
925 if (curDim_) {
926 // Compute AV
927 lclAV = MVT::Clone(*V_,curDim_);
928 lp_->apply(*lclV,*lclAV);
929
930 // Compute AV - VH
931 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
932 Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
933 MVT::MvTimesMatAddMv( -one, *lclV, subH, one, *lclAV );
934
935 // Compute FB_k^T - (AV-VH)
936 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
937 blockSize_,curDim_, curDim_ );
938 MVT::MvTimesMatAddMv( -one, *lclF, curB, one, *lclAV );
939
940 // Compute || FE_k^T - (AV-VH) ||
941 std::vector<MagnitudeType> arnNorms( curDim_ );
942 ortho_->norm( *lclAV, arnNorms );
943
944 for (int i=0; i<curDim_; i++) {
945 os << " >> Error in Krylov factorization (R = AV-VH-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
946 }
947 }
948 }
949
950 os << std::endl;
951
952 return os.str();
953 }
954
955} // end Belos namespace
956
957#endif /* BELOS_BLOCK_GMRES_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration.
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.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
Teuchos::SerialDenseVector< int, ScalarType > beta
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > R_
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SCT::magnitudeType MagnitudeType
Teuchos::SerialDenseVector< int, ScalarType > sn
void setStateSize()
Method for initalizing the state storage needed by block GMRES.
virtual ~BlockGmresIter()
Destructor.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this linear prob...
std::string accuracyCheck(const CheckList &chk, const std::string &where) const
Check accuracy of Arnoldi factorization.
Teuchos::SerialDenseVector< int, MagnitudeType > cs
const Teuchos::RCP< OutputManager< ScalarType > > om_
BlockGmresIter(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)
BlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver option...
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
MultiVecTraits< ScalarType, MV > MVT
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
bool isInitialized()
States whether the solver has been initialized or not.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > z_
void setBlockSize(int blockSize)
Set the blocksize.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
void initializeGmres(GmresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
OperatorTraits< ScalarType, MV, OP > OPT
int getNumIters() const
Get the current iteration count.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
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.
@ OrthoDetails
Structure to contain pointers to GmresIteration state variables.
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.