Belos Version of the Day
Loading...
Searching...
No Matches
BelosPseudoBlockGmresIter.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_PSEUDO_BLOCK_GMRES_ITER_HPP
43#define BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51#include "BelosIteration.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
83
84
89 template <class ScalarType, class MV>
91
92 typedef Teuchos::ScalarTraits<ScalarType> SCT;
93 typedef typename SCT::magnitudeType MagnitudeType;
94
99 int curDim;
101 std::vector<Teuchos::RCP<const MV> > V;
107 std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > H;
109 std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > R;
111 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > Z;
113 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > sn;
114 std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > > cs;
115
117 H(0), R(0), Z(0),
118 sn(0), cs(0)
119 {}
120 };
121
123
124
132 PseudoBlockGmresIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
133 {}};
134
136
137
138 template<class ScalarType, class MV, class OP>
139 class PseudoBlockGmresIter : virtual public Iteration<ScalarType,MV,OP> {
140
141 public:
142
143 //
144 // Convenience typedefs
145 //
148 typedef Teuchos::ScalarTraits<ScalarType> SCT;
149 typedef typename SCT::magnitudeType MagnitudeType;
150
152
153
162 PseudoBlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
163 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
164 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
165 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
166 Teuchos::ParameterList &params );
167
171
172
174
175
197 void iterate();
198
221
226 {
228 initialize(empty);
229 }
230
240 state.curDim = curDim_;
241 state.V.resize(numRHS_);
242 state.H.resize(numRHS_);
243 state.Z.resize(numRHS_);
244 state.sn.resize(numRHS_);
245 state.cs.resize(numRHS_);
246 for (int i=0; i<numRHS_; ++i) {
247 state.V[i] = V_[i];
248 state.H[i] = H_[i];
249 state.Z[i] = Z_[i];
250 state.sn[i] = sn_[i];
251 state.cs[i] = cs_[i];
252 }
253 return state;
254 }
255
257
258
260
261
263 int getNumIters() const { return iter_; }
264
266 void resetNumIters( int iter = 0 ) { iter_ = iter; }
267
285 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
286
288
293 Teuchos::RCP<MV> getCurrentUpdate() const;
294
296
299 void updateLSQR( int dim = -1 );
300
302 int getCurSubspaceDim() const {
303 if (!initialized_) return 0;
304 return curDim_;
305 };
306
308 int getMaxSubspaceDim() const { return numBlocks_; }
309
311
312
314
315
317 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
318
320 int getBlockSize() const { return 1; }
321
323 void setBlockSize(int blockSize) {
324 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
325 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
326 }
327
329 int getNumBlocks() const { return numBlocks_; }
330
332 void setNumBlocks(int numBlocks);
333
335 bool isInitialized() { return initialized_; }
336
338
339 private:
340
341 //
342 // Classes inputed through constructor that define the linear problem to be solved.
343 //
344 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
345 const Teuchos::RCP<OutputManager<ScalarType> > om_;
346 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
347 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
348
349 //
350 // Algorithmic parameters
351 //
352 // numRHS_ is the current number of linear systems being solved.
353 int numRHS_;
354 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
355 int numBlocks_;
356
357 // Storage for QR factorization of the least squares system.
358 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
359 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
360
361 // Pointers to a work vector used to improve aggregate performance.
362 Teuchos::RCP<MV> U_vec_, AU_vec_;
363
364 // Pointers to the current right-hand side and solution multivecs being solved for.
365 Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
366
367 //
368 // Current solver state
369 //
370 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
371 // is capable of running; _initialize is controlled by the initialize() member method
372 // For the implications of the state of initialized_, please see documentation for initialize()
373 bool initialized_;
374
375 // Current subspace dimension, and number of iterations performed.
376 int curDim_, iter_;
377
378 //
379 // State Storage
380 //
381 std::vector<Teuchos::RCP<MV> > V_;
382 //
383 // Projected matrices
384 // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
385 //
386 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
387 //
388 // QR decomposition of Projected matrices for solving the least squares system HY = B.
389 // R_: Upper triangular reduction of H
390 // Z_: Q applied to right-hand side of the least squares system
391 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
392 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
393 };
394
396 // Constructor.
397 template<class ScalarType, class MV, class OP>
399 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
400 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
401 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
402 Teuchos::ParameterList &params ):
403 lp_(problem),
404 om_(printer),
405 stest_(tester),
406 ortho_(ortho),
407 numRHS_(0),
408 numBlocks_(0),
409 initialized_(false),
410 curDim_(0),
411 iter_(0)
412 {
413 // Get the maximum number of blocks allowed for each Krylov subspace
414 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
415 "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
416 int nb = Teuchos::getParameter<int>(params, "Num Blocks");
417
418 setNumBlocks( nb );
419 }
420
422 // Set the block size and make necessary adjustments.
423 template <class ScalarType, class MV, class OP>
425 {
426 // This routine only allocates space; it doesn't not perform any computation
427 // any change in size will invalidate the state of the solver.
428
429 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
430
431 numBlocks_ = numBlocks;
432 curDim_ = 0;
433
434 initialized_ = false;
435 }
436
438 // Get the current update from this subspace.
439 template <class ScalarType, class MV, class OP>
441 {
442 //
443 // If this is the first iteration of the Arnoldi factorization,
444 // there is no update, so return Teuchos::null.
445 //
446 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
447 if (curDim_==0) {
448 return currentUpdate;
449 } else {
450 currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
451 std::vector<int> index(1), index2(curDim_);
452 for (int i=0; i<curDim_; ++i) {
453 index2[i] = i;
454 }
455 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
456 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
457 Teuchos::BLAS<int,ScalarType> blas;
458
459 for (int i=0; i<numRHS_; ++i) {
460 index[0] = i;
461 Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
462 //
463 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
464 //
465 Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
466 //
467 // Solve the least squares problem and compute current solutions.
468 //
469 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
470 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
471 H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
472
473 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
474 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
475 }
476 }
477 return currentUpdate;
478 }
479
480
482 // Get the native residuals stored in this iteration.
483 // Note: No residual vector will be returned by Gmres.
484 template <class ScalarType, class MV, class OP>
485 Teuchos::RCP<const MV>
487 getNativeResiduals (std::vector<MagnitudeType> *norms) const
488 {
489 typedef typename Teuchos::ScalarTraits<ScalarType> STS;
490
491 if (norms)
492 { // Resize the incoming std::vector if necessary. The type
493 // cast avoids the compiler warning resulting from a signed /
494 // unsigned integer comparison.
495 if (static_cast<int> (norms->size()) < numRHS_)
496 norms->resize (numRHS_);
497
498 Teuchos::BLAS<int, ScalarType> blas;
499 for (int j = 0; j < numRHS_; ++j)
500 {
501 const ScalarType curNativeResid = (*Z_[j])(curDim_);
502 (*norms)[j] = STS::magnitude (curNativeResid);
503 }
504 }
505 return Teuchos::null;
506 }
507
508
509 template <class ScalarType, class MV, class OP>
510 void
513 {
514 using Teuchos::RCP;
515
516 // (Re)set the number of right-hand sides, by interrogating the
517 // current LinearProblem to solve.
518 this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
519
520 // NOTE: In PseudoBlockGmresIter, V and Z are required!!!
521 // Inconsistent multivectors widths and lengths will not be tolerated, and
522 // will be treated with exceptions.
523 //
524 std::string errstr ("Belos::PseudoBlockGmresIter::initialize(): "
525 "Specified multivectors must have a consistent "
526 "length and width.");
527
528 // Check that newstate has V and Z arrays with nonzero length.
529 TEUCHOS_TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0,
530 std::invalid_argument,
531 "Belos::PseudoBlockGmresIter::initialize(): "
532 "V and/or Z was not specified in the input state; "
533 "the V and/or Z arrays have length zero.");
534
535 // In order to create basis multivectors, we have to clone them
536 // from some already existing multivector. We require that at
537 // least one of the right-hand side B and left-hand side X in the
538 // LinearProblem be non-null. Thus, we can clone from either B or
539 // X. We prefer to close from B, since B is in the range of the
540 // operator A and the basis vectors should also be in the range of
541 // A (the first basis vector is a scaled residual vector).
542 // However, if B is not given, we will try our best by cloning
543 // from X.
544 RCP<const MV> lhsMV = lp_->getLHS();
545 RCP<const MV> rhsMV = lp_->getRHS();
546
547 // If the right-hand side is null, we make do with the left-hand
548 // side, otherwise we use the right-hand side.
549 RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
550 //RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
551
552 TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
553 std::invalid_argument,
554 "Belos::PseudoBlockGmresIter::initialize(): "
555 "The linear problem to solve does not specify multi"
556 "vectors from which we can clone basis vectors. The "
557 "right-hand side(s), left-hand side(s), or both should "
558 "be nonnull.");
559
560 // Check the new dimension is not more that the maximum number of
561 // allowable blocks.
562 TEUCHOS_TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1,
563 std::invalid_argument,
564 errstr);
565 curDim_ = newstate.curDim;
566
567 // Initialize the state storage. If the subspace has not be
568 // initialized before, generate it using the right-hand side or
569 // left-hand side from the LinearProblem lp_ to solve.
570 V_.resize(numRHS_);
571 for (int i=0; i<numRHS_; ++i) {
572 // Create a new vector if we need to. We "need to" if the
573 // current vector V_[i] is null, or if it doesn't have enough
574 // columns.
575 if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
576 V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
577 }
578 // Check that the newstate vector newstate.V[i] has dimensions
579 // consistent with those of V_[i].
580 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V[i]) != MVT::GetGlobalLength(*V_[i]),
581 std::invalid_argument, errstr );
582 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
583 std::invalid_argument, errstr );
584 //
585 // If newstate.V[i] and V_[i] are not identically the same
586 // vector, then copy newstate.V[i] into V_[i].
587 //
588 int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
589 if (newstate.V[i] != V_[i]) {
590 // Only copy over the first block and print a warning.
591 if (curDim_ == 0 && lclDim > 1) {
592 om_->stream(Warnings)
593 << "Belos::PseudoBlockGmresIter::initialize(): the solver was "
594 << "initialized with a kernel of " << lclDim
595 << std::endl
596 << "The block size however is only " << 1
597 << std::endl
598 << "The last " << lclDim - 1 << " vectors will be discarded."
599 << std::endl;
600 }
601 std::vector<int> nevind (curDim_ + 1);
602 for (int j = 0; j < curDim_ + 1; ++j)
603 nevind[j] = j;
604
605 RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
606 RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
607 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
608 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
609 MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
610
611 // Done with local pointers
612 lclV = Teuchos::null;
613 }
614 }
615
616
617 // Check size of Z
618 Z_.resize(numRHS_);
619 for (int i=0; i<numRHS_; ++i) {
620 // Create a vector if we need to.
621 if (Z_[i] == Teuchos::null) {
622 Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
623 }
624 if (Z_[i]->length() < numBlocks_+1) {
625 Z_[i]->shapeUninitialized(numBlocks_+1, 1);
626 }
627
628 // Check that the newstate vector is consistent.
629 TEUCHOS_TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
630
631 // Put data into Z_, make sure old information is not still hanging around.
632 if (newstate.Z[i] != Z_[i]) {
633 if (curDim_==0)
634 Z_[i]->putScalar();
635
636 Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
637 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
638 lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
639 lclZ->assign(newZ);
640
641 // Done with local pointers
642 lclZ = Teuchos::null;
643 }
644 }
645
646
647 // Check size of H
648 H_.resize(numRHS_);
649 for (int i=0; i<numRHS_; ++i) {
650 // Create a matrix if we need to.
651 if (H_[i] == Teuchos::null) {
652 H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
653 }
654 if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
655 H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
656 }
657
658 // Put data into H_ if it exists, make sure old information is not still hanging around.
659 if ((int)newstate.H.size() == numRHS_) {
660
661 // Check that the newstate matrix is consistent.
662 TEUCHOS_TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument,
663 "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
664
665 if (newstate.H[i] != H_[i]) {
666 // H_[i]->putScalar();
667
668 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
669 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
670 lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
671 lclH->assign(newH);
672
673 // Done with local pointers
674 lclH = Teuchos::null;
675 }
676 }
677 }
678
680 // Reinitialize storage for least squares solve
681 //
682 cs_.resize(numRHS_);
683 sn_.resize(numRHS_);
684
685 // Copy over rotation angles if they exist
686 if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
687 for (int i=0; i<numRHS_; ++i) {
688 if (cs_[i] != newstate.cs[i])
689 cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
690 if (sn_[i] != newstate.sn[i])
691 sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
692 }
693 }
694
695 // Resize or create the vectors as necessary
696 for (int i=0; i<numRHS_; ++i) {
697 if (cs_[i] == Teuchos::null)
698 cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
699 else
700 cs_[i]->resize(numBlocks_+1);
701 if (sn_[i] == Teuchos::null)
702 sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
703 else
704 sn_[i]->resize(numBlocks_+1);
705 }
706
707 // the solver is initialized
708 initialized_ = true;
709
710 }
711
712
714 // Iterate until the status test informs us we should stop.
715 template <class ScalarType, class MV, class OP>
717 {
718 //
719 // Allocate/initialize data structures
720 //
721 if (initialized_ == false) {
722 initialize();
723 }
724
725 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
726 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
727
728 // Compute the current search dimension.
729 int searchDim = numBlocks_;
730 //
731 // Associate each initial block of V_[i] with U_vec[i]
732 // Reset the index vector (this might have been changed if there was a restart)
733 //
734 std::vector<int> index(1);
735 std::vector<int> index2(1);
736 index[0] = curDim_;
737 Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
738
739 // Create AU_vec to hold A*U_vec.
740 Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
741
742 for (int i=0; i<numRHS_; ++i) {
743 index2[0] = i;
744 Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
745 Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
746 MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
747 }
748
750 // iterate until the status test tells us to stop.
751 //
752 // also break if our basis is full
753 //
754 while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
755
756 iter_++;
757 //
758 // Apply the operator to _work_vector
759 //
760 lp_->apply( *U_vec, *AU_vec );
761 //
762 //
763 // Resize index.
764 //
765 int num_prev = curDim_+1;
766 index.resize( num_prev );
767 for (int i=0; i<num_prev; ++i) {
768 index[i] = i;
769 }
770 //
771 // Orthogonalize next Krylov vector for each right-hand side.
772 //
773 for (int i=0; i<numRHS_; ++i) {
774 //
775 // Get previous Krylov vectors.
776 //
777 Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
778 Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
779 //
780 // Get a view of the new candidate std::vector.
781 //
782 index2[0] = i;
783 Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
784 //
785 // Get a view of the current part of the upper-hessenberg matrix.
786 //
787 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new
788 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
789 ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
790 Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
791
792 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
793 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
794 ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
795 //
796 // Orthonormalize the new block of the Krylov expansion
797 // NOTE: Rank deficiencies are not checked because this is a single-std::vector Krylov method.
798 //
799 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
800 //
801 // NOTE: V_new is a copy of the iter+1 vector in V_[i], so the normalized vector has to be
802 // be copied back in when V_new is changed.
803 //
804 index2[0] = curDim_+1;
805 Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
806 MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
807 }
808 //
809 // Now _AU_vec is the new _U_vec, so swap these two vectors.
810 // NOTE: This alleviates the need for allocating a vector for AU_vec each iteration.
811 //
812 Teuchos::RCP<MV> tmp_AU_vec = U_vec;
813 U_vec = AU_vec;
814 AU_vec = tmp_AU_vec;
815 //
816 // V has been extended, and H has been extended.
817 //
818 // Update the QR factorization of the upper Hessenberg matrix
819 //
820 updateLSQR();
821 //
822 // Update basis dim and release all pointers.
823 //
824 curDim_ += 1;
825 //
826 } // end while (statusTest == false)
827
828 }
829
831 // Update the least squares solution for each right-hand side.
832 template<class ScalarType, class MV, class OP>
834 {
835 // Get correct dimension based on input "dim"
836 // Remember that ortho failures result in an exit before updateLSQR() is called.
837 // Therefore, it is possible that dim == curDim_.
838 int curDim = curDim_;
839 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
840 curDim = dim;
841 }
842
843 int i, j;
844 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
845
846 Teuchos::BLAS<int, ScalarType> blas;
847
848 for (i=0; i<numRHS_; ++i) {
849 //
850 // Update the least-squares QR for each linear system.
851 //
852 // QR factorization of Least-Squares system with Givens rotations
853 //
854 for (j=0; j<curDim; j++) {
855 //
856 // Apply previous Givens rotations to new column of Hessenberg matrix
857 //
858 blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
859 }
860 //
861 // Calculate new Givens rotation
862 //
863 blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
864 (*H_[i])(curDim+1,curDim) = zero;
865 //
866 // Update RHS w/ new transformation
867 //
868 blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
869 }
870
871 } // end updateLSQR()
872
873} // end Belos namespace
874
875#endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which describes the basic interface to the linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
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.
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...
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
PseudoBlockGmresIterOrthoFailure(const std::string &what_arg)
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
MultiVecTraits< ScalarType, MV > MVT
bool isInitialized()
States whether the solver has been initialized or not.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
OperatorTraits< ScalarType, MV, OP > OPT
PseudoBlockGmresIter(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)
PseudoBlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver ...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the "native" residual vectors.
void setBlockSize(int blockSize)
Set the blocksize.
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.
void resetNumIters(int iter=0)
Reset the iteration count.
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
int getNumIters() const
Get the current iteration count.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > R
The current upper-triangular matrix from the QR reduction of H.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
Teuchos::ScalarTraits< ScalarType > SCT
int curDim
The current dimension of the reduction.

Generated for Belos by doxygen 1.10.0