Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBlockGmresIter.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_GMRES_ITER_HPP
43 #define BELOS_BLOCK_GMRES_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 BlockGmresIter : 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 ~BlockGmresIter() {};
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.H = H_;
188  state.R = R_;
189  state.z = z_;
190  return state;
191  }
192 
194 
195 
197 
198 
200  int getNumIters() const { return iter_; }
201 
203  void resetNumIters( int iter = 0 ) { iter_ = iter; }
204 
207  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
208 
210 
216 
218 
221  void updateLSQR( int dim = -1 );
222 
224  int getCurSubspaceDim() const {
225  if (!initialized_) return 0;
226  return curDim_;
227  };
228 
230  int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }
231 
233 
234 
236 
237 
239  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
240 
242  int getBlockSize() const { return blockSize_; }
243 
245  void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }
246 
248  int getNumBlocks() const { return numBlocks_; }
249 
251  void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }
252 
259  void setSize(int blockSize, int numBlocks);
260 
262  bool isInitialized() { return initialized_; }
263 
265 
266  private:
267 
268  //
269  // Internal methods
270  //
272  void setStateSize();
273 
274  //
275  // Classes inputed through constructor that define the linear problem to be solved.
276  //
281 
282  //
283  // Algorithmic parameters
284  //
285  // blockSize_ is the solver block size.
286  // It controls the number of vectors added to the basis on each iteration.
287  int blockSize_;
288  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
289  int numBlocks_;
290 
291  // Storage for QR factorization of the least squares system.
294 
295  //
296  // Current solver state
297  //
298  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
299  // is capable of running; _initialize is controlled by the initialize() member method
300  // For the implications of the state of initialized_, please see documentation for initialize()
301  bool initialized_;
302 
303  // stateStorageInitialized_ specifies that the state storage has be initialized to the current
304  // blockSize_ and numBlocks_. This initialization may be postponed if the linear problem was
305  // generated without the right-hand side or solution vectors.
306  bool stateStorageInitialized_;
307 
308  // keepHessenberg_ specifies that the iteration must keep the Hessenberg matrix formed via the
309  // Arnoldi factorization and the upper triangular matrix that is the Hessenberg matrix reduced via
310  // QR factorization separate.
311  bool keepHessenberg_;
312 
313  // initHessenberg_ specifies that the iteration should reinitialize the Hessenberg matrix by zeroing
314  // out all entries before an iteration is started.
315  bool initHessenberg_;
316 
317  // Current subspace dimension, and number of iterations performed.
318  int curDim_, iter_;
319 
320  //
321  // State Storage
322  //
323  Teuchos::RCP<MV> V_;
324  //
325  // Projected matrices
326  // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
327  //
329  //
330  // QR decomposition of Projected matrices for solving the least squares system HY = B.
331  // R_: Upper triangular reduction of H
332  // z_: Q applied to right-hand side of the least squares system
335 };
336 
338  // Constructor.
339  template<class ScalarType, class MV, class OP>
341  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
344  Teuchos::ParameterList &params ):
345  lp_(problem),
346  om_(printer),
347  stest_(tester),
348  ortho_(ortho),
349  blockSize_(0),
350  numBlocks_(0),
351  initialized_(false),
352  stateStorageInitialized_(false),
353  keepHessenberg_(false),
354  initHessenberg_(false),
355  curDim_(0),
356  iter_(0)
357  {
358  // Find out whether we are saving the Hessenberg matrix.
359  keepHessenberg_ = params.get("Keep Hessenberg", false);
360 
361  // Find out whether we are initializing the Hessenberg matrix.
362  initHessenberg_ = params.get("Initialize Hessenberg", false);
363 
364  // Get the maximum number of blocks allowed for this Krylov subspace
365  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
366  "Belos::BlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
367  int nb = Teuchos::getParameter<int>(params, "Num Blocks");
368 
369  // Set the block size and allocate data
370  int bs = params.get("Block Size", 1);
371  setSize( bs, nb );
372  }
373 
375  // Set the block size and make necessary adjustments.
376  template <class ScalarType, class MV, class OP>
377  void BlockGmresIter<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
378  {
379  // This routine only allocates space; it doesn't not perform any computation
380  // any change in size will invalidate the state of the solver.
381 
382  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setSize was passed a non-positive argument.");
383  if (blockSize == blockSize_ && numBlocks == numBlocks_) {
384  // do nothing
385  return;
386  }
387 
388  if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
389  stateStorageInitialized_ = false;
390 
391  blockSize_ = blockSize;
392  numBlocks_ = numBlocks;
393 
394  initialized_ = false;
395  curDim_ = 0;
396 
397  // Use the current blockSize_ and numBlocks_ to initialize the state storage.
398  setStateSize();
399 
400  }
401 
403  // Setup the state storage.
404  template <class ScalarType, class MV, class OP>
406  {
407  if (!stateStorageInitialized_) {
408 
409  // Check if there is any multivector to clone from.
410  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
411  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
412  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
413  stateStorageInitialized_ = false;
414  return;
415  }
416  else {
417 
419  // blockSize*numBlocks dependent
420  //
421  int newsd = blockSize_*(numBlocks_+1);
422 
423  if (blockSize_==1) {
424  cs.resize( newsd );
425  sn.resize( newsd );
426  }
427  else {
428  beta.resize( newsd );
429  }
430 
431  // Initialize the state storage
432  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_*static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
433  "Belos::BlockGmresIter::setStateSize(): Cannot generate a Krylov basis with dimension larger the operator!");
434 
435  // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
436  if (V_ == Teuchos::null) {
437  // Get the multivector that is not null.
438  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
439  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
440  "Belos::BlockGmresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
441  V_ = MVT::Clone( *tmp, newsd );
442  }
443  else {
444  // Generate V_ by cloning itself ONLY if more space is needed.
445  if (MVT::GetNumberVecs(*V_) < newsd) {
446  Teuchos::RCP<const MV> tmp = V_;
447  V_ = MVT::Clone( *tmp, newsd );
448  }
449  }
450 
451  // Generate R_ only if it doesn't exist, otherwise resize it.
452  if (R_ == Teuchos::null) {
454  }
455  if (initHessenberg_) {
456  R_->shape( newsd, newsd-blockSize_ );
457  }
458  else {
459  if (R_->numRows() < newsd || R_->numCols() < newsd-blockSize_) {
460  R_->shapeUninitialized( newsd, newsd-blockSize_ );
461  }
462  }
463 
464  // Generate H_ only if it doesn't exist, and we are keeping the upper Hessenberg matrix.
465  if (keepHessenberg_) {
466  if (H_ == Teuchos::null) {
468  }
469  if (initHessenberg_) {
470  H_->shape( newsd, newsd-blockSize_ );
471  }
472  else {
473  if (H_->numRows() < newsd || H_->numCols() < newsd-blockSize_) {
474  H_->shapeUninitialized( newsd, newsd-blockSize_ );
475  }
476  }
477  }
478  else {
479  // Point H_ and R_ at the same object.
480  H_ = R_;
481  }
482 
483  // Generate z_ only if it doesn't exist, otherwise resize it.
484  if (z_ == Teuchos::null) {
486  }
487  if (z_-> numRows() < newsd || z_->numCols() < blockSize_) {
488  z_->shapeUninitialized( newsd, blockSize_ );
489  }
490 
491  // State storage has now been initialized.
492  stateStorageInitialized_ = true;
493  }
494  }
495  }
496 
498  // Get the current update from this subspace.
499  template <class ScalarType, class MV, class OP>
501  {
502  //
503  // If this is the first iteration of the Arnoldi factorization,
504  // there is no update, so return Teuchos::null.
505  //
506  Teuchos::RCP<MV> currentUpdate = Teuchos::null;
507  if (curDim_==0) {
508  return currentUpdate;
509  } else {
510  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
511  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
513  currentUpdate = MVT::Clone( *V_, blockSize_ );
514  //
515  // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
516  //
517  Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, *z_, curDim_, blockSize_ );
518  //
519  // Solve the least squares problem.
520  //
522  Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
523  R_->values(), R_->stride(), y.values(), y.stride() );
524  //
525  // Compute the current update.
526  //
527  std::vector<int> index(curDim_);
528  for ( int i=0; i<curDim_; i++ ) {
529  index[i] = i;
530  }
531  Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
532  MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
533  }
534  return currentUpdate;
535  }
536 
537 
539  // Get the native residuals stored in this iteration.
540  // Note: No residual std::vector will be returned by Gmres.
541  template <class ScalarType, class MV, class OP>
543  {
544  //
545  // NOTE: Make sure the incoming std::vector is the correct size!
546  //
547  if ( norms && (int)norms->size() < blockSize_ )
548  norms->resize( blockSize_ );
549 
550  if (norms) {
552  for (int j=0; j<blockSize_; j++) {
553  (*norms)[j] = blas.NRM2( blockSize_, &(*z_)(curDim_, j), 1);
554  }
555  }
556  return Teuchos::null;
557  }
558 
559 
560 
562  // Initialize this iteration object
563  template <class ScalarType, class MV, class OP>
565  {
566  // Initialize the state storage if it isn't already.
567  if (!stateStorageInitialized_)
568  setStateSize();
569 
570  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
571  "Belos::BlockGmresIter::initialize(): Cannot initialize state storage!");
572 
574 
575  // NOTE: In BlockGmresIter, V and Z are required!!!
576  // inconsitent multivectors widths and lengths will not be tolerated, and
577  // will be treated with exceptions.
578  //
579  std::string errstr("Belos::BlockGmresIter::initialize(): Specified multivectors must have a consistent length and width.");
580 
581  if (newstate.V != Teuchos::null && newstate.z != Teuchos::null) {
582 
583  // initialize V_,z_, and curDim_
584 
585  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
586  std::invalid_argument, errstr );
587  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
588  std::invalid_argument, errstr );
589  TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*(numBlocks_+1),
590  std::invalid_argument, errstr );
591 
592  curDim_ = newstate.curDim;
593  int lclDim = MVT::GetNumberVecs(*newstate.V);
594 
595  // check size of Z
596  TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < blockSize_, std::invalid_argument, errstr);
597 
598 
599  // copy basis vectors from newstate into V
600  if (newstate.V != V_) {
601  // only copy over the first block and print a warning.
602  if (curDim_ == 0 && lclDim > blockSize_) {
603  om_->stream(Warnings) << "Belos::BlockGmresIter::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
604  << "The block size however is only " << blockSize_ << std::endl
605  << "The last " << lclDim - blockSize_ << " vectors will be discarded." << std::endl;
606  }
607  std::vector<int> nevind(curDim_+blockSize_);
608  for (int i=0; i<curDim_+blockSize_; i++) nevind[i] = i;
609  Teuchos::RCP<const MV> newV = MVT::CloneView( *newstate.V, nevind );
610  Teuchos::RCP<MV> lclV = MVT::CloneViewNonConst( *V_, nevind );
611  MVT::MvAddMv( one, *newV, zero, *newV, *lclV );
612 
613  // done with local pointers
614  lclV = Teuchos::null;
615  }
616 
617  // put data into z_, make sure old information is not still hanging around.
618  if (newstate.z != z_) {
619  z_->putScalar();
620  Teuchos::SerialDenseMatrix<int,ScalarType> newZ(Teuchos::View,*newstate.z,curDim_+blockSize_,blockSize_);
622  lclZ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*z_,curDim_+blockSize_,blockSize_) );
623  lclZ->assign(newZ);
624 
625  // done with local pointers
626  lclZ = Teuchos::null;
627  }
628 
629  }
630  else {
631 
632  TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,
633  "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial kernel V_0.");
634 
635  TEUCHOS_TEST_FOR_EXCEPTION(newstate.z == Teuchos::null,std::invalid_argument,
636  "Belos::BlockGmresIter::initialize(): BlockGmresStateIterState does not have initial norms z_0.");
637  }
638 
639  // the solver is initialized
640  initialized_ = true;
641 
642  /*
643  if (om_->isVerbosity( Debug ) ) {
644  // Check almost everything here
645  CheckList chk;
646  chk.checkV = true;
647  chk.checkArn = true;
648  chk.checkAux = true;
649  om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
650  }
651  */
652 
653  }
654 
655 
657  // Iterate until the status test informs us we should stop.
658  template <class ScalarType, class MV, class OP>
660  {
661  //
662  // Allocate/initialize data structures
663  //
664  if (initialized_ == false) {
665  initialize();
666  }
667 
668  // Compute the current search dimension.
669  int searchDim = blockSize_*numBlocks_;
670 
672  // iterate until the status test tells us to stop.
673  //
674  // also break if our basis is full
675  //
676  while (stest_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
677 
678  iter_++;
679 
680  // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
681  int lclDim = curDim_ + blockSize_;
682 
683  // Get the current part of the basis.
684  std::vector<int> curind(blockSize_);
685  for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
686  Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
687 
688  // Get a view of the previous vectors.
689  // This is used for orthogonalization and for computing V^H K H
690  for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
691  Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
692 
693  // Compute the next std::vector in the Krylov basis: Vnext = Op*Vprev
694  lp_->apply(*Vprev,*Vnext);
695  Vprev = Teuchos::null;
696 
697  // Remove all previous Krylov basis vectors from Vnext
698  // Get a view of all the previous vectors
699  std::vector<int> prevind(lclDim);
700  for (int i=0; i<lclDim; i++) { prevind[i] = i; }
701  Vprev = MVT::CloneView(*V_,prevind);
702  Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
703 
704  // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
707  ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
709  AsubH.append( subH );
710 
711  // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
714  ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
715  subH2->putScalar(); // Initialize subdiagonal to zero
716  int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subH2,AVprev);
717 
718  // Copy over the coefficients if we are saving the upper Hessenberg matrix,
719  // just in case we run into an error.
720  if (keepHessenberg_) {
721  // Copy over the orthogonalization coefficients.
724  ( Teuchos::View,*R_,lclDim,blockSize_,0,curDim_ ) );
725  subR->assign(*subH);
726 
727  // Copy over the lower diagonal block of the Hessenberg matrix.
730  ( Teuchos::View,*R_,blockSize_,blockSize_,lclDim,curDim_ ) );
731  subR2->assign(*subH2);
732  }
733 
735  "Belos::BlockGmresIter::iterate(): couldn't generate basis of full rank.");
736  //
737  // V has been extended, and H has been extended.
738  //
739  // Update the QR factorization of the upper Hessenberg matrix
740  //
741  updateLSQR();
742  //
743  // Update basis dim and release all pointers.
744  //
745  Vnext = Teuchos::null;
746  curDim_ += blockSize_;
747  //
748  /*
749  // When required, monitor some orthogonalities
750  if (om_->isVerbosity( Debug ) ) {
751  // Check almost everything here
752  CheckList chk;
753  chk.checkV = true;
754  chk.checkArn = true;
755  om_->print( Debug, accuracyCheck(chk, ": after local update") );
756  }
757  else if (om_->isVerbosity( OrthoDetails ) ) {
758  CheckList chk;
759  chk.checkV = true;
760  om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
761  }
762  */
763 
764  } // end while (statusTest == false)
765 
766  }
767 
768 
769  template<class ScalarType, class MV, class OP>
771  {
772  int i, j, maxidx;
773  ScalarType sigma, mu, vscale, maxelem;
775 
776  // Get correct dimension based on input "dim"
777  // Remember that ortho failures result in an exit before updateLSQR() is called.
778  // Therefore, it is possible that dim == curDim_.
779  int curDim = curDim_;
780  if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
781  curDim = dim;
782  }
783 
785  //
786  // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
787  // system to upper-triangular form.
788  //
789  if (blockSize_ == 1) {
790  //
791  // QR factorization of Least-Squares system with Givens rotations
792  //
793  for (i=0; i<curDim; i++) {
794  //
795  // Apply previous Givens rotations to new column of Hessenberg matrix
796  //
797  blas.ROT( 1, &(*R_)(i,curDim), 1, &(*R_)(i+1, curDim), 1, &cs[i], &sn[i] );
798  }
799  //
800  // Calculate new Givens rotation
801  //
802  blas.ROTG( &(*R_)(curDim,curDim), &(*R_)(curDim+1,curDim), &cs[curDim], &sn[curDim] );
803  (*R_)(curDim+1,curDim) = zero;
804  //
805  // Update RHS w/ new transformation
806  //
807  blas.ROT( 1, &(*z_)(curDim,0), 1, &(*z_)(curDim+1,0), 1, &cs[curDim], &sn[curDim] );
808  }
809  else {
810  //
811  // QR factorization of Least-Squares system with Householder reflectors
812  //
813  for (j=0; j<blockSize_; j++) {
814  //
815  // Apply previous Householder reflectors to new block of Hessenberg matrix
816  //
817  for (i=0; i<curDim+j; i++) {
818  sigma = blas.DOT( blockSize_, &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
819  sigma += (*R_)(i,curDim+j);
820  sigma *= beta[i];
821  blas.AXPY(blockSize_, ScalarType(-sigma), &(*R_)(i+1,i), 1, &(*R_)(i+1,curDim+j), 1);
822  (*R_)(i,curDim+j) -= sigma;
823  }
824  //
825  // Compute new Householder reflector
826  //
827  maxidx = blas.IAMAX( blockSize_+1, &(*R_)(curDim+j,curDim+j), 1 );
828  maxelem = (*R_)(curDim+j+maxidx-1,curDim+j);
829  for (i=0; i<blockSize_+1; i++)
830  (*R_)(curDim+j+i,curDim+j) /= maxelem;
831  sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j), 1,
832  &(*R_)(curDim+j+1,curDim+j), 1 );
833  if (sigma == zero) {
834  beta[curDim + j] = zero;
835  } else {
836  mu = Teuchos::ScalarTraits<ScalarType>::squareroot((*R_)(curDim+j,curDim+j)*(*R_)(curDim+j,curDim+j)+sigma);
837  if ( Teuchos::ScalarTraits<ScalarType>::real((*R_)(curDim+j,curDim+j))
839  vscale = (*R_)(curDim+j,curDim+j) - mu;
840  } else {
841  vscale = -sigma / ((*R_)(curDim+j,curDim+j) + mu);
842  }
843  beta[curDim+j] = Teuchos::as<ScalarType>(2)*one*vscale*vscale/(sigma + vscale*vscale);
844  (*R_)(curDim+j,curDim+j) = maxelem*mu;
845  for (i=0; i<blockSize_; i++)
846  (*R_)(curDim+j+1+i,curDim+j) /= vscale;
847  }
848  //
849  // Apply new Householder reflector to rhs
850  //
851  for (i=0; i<blockSize_; i++) {
852  sigma = blas.DOT( blockSize_, &(*R_)(curDim+j+1,curDim+j),
853  1, &(*z_)(curDim+j+1,i), 1);
854  sigma += (*z_)(curDim+j,i);
855  sigma *= beta[curDim+j];
856  blas.AXPY(blockSize_, ScalarType(-sigma), &(*R_)(curDim+j+1,curDim+j),
857  1, &(*z_)(curDim+j+1,i), 1);
858  (*z_)(curDim+j,i) -= sigma;
859  }
860  }
861  } // end if (blockSize_ == 1)
862 
863  // If the least-squares problem is updated wrt "dim" then update the curDim_.
864  if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
865  curDim_ = dim + blockSize_;
866  }
867 
868  } // end updateLSQR()
869 
870 } // end Belos namespace
871 
872 #endif /* BELOS_BLOCK_GMRES_ITER_HPP */
ScalarType * values() const
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...
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Array< T > & append(const T &x)
Class which manages the output and verbosity of the Belos solvers.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
T & get(const std::string &name, T def_value)
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
OrdinalType IAMAX(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Structure to contain pointers to GmresIteration state variables.
void resetNumIters(int iter=0)
Reset the iteration count.
SCT::magnitudeType MagnitudeType
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
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...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
virtual ~BlockGmresIter()
Destructor.
Traits class which defines basic operations on multivectors.
bool isParameter(const std::string &name) const
bool isInitialized()
States whether the solver has been initialized or not.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int curDim
The current dimension of the reduction.
OperatorTraits< ScalarType, MV, OP > OPT
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.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
ScalarType DOT(const OrdinalType &n, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy) const
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Teuchos::ScalarTraits< ScalarType > SCT
GmresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
OrdinalType numCols() const
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Class which defines basic traits for the operator type.
void initializeGmres(GmresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
BlockGmresIter(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)
BlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver option...
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...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
void setBlockSize(int blockSize)
Set the blocksize.
OrdinalType stride() const
OrdinalType numRows() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.

Generated on Thu Mar 28 2024 09:27:40 for Belos by doxygen 1.8.5