Belos Version of the Day
Loading...
Searching...
No Matches
BelosBlockFGmresIter.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_FGMRES_ITER_HPP
43#define BELOS_BLOCK_FGMRES_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_SerialDenseMatrix.hpp"
62#include "Teuchos_SerialDenseVector.hpp"
63#include "Teuchos_ScalarTraits.hpp"
64#include "Teuchos_ParameterList.hpp"
65#include "Teuchos_TimeMonitor.hpp"
66
80namespace Belos {
81
82template<class ScalarType, class MV, class OP>
83class BlockFGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
84
85 public:
86
87 //
88 // Convenience typedefs
89 //
92 typedef Teuchos::ScalarTraits<ScalarType> SCT;
93 typedef typename SCT::magnitudeType MagnitudeType;
94
96
97
107 BlockFGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
108 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
109 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
110 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
111 Teuchos::ParameterList &params );
112
114 virtual ~BlockFGmresIter() {};
116
117
119
120
142 void iterate();
143
166
171 {
173 initializeGmres(empty);
174 }
175
185 state.curDim = curDim_;
186 state.V = V_;
187 state.Z = Z_;
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
231 int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }
232
234
235
237
238
240 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
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 methods
271 //
273 void setStateSize();
274
275 //
276 // Classes inputed through constructor that define the linear problem to be solved.
277 //
278 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
279 const Teuchos::RCP<OutputManager<ScalarType> > om_;
280 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
281 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
282
283 //
284 // Algorithmic parameters
285 //
286 // blockSize_ is the solver block size.
287 // It controls the number of vectors added to the basis on each iteration.
288 int blockSize_;
289 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
290 int numBlocks_;
291
292 // Storage for QR factorization of the least squares system.
293 Teuchos::SerialDenseVector<int,ScalarType> beta, sn;
294 Teuchos::SerialDenseVector<int,MagnitudeType> cs;
295
296 //
297 // Current solver state
298 //
299 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
300 // is capable of running; _initialize is controlled by the initialize() member method
301 // For the implications of the state of initialized_, please see documentation for initialize()
302 bool initialized_;
303
304 // stateStorageInitialized_ specified that the state storage has be initialized to the current
305 // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
306 // generated without the right-hand side or solution vectors.
307 bool stateStorageInitialized_;
308
309 // Current subspace dimension, and number of iterations performed.
310 int curDim_, iter_;
311
312 //
313 // State Storage
314 //
315 Teuchos::RCP<MV> V_;
316 Teuchos::RCP<MV> Z_;
317 //
318 // Projected matrices
319 // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
320 //
321 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
322 //
323 // QR decomposition of Projected matrices for solving the least squares system HY = B.
324 // R_: Upper triangular reduction of H
325 // z_: Q applied to right-hand side of the least squares system
326 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
327 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_;
328};
329
331 // Constructor.
332 template<class ScalarType, class MV, class OP>
334 BlockFGmresIter (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
335 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
336 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
337 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
338 Teuchos::ParameterList &params ):
339 lp_(problem),
340 om_(printer),
341 stest_(tester),
342 ortho_(ortho),
343 blockSize_(0),
344 numBlocks_(0),
345 initialized_(false),
346 stateStorageInitialized_(false),
347 curDim_(0),
348 iter_(0)
349 {
350 // Get the maximum number of blocks allowed for this Krylov subspace
351 TEUCHOS_TEST_FOR_EXCEPTION(
352 ! params.isParameter ("Num Blocks"), std::invalid_argument,
353 "Belos::BlockFGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
354 const int nb = params.get<int> ("Num Blocks");
355
356 // Set the block size and allocate data.
357 const int bs = params.get ("Block Size", 1);
358 setSize (bs, nb);
359 }
360
362 // Set the block size and make necessary adjustments.
363 template <class ScalarType, class MV, class OP>
364 void BlockFGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
365 {
366 // This routine only allocates space; it doesn't not perform any computation
367 // any change in size will invalidate the state of the solver.
368
369 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockFGmresIter::setSize was passed a non-positive argument.");
370 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
371 // do nothing
372 return;
373 }
374
375 if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
376 stateStorageInitialized_ = false;
377
378 blockSize_ = blockSize;
379 numBlocks_ = numBlocks;
380
381 initialized_ = false;
382 curDim_ = 0;
383
384 // Use the current blockSize_ and numBlocks_ to initialize the state storage.
385 setStateSize();
386
387 }
388
390 // Setup the state storage.
391 template <class ScalarType, class MV, class OP>
393 {
394 using Teuchos::RCP;
395 using Teuchos::rcp;
396 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
397
398 if (! stateStorageInitialized_) {
399 // Check if there is any multivector to clone from.
400 RCP<const MV> lhsMV = lp_->getLHS();
401 RCP<const MV> rhsMV = lp_->getRHS();
402 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
403 stateStorageInitialized_ = false;
404 return;
405 }
406 else {
408 // blockSize*numBlocks dependent
409 //
410 int newsd = blockSize_*(numBlocks_+1);
411
412 if (blockSize_==1) {
413 cs.resize (newsd);
414 sn.resize (newsd);
415 }
416 else {
417 beta.resize (newsd);
418 }
419
420 // Initialize the state storage
421 TEUCHOS_TEST_FOR_EXCEPTION(
422 blockSize_ * static_cast<ptrdiff_t> (numBlocks_) > MVT::GetGlobalLength (*rhsMV),
423 std::invalid_argument, "Belos::BlockFGmresIter::setStateSize(): "
424 "Cannot generate a Krylov basis with dimension larger the operator!");
425
426 // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
427 if (V_ == Teuchos::null) {
428 // Get the multivector that is not null.
429 RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
430 TEUCHOS_TEST_FOR_EXCEPTION(
431 tmp == Teuchos::null, std::invalid_argument,
432 "Belos::BlockFGmresIter::setStateSize(): "
433 "linear problem does not specify multivectors to clone from.");
434 V_ = MVT::Clone (*tmp, newsd);
435 }
436 else {
437 // Generate V_ by cloning itself ONLY if more space is needed.
438 if (MVT::GetNumberVecs (*V_) < newsd) {
439 RCP<const MV> tmp = V_;
440 V_ = MVT::Clone (*tmp, newsd);
441 }
442 }
443
444 if (Z_ == Teuchos::null) {
445 // Get the multivector that is not null.
446 RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
447 TEUCHOS_TEST_FOR_EXCEPTION(
448 tmp == Teuchos::null, std::invalid_argument,
449 "Belos::BlockFGmresIter::setStateSize(): "
450 "linear problem does not specify multivectors to clone from.");
451 Z_ = MVT::Clone (*tmp, newsd);
452 }
453 else {
454 // Generate Z_ by cloning itself ONLY if more space is needed.
455 if (MVT::GetNumberVecs (*Z_) < newsd) {
456 RCP<const MV> tmp = Z_;
457 Z_ = MVT::Clone (*tmp, newsd);
458 }
459 }
460
461 // Generate H_ only if it doesn't exist, otherwise resize it.
462 if (H_ == Teuchos::null) {
463 H_ = rcp (new SDM (newsd, newsd-blockSize_));
464 }
465 else {
466 H_->shapeUninitialized (newsd, newsd - blockSize_);
467 }
468
469 // TODO: Insert logic so that Hessenberg matrix can be saved and reduced matrix is stored in R_
470 //R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( newsd, newsd-blockSize_ ) );
471 // Generate z_ only if it doesn't exist, otherwise resize it.
472 if (z_ == Teuchos::null) {
473 z_ = rcp (new SDM (newsd, blockSize_));
474 }
475 else {
476 z_->shapeUninitialized (newsd, blockSize_);
477 }
478
479 // State storage has now been initialized.
480 stateStorageInitialized_ = true;
481 }
482 }
483 }
484
485
486 template <class ScalarType, class MV, class OP>
487 Teuchos::RCP<MV>
489 {
490 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
491
492 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
493 if (curDim_ == 0) {
494 // If this is the first iteration of the Arnoldi factorization,
495 // then there is no update, so return Teuchos::null.
496 return currentUpdate;
497 }
498 else {
499 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero ();
500 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one ();
501 Teuchos::BLAS<int,ScalarType> blas;
502
503 currentUpdate = MVT::Clone (*Z_, blockSize_);
504
505 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
506 SDM y (Teuchos::Copy, *z_, curDim_, blockSize_);
507
508 // Solve the least squares problem.
509 blas.TRSM (Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
510 Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
511 H_->values (), H_->stride (), y.values (), y.stride ());
512
513 // Compute the current update.
514 std::vector<int> index (curDim_);
515 for (int i = 0; i < curDim_; ++i) {
516 index[i] = i;
517 }
518 Teuchos::RCP<const MV> Zjp1 = MVT::CloneView (*Z_, index);
519 MVT::MvTimesMatAddMv (one, *Zjp1, y, zero, *currentUpdate);
520 }
521 return currentUpdate;
522 }
523
524
525 template <class ScalarType, class MV, class OP>
526 Teuchos::RCP<const MV>
528 getNativeResiduals (std::vector<MagnitudeType> *norms) const
529 {
530 // NOTE: Make sure the incoming std::vector is the correct size!
531 if (norms != NULL && (int)norms->size() < blockSize_) {
532 norms->resize (blockSize_);
533 }
534
535 if (norms != NULL) {
536 Teuchos::BLAS<int, ScalarType> blas;
537 for (int j = 0; j < blockSize_; ++j) {
538 (*norms)[j] = blas.NRM2 (blockSize_, &(*z_)(curDim_, j), 1);
539 }
540 }
541
542 // FGmres does not return a residual (multi)vector.
543 return Teuchos::null;
544 }
545
546
547 template <class ScalarType, class MV, class OP>
550 {
551 using Teuchos::RCP;
552 using Teuchos::rcp;
553 using std::endl;
554 typedef Teuchos::ScalarTraits<ScalarType> STS;
555 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
556 const ScalarType ZERO = STS::zero ();
557 const ScalarType ONE = STS::one ();
558
559 // Initialize the state storage if it isn't already.
560 if (! stateStorageInitialized_) {
561 setStateSize ();
562 }
563
564 TEUCHOS_TEST_FOR_EXCEPTION(
565 ! stateStorageInitialized_, std::invalid_argument,
566 "Belos::BlockFGmresIter::initialize(): Cannot initialize state storage!");
567
568 // NOTE: In BlockFGmresIter, V and Z are required!!! Inconsistent
569 // multivectors widths and lengths will not be tolerated, and will
570 // be treated with exceptions.
571 const char errstr[] = "Belos::BlockFGmresIter::initialize(): The given "
572 "multivectors must have a consistent length and width.";
573
574 if (! newstate.V.is_null () && ! newstate.z.is_null ()) {
575
576 // initialize V_,z_, and curDim_
577
578 TEUCHOS_TEST_FOR_EXCEPTION(
579 MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
580 std::invalid_argument, errstr );
581 TEUCHOS_TEST_FOR_EXCEPTION(
582 MVT::GetNumberVecs(*newstate.V) < blockSize_,
583 std::invalid_argument, errstr );
584 TEUCHOS_TEST_FOR_EXCEPTION(
585 newstate.curDim > blockSize_*(numBlocks_+1),
586 std::invalid_argument, errstr );
587
588 curDim_ = newstate.curDim;
589 const int lclDim = MVT::GetNumberVecs(*newstate.V);
590
591 // check size of Z
592 TEUCHOS_TEST_FOR_EXCEPTION(
593 newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_,
594 std::invalid_argument, errstr);
595
596 // copy basis vectors from newstate into V
597 if (newstate.V != V_) {
598 // only copy over the first block and print a warning.
599 if (curDim_ == 0 && lclDim > blockSize_) {
600 std::ostream& warn = om_->stream (Warnings);
601 warn << "Belos::BlockFGmresIter::initialize(): the solver was "
602 << "initialized with a kernel of " << lclDim << endl
603 << "The block size however is only " << blockSize_ << endl
604 << "The last " << lclDim - blockSize_
605 << " vectors will be discarded." << endl;
606 }
607 std::vector<int> nevind (curDim_ + blockSize_);
608 for (int i = 0; i < curDim_ + blockSize_; ++i) {
609 nevind[i] = i;
610 }
611 RCP<const MV> newV = MVT::CloneView (*newstate.V, nevind);
612 RCP<MV> lclV = MVT::CloneViewNonConst (*V_, nevind);
613 MVT::MvAddMv (ONE, *newV, ZERO, *newV, *lclV);
614
615 // done with local pointers
616 lclV = Teuchos::null;
617 }
618
619 // put data into z_, make sure old information is not still hanging around.
620 if (newstate.z != z_) {
621 z_->putScalar();
622 SDM newZ (Teuchos::View, *newstate.z, curDim_ + blockSize_, blockSize_);
623 RCP<SDM> lclz;
624 lclz = rcp (new SDM (Teuchos::View, *z_, curDim_ + blockSize_, blockSize_));
625 lclz->assign (newZ);
626 lclz = Teuchos::null; // done with local pointers
627 }
628 }
629 else {
630 TEUCHOS_TEST_FOR_EXCEPTION(
631 newstate.V == Teuchos::null,std::invalid_argument,
632 "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial kernel V_0.");
633
634 TEUCHOS_TEST_FOR_EXCEPTION(
635 newstate.z == Teuchos::null,std::invalid_argument,
636 "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial norms z_0.");
637 }
638
639 // the solver is initialized
640 initialized_ = true;
641 }
642
643
644 template <class ScalarType, class MV, class OP>
646 {
647 using Teuchos::Array;
648 using Teuchos::null;
649 using Teuchos::RCP;
650 using Teuchos::rcp;
651 using Teuchos::View;
652 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
653
654 // Allocate/initialize data structures
655 if (initialized_ == false) {
656 initialize();
657 }
658
659 // Compute the current search dimension.
660 const int searchDim = blockSize_ * numBlocks_;
661
662 // Iterate until the status test tells us to stop.
663 // Raise an exception if a computed block is not full rank.
664 while (stest_->checkStatus (this) != Passed && curDim_+blockSize_ <= searchDim) {
665 ++iter_;
666
667 // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
668 const int lclDim = curDim_ + blockSize_;
669
670 // Get the current part of the basis.
671 std::vector<int> curind (blockSize_);
672 for (int i = 0; i < blockSize_; ++i) {
673 curind[i] = lclDim + i;
674 }
675 RCP<MV> Vnext = MVT::CloneViewNonConst (*V_, curind);
676
677 // Get a view of the previous vectors.
678 // This is used for orthogonalization and for computing V^H K H.
679 for (int i = 0; i < blockSize_; ++i) {
680 curind[i] = curDim_ + i;
681 }
682 RCP<const MV> Vprev = MVT::CloneView (*V_, curind);
683 RCP<MV> Znext = MVT::CloneViewNonConst (*Z_, curind);
684
685 // Compute the next (multi)vector in the Krylov basis: Znext = M*Vprev
686 lp_->applyRightPrec (*Vprev, *Znext);
687 Vprev = null;
688
689 // Compute the next (multi)vector in the Krylov basis: Vnext = A*Znext
690 lp_->applyOp (*Znext, *Vnext);
691 Znext = null;
692
693 // Remove all previous Krylov basis vectors from Vnext
694 // Get a view of all the previous vectors
695 std::vector<int> prevind (lclDim);
696 for (int i = 0; i < lclDim; ++i) {
697 prevind[i] = i;
698 }
699 Vprev = MVT::CloneView (*V_, prevind);
700 Array<RCP<const MV> > AVprev (1, Vprev);
701
702 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
703 RCP<SDM> subH = rcp (new SDM (View, *H_, lclDim, blockSize_, 0, curDim_));
704 Array<RCP<SDM> > AsubH;
705 AsubH.append (subH);
706
707 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
708 RCP<SDM> subR = rcp (new SDM (View, *H_, blockSize_, blockSize_, lclDim, curDim_));
709 const int rank = ortho_->projectAndNormalize (*Vnext, AsubH, subR, AVprev);
710 TEUCHOS_TEST_FOR_EXCEPTION(
711 rank != blockSize_, GmresIterationOrthoFailure,
712 "Belos::BlockFGmresIter::iterate(): After orthogonalization, the new "
713 "basis block does not have full rank. It contains " << blockSize_
714 << " vector" << (blockSize_ != 1 ? "s" : "")
715 << ", but its rank is " << rank << ".");
716
717 //
718 // V has been extended, and H has been extended.
719 //
720 // Update the QR factorization of the upper Hessenberg matrix
721 //
722 updateLSQR ();
723 //
724 // Update basis dim and release all pointers.
725 //
726 Vnext = null;
727 curDim_ += blockSize_;
728 } // end while (statusTest == false)
729 }
730
731
732 template<class ScalarType, class MV, class OP>
734 {
735 typedef Teuchos::ScalarTraits<ScalarType> STS;
736 typedef Teuchos::ScalarTraits<MagnitudeType> STM;
737
738 const ScalarType zero = STS::zero ();
739 const ScalarType two = (STS::one () + STS::one());
740 ScalarType sigma, mu, vscale, maxelem;
741 Teuchos::BLAS<int, ScalarType> blas;
742
743 // Get correct dimension based on input 'dim'. Remember that
744 // orthogonalization failures result in an exit before
745 // updateLSQR() is called. Therefore, it is possible that dim ==
746 // curDim_.
747 int curDim = curDim_;
748 if (dim >= curDim_ && dim < getMaxSubspaceDim ()) {
749 curDim = dim;
750 }
751
752 // Apply previous transformations, and compute new transformation
753 // to reduce upper Hessenberg system to upper triangular form.
754 // The type of transformation we use depends the block size. We
755 // use Givens rotations for a block size of 1, and Householder
756 // reflectors otherwise.
757 if (blockSize_ == 1) {
758 // QR factorization of upper Hessenberg matrix using Givens rotations
759 for (int i = 0; i < curDim; ++i) {
760 // Apply previous Givens rotations to new column of Hessenberg matrix
761 blas.ROT (1, &(*H_)(i, curDim), 1, &(*H_)(i+1, curDim), 1, &cs[i], &sn[i]);
762 }
763 // Calculate new Givens rotation
764 blas.ROTG (&(*H_)(curDim, curDim), &(*H_)(curDim+1, curDim), &cs[curDim], &sn[curDim]);
765 (*H_)(curDim+1, curDim) = zero;
766
767 // Update RHS w/ new transformation
768 blas.ROT (1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim]);
769 }
770 else {
771 // QR factorization of least-squares system using Householder reflectors.
772 for (int j = 0; j < blockSize_; ++j) {
773 // Apply previous Householder reflectors to new block of Hessenberg matrix
774 for (int i = 0; i < curDim + j; ++i) {
775 sigma = blas.DOT (blockSize_, &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
776 sigma += (*H_)(i,curDim+j);
777 sigma *= beta[i];
778 blas.AXPY (blockSize_, ScalarType(-sigma), &(*H_)(i+1,i), 1, &(*H_)(i+1,curDim+j), 1);
779 (*H_)(i,curDim+j) -= sigma;
780 }
781
782 // Compute new Householder reflector
783 const int maxidx = blas.IAMAX (blockSize_+1, &(*H_)(curDim+j,curDim+j), 1);
784 maxelem = (*H_)(curDim + j + maxidx - 1, curDim + j);
785 for (int i = 0; i < blockSize_ + 1; ++i) {
786 (*H_)(curDim+j+i,curDim+j) /= maxelem;
787 }
788 sigma = blas.DOT (blockSize_, &(*H_)(curDim + j + 1, curDim + j), 1,
789 &(*H_)(curDim + j + 1, curDim + j), 1);
790 if (sigma == zero) {
791 beta[curDim + j] = zero;
792 } else {
793 mu = STS::squareroot ((*H_)(curDim+j,curDim+j)*(*H_)(curDim+j,curDim+j)+sigma);
794 if (STS::real ((*H_)(curDim + j, curDim + j)) < STM::zero ()) {
795 vscale = (*H_)(curDim+j,curDim+j) - mu;
796 } else {
797 vscale = -sigma / ((*H_)(curDim+j, curDim+j) + mu);
798 }
799 beta[curDim+j] = two * vscale * vscale / (sigma + vscale*vscale);
800 (*H_)(curDim+j, curDim+j) = maxelem*mu;
801 for (int i = 0; i < blockSize_; ++i) {
802 (*H_)(curDim+j+1+i,curDim+j) /= vscale;
803 }
804 }
805
806 // Apply new Householder reflector to the right-hand side.
807 for (int i = 0; i < blockSize_; ++i) {
808 sigma = blas.DOT (blockSize_, &(*H_)(curDim+j+1,curDim+j),
809 1, &(*z_)(curDim+j+1,i), 1);
810 sigma += (*z_)(curDim+j,i);
811 sigma *= beta[curDim+j];
812 blas.AXPY (blockSize_, ScalarType(-sigma), &(*H_)(curDim+j+1,curDim+j),
813 1, &(*z_)(curDim+j+1,i), 1);
814 (*z_)(curDim+j,i) -= sigma;
815 }
816 }
817 } // end if (blockSize_ == 1)
818
819 // If the least-squares problem is updated wrt "dim" then update curDim_.
820 if (dim >= curDim_ && dim < getMaxSubspaceDim ()) {
821 curDim_ = dim + blockSize_;
822 }
823 } // end updateLSQR()
824
825} // namespace Belos
826
827#endif /* BELOS_BLOCK_FGMRES_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 flexible GMRES iteration, where a block Krylov subspace is constructe...
int getNumIters() const
Get the current iteration count.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void iterate()
This method performs block FGmres iterations until the status test indicates the need to stop or an e...
MultiVecTraits< ScalarType, MV > MVT
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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.
SCT::magnitudeType MagnitudeType
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.
void resetNumIters(int iter=0)
Reset the iteration count.
virtual ~BlockFGmresIter()
Destructor.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
void setBlockSize(int blockSize)
Set the blocksize.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
BlockFGmresIter(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)
BlockFGmresIter constructor with linear problem, solver utilities, and parameter list of solver optio...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
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...
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Teuchos::ScalarTraits< ScalarType > SCT
bool isInitialized()
States whether the solver has been initialized or not.
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.
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.
Teuchos::RCP< const MV > Z
The current preconditioned Krylov basis (only used in flexible GMRES).
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.

Generated for Belos by doxygen 1.10.0