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 
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  //
150 
152 
153 
163  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
166  Teuchos::ParameterList &params );
167 
169  virtual ~PseudoBlockGmresIter() {};
171 
172 
174 
175 
197  void iterate();
198 
221 
225  void initialize()
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 
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  //
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,
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();
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  //
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>
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 
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) {
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);
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) {
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_);
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  //
789  ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
791 
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 
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 */
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.
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 ...

Generated on Wed Mar 27 2024 09:25:07 for Belos by doxygen 1.8.5