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

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11