Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
53 #include "BelosLinearProblem.hpp"
54 #include "BelosMatOrthoManager.hpp"
55 #include "BelosOutputManager.hpp"
56 #include "BelosStatusTest.hpp"
57 #include "BelosOperatorTraits.hpp"
58 #include "BelosMultiVecTraits.hpp"
59 
60 #include "Teuchos_BLAS.hpp"
63 #include "Teuchos_ScalarTraits.hpp"
65 #include "Teuchos_TimeMonitor.hpp"
66 
80 namespace Belos {
81 
83 
84 
89  template <class ScalarType, class MV>
91 
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 
138  PseudoBlockGmresIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
139  {}};
140 
148  PseudoBlockGmresIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
149  {}};
150 
152 
153 
154  template<class ScalarType, class MV, class OP>
155  class PseudoBlockGmresIter : virtual public Iteration<ScalarType,MV,OP> {
156 
157  public:
158 
159  //
160  // Convenience typedefs
161  //
166 
168 
169 
179  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
182  Teuchos::ParameterList &params );
183 
185  virtual ~PseudoBlockGmresIter() {};
187 
188 
190 
191 
213  void iterate();
214 
237 
241  void initialize()
242  {
244  initialize(empty);
245  }
246 
256  state.curDim = curDim_;
257  state.V.resize(numRHS_);
258  state.H.resize(numRHS_);
259  state.Z.resize(numRHS_);
260  state.sn.resize(numRHS_);
261  state.cs.resize(numRHS_);
262  for (int i=0; i<numRHS_; ++i) {
263  state.V[i] = V_[i];
264  state.H[i] = H_[i];
265  state.Z[i] = Z_[i];
266  state.sn[i] = sn_[i];
267  state.cs[i] = cs_[i];
268  }
269  return state;
270  }
271 
273 
274 
276 
277 
279  int getNumIters() const { return iter_; }
280 
282  void resetNumIters( int iter = 0 ) { iter_ = iter; }
283 
301  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
302 
304 
310 
312 
315  void updateLSQR( int dim = -1 );
316 
318  int getCurSubspaceDim() const {
319  if (!initialized_) return 0;
320  return curDim_;
321  };
322 
324  int getMaxSubspaceDim() const { return numBlocks_; }
325 
327 
328 
330 
331 
333  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
334 
336  int getBlockSize() const { return 1; }
337 
339  void setBlockSize(int blockSize) {
340  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
341  "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
342  }
343 
345  int getNumBlocks() const { return numBlocks_; }
346 
348  void setNumBlocks(int numBlocks);
349 
351  bool isInitialized() { return initialized_; }
352 
354 
355  private:
356 
357  //
358  // Classes inputed through constructor that define the linear problem to be solved.
359  //
364 
365  //
366  // Algorithmic parameters
367  //
368  // numRHS_ is the current number of linear systems being solved.
369  int numRHS_;
370  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
371  int numBlocks_;
372 
373  // Storage for QR factorization of the least squares system.
374  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
375  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
376 
377  // Pointers to a work vector used to improve aggregate performance.
378  Teuchos::RCP<MV> U_vec_, AU_vec_;
379 
380  // Pointers to the current right-hand side and solution multivecs being solved for.
381  Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
382 
383  //
384  // Current solver state
385  //
386  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
387  // is capable of running; _initialize is controlled by the initialize() member method
388  // For the implications of the state of initialized_, please see documentation for initialize()
389  bool initialized_;
390 
391  // Current subspace dimension, and number of iterations performed.
392  int curDim_, iter_;
393 
394  //
395  // State Storage
396  //
397  std::vector<Teuchos::RCP<MV> > V_;
398  //
399  // Projected matrices
400  // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
401  //
402  std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
403  //
404  // QR decomposition of Projected matrices for solving the least squares system HY = B.
405  // R_: Upper triangular reduction of H
406  // Z_: Q applied to right-hand side of the least squares system
407  std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
408  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
409  };
410 
412  // Constructor.
413  template<class ScalarType, class MV, class OP>
415  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
418  Teuchos::ParameterList &params ):
419  lp_(problem),
420  om_(printer),
421  stest_(tester),
422  ortho_(ortho),
423  numRHS_(0),
424  numBlocks_(0),
425  initialized_(false),
426  curDim_(0),
427  iter_(0)
428  {
429  // Get the maximum number of blocks allowed for each Krylov subspace
430  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
431  "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
432  int nb = Teuchos::getParameter<int>(params, "Num Blocks");
433 
434  setNumBlocks( nb );
435  }
436 
438  // Set the block size and make necessary adjustments.
439  template <class ScalarType, class MV, class OP>
441  {
442  // This routine only allocates space; it doesn't not perform any computation
443  // any change in size will invalidate the state of the solver.
444 
445  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
446 
447  numBlocks_ = numBlocks;
448  curDim_ = 0;
449 
450  initialized_ = false;
451  }
452 
454  // Get the current update from this subspace.
455  template <class ScalarType, class MV, class OP>
457  {
458  //
459  // If this is the first iteration of the Arnoldi factorization,
460  // there is no update, so return Teuchos::null.
461  //
462  Teuchos::RCP<MV> currentUpdate = Teuchos::null;
463  if (curDim_==0) {
464  return currentUpdate;
465  } else {
466  currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
467  std::vector<int> index(1), index2(curDim_);
468  for (int i=0; i<curDim_; ++i) {
469  index2[i] = i;
470  }
471  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
472  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
474 
475  for (int i=0; i<numRHS_; ++i) {
476  index[0] = i;
477  Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
478  //
479  // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
480  //
481  Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
482  //
483  // Solve the least squares problem and compute current solutions.
484  //
486  Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
487  H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
488 
489  Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
490  MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
491  }
492  }
493  return currentUpdate;
494  }
495 
496 
498  // Get the native residuals stored in this iteration.
499  // Note: No residual vector will be returned by Gmres.
500  template <class ScalarType, class MV, class OP>
503  getNativeResiduals (std::vector<MagnitudeType> *norms) const
504  {
505  typedef typename Teuchos::ScalarTraits<ScalarType> STS;
506 
507  if (norms)
508  { // Resize the incoming std::vector if necessary. The type
509  // cast avoids the compiler warning resulting from a signed /
510  // unsigned integer comparison.
511  if (static_cast<int> (norms->size()) < numRHS_)
512  norms->resize (numRHS_);
513 
515  for (int j = 0; j < numRHS_; ++j)
516  {
517  const ScalarType curNativeResid = (*Z_[j])(curDim_);
518  (*norms)[j] = STS::magnitude (curNativeResid);
519  }
520  }
521  return Teuchos::null;
522  }
523 
524 
525  template <class ScalarType, class MV, class OP>
526  void
529  {
530  using Teuchos::RCP;
531 
532  // (Re)set the number of right-hand sides, by interrogating the
533  // current LinearProblem to solve.
534  this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
535 
536  // NOTE: In PseudoBlockGmresIter, V and Z are required!!!
537  // Inconsistent multivectors widths and lengths will not be tolerated, and
538  // will be treated with exceptions.
539  //
540  std::string errstr ("Belos::PseudoBlockGmresIter::initialize(): "
541  "Specified multivectors must have a consistent "
542  "length and width.");
543 
544  // Check that newstate has V and Z arrays with nonzero length.
545  TEUCHOS_TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0,
546  std::invalid_argument,
547  "Belos::PseudoBlockGmresIter::initialize(): "
548  "V and/or Z was not specified in the input state; "
549  "the V and/or Z arrays have length zero.");
550 
551  // In order to create basis multivectors, we have to clone them
552  // from some already existing multivector. We require that at
553  // least one of the right-hand side B and left-hand side X in the
554  // LinearProblem be non-null. Thus, we can clone from either B or
555  // X. We prefer to close from B, since B is in the range of the
556  // operator A and the basis vectors should also be in the range of
557  // A (the first basis vector is a scaled residual vector).
558  // However, if B is not given, we will try our best by cloning
559  // from X.
560  RCP<const MV> lhsMV = lp_->getLHS();
561  RCP<const MV> rhsMV = lp_->getRHS();
562 
563  // If the right-hand side is null, we make do with the left-hand
564  // side, otherwise we use the right-hand side.
565  RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
566  //RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
567 
568  TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
569  std::invalid_argument,
570  "Belos::PseudoBlockGmresIter::initialize(): "
571  "The linear problem to solve does not specify multi"
572  "vectors from which we can clone basis vectors. The "
573  "right-hand side(s), left-hand side(s), or both should "
574  "be nonnull.");
575 
576  // Check the new dimension is not more that the maximum number of
577  // allowable blocks.
578  TEUCHOS_TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1,
579  std::invalid_argument,
580  errstr);
581  curDim_ = newstate.curDim;
582 
583  // Initialize the state storage. If the subspace has not be
584  // initialized before, generate it using the right-hand side or
585  // left-hand side from the LinearProblem lp_ to solve.
586  V_.resize(numRHS_);
587  for (int i=0; i<numRHS_; ++i) {
588  // Create a new vector if we need to. We "need to" if the
589  // current vector V_[i] is null, or if it doesn't have enough
590  // columns.
591  if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
592  V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
593  }
594  // Check that the newstate vector newstate.V[i] has dimensions
595  // consistent with those of V_[i].
596  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V[i]) != MVT::GetGlobalLength(*V_[i]),
597  std::invalid_argument, errstr );
598  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
599  std::invalid_argument, errstr );
600  //
601  // If newstate.V[i] and V_[i] are not identically the same
602  // vector, then copy newstate.V[i] into V_[i].
603  //
604  int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
605  if (newstate.V[i] != V_[i]) {
606  // Only copy over the first block and print a warning.
607  if (curDim_ == 0 && lclDim > 1) {
608  om_->stream(Warnings)
609  << "Belos::PseudoBlockGmresIter::initialize(): the solver was "
610  << "initialized with a kernel of " << lclDim
611  << std::endl
612  << "The block size however is only " << 1
613  << std::endl
614  << "The last " << lclDim - 1 << " vectors will be discarded."
615  << std::endl;
616  }
617  std::vector<int> nevind (curDim_ + 1);
618  for (int j = 0; j < curDim_ + 1; ++j)
619  nevind[j] = j;
620 
621  RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
622  RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
623  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
624  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
625  MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
626 
627  // Done with local pointers
628  lclV = Teuchos::null;
629  }
630  }
631 
632 
633  // Check size of Z
634  Z_.resize(numRHS_);
635  for (int i=0; i<numRHS_; ++i) {
636  // Create a vector if we need to.
637  if (Z_[i] == Teuchos::null) {
639  }
640  if (Z_[i]->length() < numBlocks_+1) {
641  Z_[i]->shapeUninitialized(numBlocks_+1, 1);
642  }
643 
644  // Check that the newstate vector is consistent.
645  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
646 
647  // Put data into Z_, make sure old information is not still hanging around.
648  if (newstate.Z[i] != Z_[i]) {
649  if (curDim_==0)
650  Z_[i]->putScalar();
651 
652  Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
654  lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
655  lclZ->assign(newZ);
656 
657  // Done with local pointers
658  lclZ = Teuchos::null;
659  }
660  }
661 
662 
663  // Check size of H
664  H_.resize(numRHS_);
665  for (int i=0; i<numRHS_; ++i) {
666  // Create a matrix if we need to.
667  if (H_[i] == Teuchos::null) {
669  }
670  if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
671  H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
672  }
673 
674  // Put data into H_ if it exists, make sure old information is not still hanging around.
675  if ((int)newstate.H.size() == numRHS_) {
676 
677  // Check that the newstate matrix is consistent.
678  TEUCHOS_TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument,
679  "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
680 
681  if (newstate.H[i] != H_[i]) {
682  // H_[i]->putScalar();
683 
684  Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
686  lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
687  lclH->assign(newH);
688 
689  // Done with local pointers
690  lclH = Teuchos::null;
691  }
692  }
693  }
694 
696  // Reinitialize storage for least squares solve
697  //
698  cs_.resize(numRHS_);
699  sn_.resize(numRHS_);
700 
701  // Copy over rotation angles if they exist
702  if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
703  for (int i=0; i<numRHS_; ++i) {
704  if (cs_[i] != newstate.cs[i])
705  cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
706  if (sn_[i] != newstate.sn[i])
707  sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
708  }
709  }
710 
711  // Resize or create the vectors as necessary
712  for (int i=0; i<numRHS_; ++i) {
713  if (cs_[i] == Teuchos::null)
714  cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
715  else
716  cs_[i]->resize(numBlocks_+1);
717  if (sn_[i] == Teuchos::null)
718  sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
719  else
720  sn_[i]->resize(numBlocks_+1);
721  }
722 
723  // the solver is initialized
724  initialized_ = true;
725 
726  /*
727  if (om_->isVerbosity( Debug ) ) {
728  // Check almost everything here
729  CheckList chk;
730  chk.checkV = true;
731  chk.checkArn = true;
732  chk.checkAux = true;
733  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
734  }
735  */
736 
737  }
738 
739 
741  // Iterate until the status test informs us we should stop.
742  template <class ScalarType, class MV, class OP>
744  {
745  //
746  // Allocate/initialize data structures
747  //
748  if (initialized_ == false) {
749  initialize();
750  }
751 
752  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
753  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
754 
755  // Compute the current search dimension.
756  int searchDim = numBlocks_;
757  //
758  // Associate each initial block of V_[i] with U_vec[i]
759  // Reset the index vector (this might have been changed if there was a restart)
760  //
761  std::vector<int> index(1);
762  std::vector<int> index2(1);
763  index[0] = curDim_;
764  Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
765 
766  // Create AU_vec to hold A*U_vec.
767  Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
768 
769  for (int i=0; i<numRHS_; ++i) {
770  index2[0] = i;
771  Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
772  Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
773  MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
774  }
775 
777  // iterate until the status test tells us to stop.
778  //
779  // also break if our basis is full
780  //
781  while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
782 
783  iter_++;
784  //
785  // Apply the operator to _work_vector
786  //
787  lp_->apply( *U_vec, *AU_vec );
788  //
789  //
790  // Resize index.
791  //
792  int num_prev = curDim_+1;
793  index.resize( num_prev );
794  for (int i=0; i<num_prev; ++i) {
795  index[i] = i;
796  }
797  //
798  // Orthogonalize next Krylov vector for each right-hand side.
799  //
800  for (int i=0; i<numRHS_; ++i) {
801  //
802  // Get previous Krylov vectors.
803  //
804  Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
805  Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
806  //
807  // Get a view of the new candidate std::vector.
808  //
809  index2[0] = i;
810  Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
811  //
812  // Get a view of the current part of the upper-hessenberg matrix.
813  //
816  ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
818 
821  ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
822  //
823  // Orthonormalize the new block of the Krylov expansion
824  // NOTE: Rank deficiencies are not checked because this is a single-std::vector Krylov method.
825  //
826  ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
827  //
828  // NOTE: V_new is a copy of the iter+1 vector in V_[i], so the normalized vector has to be
829  // be copied back in when V_new is changed.
830  //
831  index2[0] = curDim_+1;
832  Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
833  MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
834  }
835  //
836  // Now _AU_vec is the new _U_vec, so swap these two vectors.
837  // NOTE: This alleviates the need for allocating a vector for AU_vec each iteration.
838  //
839  Teuchos::RCP<MV> tmp_AU_vec = U_vec;
840  U_vec = AU_vec;
841  AU_vec = tmp_AU_vec;
842  //
843  // V has been extended, and H has been extended.
844  //
845  // Update the QR factorization of the upper Hessenberg matrix
846  //
847  updateLSQR();
848  //
849  // Update basis dim and release all pointers.
850  //
851  curDim_ += 1;
852  //
853  /*
854  // When required, monitor some orthogonalities
855  if (om_->isVerbosity( Debug ) ) {
856  // Check almost everything here
857  CheckList chk;
858  chk.checkV = true;
859  chk.checkArn = true;
860  om_->print( Debug, accuracyCheck(chk, ": after local update") );
861  }
862  else if (om_->isVerbosity( OrthoDetails ) ) {
863  CheckList chk;
864  chk.checkV = true;
865  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
866  }
867  */
868 
869  } // end while (statusTest == false)
870 
871  }
872 
874  // Update the least squares solution for each right-hand side.
875  template<class ScalarType, class MV, class OP>
877  {
878  // Get correct dimension based on input "dim"
879  // Remember that ortho failures result in an exit before updateLSQR() is called.
880  // Therefore, it is possible that dim == curDim_.
881  int curDim = curDim_;
882  if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
883  curDim = dim;
884  }
885 
886  int i, j;
887  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
888 
890 
891  for (i=0; i<numRHS_; ++i) {
892  //
893  // Update the least-squares QR for each linear system.
894  //
895  // QR factorization of Least-Squares system with Givens rotations
896  //
897  for (j=0; j<curDim; j++) {
898  //
899  // Apply previous Givens rotations to new column of Hessenberg matrix
900  //
901  blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
902  }
903  //
904  // Calculate new Givens rotation
905  //
906  blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
907  (*H_[i])(curDim+1,curDim) = zero;
908  //
909  // Update RHS w/ new transformation
910  //
911  blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
912  }
913 
914  } // end updateLSQR()
915 
916 } // end Belos namespace
917 
918 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */
ScalarType * values() const
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
Collection of types and exceptions used within the Belos solvers.
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > R
The current upper-triangular matrix from the QR reduction of H.
Teuchos::ScalarTraits< ScalarType > SCT
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
int resize(OrdinalType length_in)
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the &quot;native&quot; residual vectors.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
Structure to contain pointers to PseudoBlockGmresIter state variables.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
bool isInitialized()
States whether the solver has been initialized or not.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
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 ...
PseudoBlockGmresIterOrthoFailure(const std::string &what_arg)
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
PseudoBlockGmresIterInitFailure(const std::string &what_arg)
Teuchos::ScalarTraits< ScalarType > SCT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void resetNumIters(int iter=0)
Reset the iteration count.
MultiVecTraits< ScalarType, MV > MVT
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
OperatorTraits< ScalarType, MV, OP > OPT
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
Class which defines basic traits for the operator type.
int curDim
The current dimension of the reduction.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given&#39;s rotation coefficients.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos header file which uses auto-configuration information to include necessary C++ headers...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
OrdinalType stride() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
PseudoBlockGmresIterInitFailure is thrown when the PseudoBlockGmresIter object is unable to generate ...

Generated on Fri Jun 5 2020 10:20:45 for Belos by doxygen 1.8.5