Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
51 #include "BelosGmresIteration.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 
82 template<class ScalarType, class MV, class OP>
83 class BlockFGmresIter : virtual public GmresIteration<ScalarType,MV,OP> {
84 
85  public:
86 
87  //
88  // Convenience typedefs
89  //
94 
96 
97 
108  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
111  Teuchos::ParameterList &params );
112 
114  virtual ~BlockFGmresIter() {};
116 
117 
119 
120 
142  void iterate();
143 
166 
170  void initialize()
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 
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  //
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.
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  //
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
328 };
329 
331  // Constructor.
332  template<class ScalarType, class MV, class OP>
335  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
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
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;
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
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;
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;
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>
489  {
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 ();
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.
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>
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) {
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;
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 
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 
579  MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
580  std::invalid_argument, errstr );
582  MVT::GetNumberVecs(*newstate.V) < blockSize_,
583  std::invalid_argument, errstr );
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
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 {
631  newstate.V == Teuchos::null,std::invalid_argument,
632  "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial kernel V_0.");
633 
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;
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);
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  {
737 
738  const ScalarType zero = STS::zero ();
739  const ScalarType two = (STS::one () + STS::one());
740  ScalarType sigma, mu, vscale, maxelem;
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 */
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
void iterate()
This method performs block FGmres iterations until the status test indicates the need to stop or an e...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
void setBlockSize(int blockSize)
Set the blocksize.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Class which manages the output and verbosity of the Belos solvers.
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...
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::ScalarTraits< ScalarType > SCT
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
Teuchos::RCP< const MV > V
The current Krylov basis.
T & get(const std::string &name, T def_value)
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Declaration of basic traits for the multivector type.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
OrdinalType IAMAX(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Structure to contain pointers to GmresIteration state variables.
A pure virtual class for defining the status tests for the Belos iterative solvers.
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 resetNumIters(int iter=0)
Reset the iteration count.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
Traits class which defines basic operations on multivectors.
bool isParameter(const std::string &name) const
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
This class implements the block flexible GMRES iteration, where a block Krylov subspace is constructe...
ScalarType DOT(const OrdinalType &n, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy) const
virtual ~BlockFGmresIter()
Destructor.
void initializeGmres(GmresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
bool isInitialized()
States whether the solver has been initialized or not.
MultiVecTraits< ScalarType, MV > MVT
OrdinalType numCols() const
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Class which defines basic traits for the operator type.
SCT::magnitudeType MagnitudeType
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.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
OrdinalType numRows() const
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
bool is_null() const

Generated on Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5