Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBlockGCRODRIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_BLOCK_GCRODR_ITER_HPP
11 #define BELOS_BLOCK_GCRODR_ITER_HPP
12 
13 
18 #include "BelosConfigDefs.hpp"
19 #include "BelosTypes.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosMatOrthoManager.hpp"
23 #include "BelosOutputManager.hpp"
24 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
28 #include "Teuchos_BLAS.hpp"
31 #include "Teuchos_ScalarTraits.hpp"
33 #include "Teuchos_TimeMonitor.hpp"
34 
35 // MLP
36 #include <unistd.h>
37 
52 namespace Belos{
53 
55 
56 
62  template <class ScalarType, class MV>
69  int curDim;
70 
73 
76 
83 
87 
88  BlockGCRODRIterState() : curDim(0), V(Teuchos::null),
89  U(Teuchos::null), C(Teuchos::null),
90  H(Teuchos::null), B(Teuchos::null)
91  {}
92 
93  };
94 
96 
97 
98 
100 
101 
115  public:
116  BlockGCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
117  };
118 
127  public:
128  BlockGCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
129  };
130 
132 
133 
134  template<class ScalarType, class MV, class OP>
135  class BlockGCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
136  public:
137 
138  //
139  //Convenience typedefs
140  //
147 
149 
150 
160  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
163  Teuchos::ParameterList &params );
164 
166  virtual ~BlockGCRODRIter() {};
168 
170 
171 
193  void iterate();
194 
217  void initialize() {
219  initialize(empty);
220  }
221 
226 
236  state.curDim = curDim_;
237  state.V = V_;
238  state.U = U_;
239  state.C = C_;
240  state.H = H_;
241  state.B = B_;
242  return state;
243  }
245 
247 
248 
250  bool isInitialized(){ return initialized_;};
251 
253  int getNumIters() const { return iter_; };
254 
256  void resetNumIters( int iter = 0 ) { iter_ = iter; };
257 
260  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
261 
263 
269 
270 
272 
273 
275 
276 
277 
279  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; };
280 
282  int getNumBlocks() const { return numBlocks_; }
283 
285  int getBlockSize() const { return blockSize_; };
286 
288  int getCurSubspaceDim() const {
289  if (!initialized_) return 0;
290  return curDim_;
291  };
292 
294  int getMaxSubspaceDim() const { return numBlocks_*blockSize_; };
295 
297  int getRecycledBlocks() const { return recycledBlocks_; };
298 
300 
301 
303 
304 
305  void updateLSQR( int dim = -1);
306 
308  void setBlockSize(int blockSize){ blockSize_ = blockSize; }
309 
311  void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
312 
314  void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
315 
317  void setSize( int recycledBlocks, int numBlocks ) {
318  // only call resize if size changed
319  if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
320  recycledBlocks_ = recycledBlocks;
321  numBlocks_ = numBlocks;
322  cs_.sizeUninitialized( numBlocks_ );
323  sn_.sizeUninitialized( numBlocks_ );
324  Z_.shapeUninitialized( numBlocks_*blockSize_, blockSize_ );
325  }
326  }
327 
329 
330  private:
331 
332  //
333  // Internal methods
334  //
335 
336 
337 
338  //Classes inputed through constructor that define the linear problem to be solved
339  //
344 
345  //
346  //Algorithmic Parameters
347  //
348 
349  //numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
350  //blockSize_ is the number of columns in each block Krylov vector.
351  int numBlocks_, blockSize_;
352 
353  //boolean vector indicating which right-hand sides we care about
354  //when we are testing residual norms. THIS IS NOT IMPLEMENTED. RIGHT NOW JUST
355  //SELECTS ALL RIGHT HANDS SIDES FOR NORM COMPUTATION.
356  std::vector<bool> trueRHSIndices_;
357 
358  // recycledBlocks_ is the size of the allocated space for the recycled subspace, in vectors.
359  int recycledBlocks_;
360 
361  //Storage for QR factorization of the least squares system if using plane rotations.
362  SDV sn_;
364 
365  //Storage for QR factorization of the least squares system if using Householder reflections
366  //Per block Krylov vector, we actually construct a 2*blockSize_ by 2*blockSize_ matrix which
367  //is the product of all Householder transformations for that block. This has been shown to yield
368  //speed ups without losing accuracy because we can apply previous Householder transformations
369  //with BLAS3 operations.
370  std::vector< SDM >House_;
371  SDV beta_;
372 
373  //
374  //Current Solver State
375  //
376  //initialized_ specifies that the basis vectors have been initialized and the iterate() routine
377  //is capable of running; _initialize is controlled by the initialize() member method
378  //For the implications of the state of initialized_, please see documentation for initialize()
379  bool initialized_;
380 
381  // Current subspace dimension, number of iterations performed, and number of iterations performed in this cycle.
382  int curDim_, iter_, lclIter_;
383 
384  //
385  // Recycled Krylov Method Storage
386  //
387 
388 
390  Teuchos::RCP<MV> V_;
391 
393  Teuchos::RCP<MV> U_, C_;
394 
395 
397 
398 
403 
408 
416 
418  SDM Z_;
419 
421 
422  // File stream variables to use Mike Parks' Matlab output codes.
423  std::ofstream ofs;
424  char filename[30];
425 
426  };//End BlockGCRODRIter Class Definition
427 
429  //Constructor.
430  template<class ScalarType, class MV, class OP>
432  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
435  Teuchos::ParameterList &params ):lp_(problem),
436  om_(printer), stest_(tester), ortho_(ortho) {
437  numBlocks_ = 0;
438  blockSize_ = 0;
439  recycledBlocks_ = 0;
440  initialized_ = false;
441  curDim_ = 0;
442  iter_ = 0;
443  lclIter_ = 0;
444  V_ = Teuchos::null;
445  U_ = Teuchos::null;
446  C_ = Teuchos::null;
447  H_ = Teuchos::null;
448  B_ = Teuchos::null;
449  R_ = Teuchos::null;
450  // Get the maximum number of blocks allowed for this Krylov subspace
451  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::BlockGCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
452  int nb = Teuchos::getParameter<int>(params, "Num Blocks");
453 
454  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
455  int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
456 
457  TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
458  TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::BlockGCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
459 
460 
461  int bs = Teuchos::getParameter<int>(params, "Block Size");
462 
463  TEUCHOS_TEST_FOR_EXCEPTION(bs <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() the block size was passed a non-postitive argument.");
464 
465 
466  numBlocks_ = nb;
467  recycledBlocks_ = rb;
468  blockSize_ = bs;
469 
470  //INITIALIZE ENTRIES OF trueRHSIndices_ TO CHECK EVERY NORM FOR NOW. LATER, THE USER
471  //SHOULD SPECIFY WHICH RIGHT HAND SIDES ARE IMPORTANT FOR CONVERGENCE TESTING
472  trueRHSIndices_.resize(blockSize_);
473  int i;
474  for(i=0; i<blockSize_; i++){
475  trueRHSIndices_[i] = true;
476  }
477 
478  //THIS MAKES SPACE FOR GIVENS ROTATIONS BUT IN REALITY WE NEED TO DO TESTING ON BLOCK SIZE
479  //AND CHOOSE BETWEEN GIVENS ROTATIONS AND HOUSEHOLDER TRANSFORMATIONS.
480  cs_.sizeUninitialized( numBlocks_+1 );
481  sn_.sizeUninitialized( numBlocks_+1 );
482  Z_.shapeUninitialized( (numBlocks_+1)*blockSize_,blockSize_ );
483 
484  House_.resize(numBlocks_);
485 
486  for(i=0; i<numBlocks_;i++){
487  House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
488  }
489  }//End Constructor Definition
490 
492  // Iterate until the status test informs us we should stop.
493  template <class ScalarType, class MV, class OP>
495  TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, BlockGCRODRIterInitFailure,"Belos::BlockGCRODRIter::iterate(): GCRODRIter class not initialized." );
496 
497 // MLP
498 //sleep(1);
499 //std::cout << "Calling setSize" << std::endl;
500  // Force call to setsize to ensure internal storage is correct dimension
501  setSize( recycledBlocks_, numBlocks_ );
502 
503  Teuchos::RCP<MV> Vnext;
505  std::vector<int> curind(blockSize_);
506 
507  // z_ must be zeroed out in order to compute Givens rotations correctly
508  Z_.putScalar(0.0);
509 
510  // Orthonormalize the new V_0
511  for(int i = 0; i<blockSize_; i++){curind[i] = i;};
512 
513 // MLP
514 //sleep(1);
515 //std::cout << "Calling normalize" << std::endl;
516  Vnext = MVT::CloneViewNonConst(*V_,curind);
517  //Orthonormalize Initial Columns
518  //Store orthogonalization coefficients in Z0
519  Teuchos::RCP<SDM > Z0 =
520  Teuchos::rcp( new SDM(blockSize_,blockSize_) );
521  int rank = ortho_->normalize(*Vnext,Z0);
522 
523 // MLP
524 //sleep(1);
525 //std::cout << "Assigning Z" << std::endl;
526  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank at the initial step.");
527  // Copy Z0 into the leading blockSize_ by blockSize_ block of Z_
528  Teuchos::RCP<SDM > Z_block = Teuchos::rcp( new SDM(Teuchos::View, Z_, blockSize_,blockSize_) );
529  Z_block->assign(*Z0);
530 
531  std::vector<int> prevind(blockSize_*(numBlocks_ + 1));
532 
534  // iterate until the status test tells us to stop.
535  //
536  // also break if the basis is full
537  //
538  while( (stest_->checkStatus(this) != Passed) && (curDim_+blockSize_-1) < (numBlocks_*blockSize_)) {
539  lclIter_++;
540  iter_++;
541 //KMS
542 //std::cout << "Iter=" << iter_ << std::endl << "lclIter=" << lclIter_ << std::endl;
543 
544  int HFirstCol = curDim_-blockSize_;//First column of H we need view of
545  int HLastCol = HFirstCol + blockSize_-1 ;//last column of H we need a view of
546  int HLastOrthRow = HLastCol;//Last row of H we will put orthog coefficients in
547  int HFirstNormRow = HLastOrthRow + 1;//First row of H where normalization matrix goes
548 //KMS
549 //std::cout << "curDim_ = " << curDim_ << ", HFirstCol = " << HFirstCol << ", HLastCol = " << HLastCol <<", HLastOrthRow = " << HLastOrthRow << ", HFirstNormRow = " << HFirstNormRow << std::endl;
550  // Get next basis indices
551  for(int i = 0; i< blockSize_; i++){
552  curind[i] = curDim_ + i;
553  }
554  Vnext = MVT::CloneViewNonConst(*V_,curind);
555 
556  //Get a view of the previous block Krylov vector.
557  //This is used for orthogonalization and for computing V^H K H
558  // Get next basis indices
559  for(int i = 0; i< blockSize_; i++){
560  curind[blockSize_ - 1 - i] = curDim_ - i - 1;
561  }
562  Vprev = MVT::CloneView(*V_,curind);
563  // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
564  lp_->apply(*Vprev,*Vnext);
565  Vprev = Teuchos::null;
566 
567  //First, remove the recycled subspace (C) from Vnext and put coefficients in B.
568 
569  //Get a view of the matrix B and put the pointer into an array
570  //Put a pointer to C in another array
572 
574  subB = Teuchos::rcp( new SDM ( Teuchos::View,*B_,recycledBlocks_,blockSize_,0, HFirstCol ) );
575 
577  AsubB.append( subB );
578  // Project out the recycled subspace.
579  ortho_->project( *Vnext, AsubB, C );
580  //Now, remove block Krylov Subspace from Vnext and store coefficients in H_ and R_
581 
582  // Get a view of all the previous vectors
583  prevind.resize(curDim_);
584  for (int i=0; i<curDim_; i++) { prevind[i] = i; }
585  Vprev = MVT::CloneView(*V_,prevind);
586  Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
587 
588  // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
589  Teuchos::RCP<SDM> subH = Teuchos::rcp( new SDM ( Teuchos::View,*H_,curDim_,blockSize_,0,HFirstCol ) );
591  AsubH.append( subH );
592  // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
593  Teuchos::RCP<SDM > subR = Teuchos::rcp( new SDM ( Teuchos::View,*H_,blockSize_,blockSize_,HFirstNormRow,HFirstCol ) );
594  // Project out the previous Krylov vectors and normalize the next vector.
595  int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
596 
597  // Copy over the coefficients to R just in case we run into an error.
598  SDM subR2( Teuchos::View,*R_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
599  SDM subH2( Teuchos::View,*H_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
600  subR2.assign(subH2);
601 
602  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank.");
603 
604  // Update the QR factorization of the upper Hessenberg matrix
605  updateLSQR();
606 
607  // Update basis dim
608  curDim_ = curDim_ + blockSize_;
609 
610 
611 
612  }//end while(stest_->checkStatus(this) ~= Passed && curDim_+1 <= numBlocks_*blockSize_)
613 
614  }//end iterate() defintition
615 
617  //Initialize this iteration object.
618  template <class ScalarType, class MV, class OP>
620  if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
621  curDim_ = newstate.curDim;
622  V_ = newstate.V;
623  U_ = newstate.U;
624  C_ = newstate.C;
625  H_ = newstate.H;
626  B_ = newstate.B;
627  lclIter_ = 0;//resets the count of local iterations for the new cycle
628  R_ = Teuchos::rcp(new SDM(H_->numRows(), H_->numCols() )); //R_ should look like H but point to separate memory
629 
630  //All Householder product matrices start out as identity matrices.
631  //We construct an identity from which to copy.
632  SDM Identity(2*blockSize_, 2*blockSize_);
633  for(int i=0;i<2*blockSize_; i++){
634  Identity[i][i] = 1;
635  }
636  for(int i=0; i<numBlocks_;i++){
637  House_[i].assign(Identity);
638  }
639  }
640  else {
641  TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have V initialized.");
642  TEUCHOS_TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have H initialized.");
643  }
644  // the solver is initialized
645  initialized_ = true;
646  }//end initialize() defintition
647 
649  //Get the native residuals stored in this iteration.
650  //This function will only compute the native residuals for
651  //right-hand sides we are interested in, as dictated by
652  //std::vector<int> trueRHSIndices_ (THIS IS NOT YET IMPLEMENTED. JUST GETS ALL RESIDUALS)
653  //A norm of -1 is entered for all residuals about which we do not care.
654  template <class ScalarType, class MV, class OP>
656  BlockGCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const
657  {
658  //
659  // NOTE: Make sure the incoming std::vector is the correct size!
660  //
661  if (norms != NULL) {
662  if (static_cast<int> (norms->size()) < blockSize_) {
663  norms->resize( blockSize_ );
664  }
666  for (int j=0; j<blockSize_; j++) {
667  if(trueRHSIndices_[j]){
668  (*norms)[j] = blas.NRM2( blockSize_, &Z_(curDim_-blockSize_+j, j), 1);
669  }
670  else{
671  (*norms)[j] = -1;
672  }
673  }
674  return Teuchos::null;
675  } else { // norms is NULL
676  // FIXME If norms is NULL, return residual vectors.
677  return Teuchos::null;
678  }
679  }//end getNativeResiduals() definition
680 
682  //Get the current update from this subspace.
683  template <class ScalarType, class MV, class OP>
685  //
686  // If this is the first iteration of the Arnoldi factorization,
687  // there is no update, so return Teuchos::null.
688  //
689  Teuchos::RCP<MV> currentUpdate = Teuchos::null;
690 //KMS if(curDim_==0) {
691  if(curDim_<=blockSize_) {
692  return currentUpdate;
693  }
694  else{
695  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
696  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
698  currentUpdate = MVT::Clone( *V_, blockSize_ );
699  //
700  // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
701  //
702  SDM Y( Teuchos::Copy, Z_, curDim_-blockSize_, blockSize_ );
703  Teuchos::RCP<SDM> Rtmp = Teuchos::rcp(new SDM(Teuchos::View, *R_, curDim_, curDim_-blockSize_));
704 //KMS
705 //sleep(1);
706 //std::cout << "Before TRSM" << std::endl;
707 //sleep(1);
708 //std::cout << "The size of Rtmp is " << Rtmp -> numRows() << " by " << Rtmp -> numCols() << std::endl;
709 //std::cout << "The size of Y is " << Y.numRows() << " by " << Y.numCols() << std::endl;
710 //std::cout << "blockSize_ = " << blockSize_ << std::endl;
711 //std::cout << "curDim_ = " << curDim_ << std::endl;
712 //std::cout << "curDim_ - blockSize_ = " << curDim_ - blockSize_ << std::endl;
713  //
714  // Solve the least squares problem.
715  // Observe that in calling TRSM, we use the value
716  // curDim_ -blockSize_. This is because curDim_ has
717  // already been incremented but the proper size of R is still
718  // based on the previous value.
719  //
721  Teuchos::NON_UNIT_DIAG, curDim_-blockSize_, blockSize_, one,
722  Rtmp->values(), Rtmp->stride(), Y.values(), Y.stride() );
723 //KMS
724 //sleep(1);
725 //std::cout << "After TRSM" << std::endl;
726  //
727  // Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
728  //
729  std::vector<int> index(curDim_-blockSize_);
730  for ( int i=0; i<curDim_-blockSize_; i++ ) index[i] = i;
731  Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
732  MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
733 
734 
735 
736  //
737  // Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
738  //
739  if (U_ != Teuchos::null) {
740  SDM z(recycledBlocks_,blockSize_);
741  SDM subB( Teuchos::View, *B_, recycledBlocks_, curDim_-blockSize_ );
742  z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, Y, zero );
743 
744  //std::cout << (*U_).MyLength() << " " << (*U_).NumVectors() << " " << subB.numRows() << " " << subB.numCols() << " " << Y.numRows() << " " << Y.numCols()<< " " << curDim_ << " " << recycledBlocks_;
745  MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
746  }
747  }
748 
749 
750 
751  return currentUpdate;
752  }//end getCurrentUpdate() definition
753 
754  template<class ScalarType, class MV, class OP>
756 
757  int i;
758  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
759  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
760 
761  using Teuchos::rcp;
762 
763  // Get correct dimension based on input "dim"
764  // Remember that ortho failures result in an exit before updateLSQR() is called.
765  // Therefore, it is possible that dim == curDim_.
766  int curDim = curDim_;
767  if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
768  curDim = dim;
769  }
770 
772 
773  if(blockSize_ == 1){//if only one right-hand side then use Givens rotations
774  //
775  // Apply previous rotations and compute new rotations to reduce upper-Hessenberg
776  // system to upper-triangular form.
777  //
778  // QR factorization of Least-Squares system with Givens rotations
779  //
780  for (i=0; i<curDim-1; i++) {
781  //
782  // Apply previous Givens rotations to new column of Hessenberg matrix
783  //
784  blas.ROT( 1, &(*R_)(i,curDim-1), 1, &(*R_)(i+1, curDim-1), 1, &cs_[i], &sn_[i] );
785  }
786  //
787  // Calculate new Givens rotation
788  //
789  blas.ROTG( &(*R_)(curDim-1,curDim-1), &(*R_)(curDim,curDim-1), &cs_[curDim-1], &sn_[curDim-1] );
790  (*R_)(curDim,curDim-1) = zero;
791  //
792  // Update RHS w/ new transformation
793  //
794  blas.ROT( 1, &Z_(curDim-1,0), 1, &Z_(curDim,0), 1, &cs_[curDim-1], &sn_[curDim-1] );
795  }
796  else{//if multiple right-hand sides then use Householder transormations
797  //
798  //apply previous reflections and compute new reflections to reduce upper-Hessenberg
799  //system to upper-triagular form.
800 
801  //In Matlab, applying the reflection to a matrix M would look like
802  // M_refl = M - 2*v_refl*(v_refl'*M)/(norm(v_refl)^2)
803 
804  //In order to take advantage of BLAS while applying reflections to a matrix, we
805  //perform it in a two step process
806  //1. workvec = M'*v_refl {using BLAS.GEMV()}
807  //2. M_refl = M_refl - 2*v_refl*workvec'/(norm(v_refl)^2) {using BLAS.GER()}
808 
809  Teuchos::RCP< SDM > workmatrix = Teuchos::null;//matrix of column vectors onto which we apply the reflection
810  Teuchos::RCP< SDV > workvec = Teuchos::null;//where we store the result of the first step of the 2-step reflection process
811  Teuchos::RCP<SDV> v_refl = Teuchos::null;//the reflection vector
812  int R_colStart = curDim_-blockSize_;
813  Teuchos::RCP< SDM >Rblock = Teuchos::null;
814 
815  //
816  //Apply previous reflections
817  //
818  for(i=0; i<lclIter_-1; i++){
819  int R_rowStart = i*blockSize_;
820  //get a view of the part of R_ effected by these reflections.
821  Teuchos::RCP< SDM > RblockCopy = rcp(new SDM (Teuchos::Copy, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
822  Teuchos::RCP< SDM > RblockView = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
823  blas.GEMM(Teuchos::NO_TRANS,Teuchos::NO_TRANS, 2*blockSize_,blockSize_,2*blockSize_,one,House_[i].values(),House_[i].stride(), RblockCopy->values(),RblockCopy -> stride(), zero, RblockView->values(),RblockView -> stride());
824 
825  }
826 
827 
828  //Get a view of last 2*blockSize entries of entire block to
829  //generate new reflections.
830  Rblock = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, curDim_-blockSize_, curDim_-blockSize_));
831 
832  //Calculate and apply the new reflections
833  for(i=0; i<blockSize_; i++){
834  //
835  //Calculating Reflection
836  //
837  int curcol = (lclIter_ - 1)*blockSize_ + i;//current column of R_
838  int lclCurcol = i;//current column of Rblock
839  ScalarType signDiag = (*R_)(curcol,curcol) / Teuchos::ScalarTraits<ScalarType>::magnitude((*R_)(curcol,curcol));
840 
841  // Norm of the vector to be reflected.
842  // BLAS returns a ScalarType, but it really should be a magnitude.
843  ScalarType nvs = blas.NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
844  ScalarType alpha = -signDiag*nvs;
845 
846  //norm of reflection vector which is just the vector being reflected
847  //i.e. v = R_(curcol:curcol+blockSize_,curcol))
848  //v_refl = v - alpha*e1
849  //norm(v_refl) = norm(v) + alpha^2 - 2*v*alpha
850  //store in v_refl
851 
852  // Beware, nvs should have a zero imaginary part (since
853  // it is a norm of a vector), but it may not due to rounding
854  // error.
855  //nvs = nvs + alpha*alpha - 2*(*R_)(curcol,curcol)*alpha;
856  //(*R_)(curcol,curcol) -= alpha;
857 
858  //Copy relevant values of the current column of R_ into the reflection vector
859  //Modify first entry
860  //Take norm of reflection vector
861  //Square the norm
862  v_refl = rcp(new SDV(Teuchos::Copy, &((*R_)(curcol,curcol)), blockSize_ + 1 ));
863  (*v_refl)[0] -= alpha;
864  nvs = blas.NRM2(blockSize_+1,v_refl -> values(),1);
865  nvs *= nvs;
866 
867  //
868  //Apply new reflector to:
869  //1. To subsequent columns of R_ in the current block
870  //2. To House[iter_] to store product of reflections for this column
871  //3. To the least-squares right-hand side.
872  //4. The current column
873  //
874  //
875 
876 
877  //
878  //1.
879  //
880  if(i < blockSize_ - 1){//only do this when there are subsquent columns in the block to apply to
881  workvec = Teuchos::rcp(new SDV(blockSize_ - i -1));
882  //workvec = Teuchos::rcp(new SDV(2*blockSize_));
883  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, *Rblock, blockSize_+1, blockSize_ - i -1, lclCurcol, lclCurcol +1 ) );
884  blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix->values(), workmatrix->stride(), v_refl->values(), 1, zero, workvec->values(), 1);
885  blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl->values(),1,workvec->values(),1,workmatrix->values(),workmatrix->stride());
886  }
887 
888 
889  //
890  //2.
891  //
892  workvec = Teuchos::rcp(new SDV(2*blockSize_));
893  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, House_[lclIter_ -1], blockSize_+1, 2*blockSize_, i, 0 ) );
894  blas.GEMV(Teuchos::TRANS,workmatrix->numRows(),workmatrix->numCols(),one,workmatrix->values(),workmatrix->stride(), v_refl->values(), 1,zero,workvec->values(),1);
895  blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl -> values(),1,workvec->values(),1,workmatrix->values(),(*workmatrix).stride());
896 
897  //
898  //3.
899  //
900  workvec = Teuchos::rcp(new SDV(blockSize_));
901  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, Z_, blockSize_+1, blockSize_, curcol, 0 ) );
902  blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix-> values(), workmatrix->stride(), v_refl -> values(), 1, zero, workvec->values(), 1);
903  blas.GER((*workmatrix).numRows(),(*workmatrix).numCols(), -2.0*one/nvs,v_refl -> values(), 1,&((*workvec)[0]),1,(*workmatrix)[0],(*workmatrix).stride());
904 
905  //
906  //4.
907  //
908  (*R_)[curcol][curcol] = alpha;
909  for(int ii=1; ii<= blockSize_; ii++){
910  (*R_)[curcol][curcol+ii] = 0;
911  }
912  }
913 
914  }
915 
916  } // end updateLSQR()
917 
918 
919 }//End Belos Namespace
920 
921 #endif /* BELOS_BLOCK_GCRODR_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
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
void GER(const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy, ScalarType *A, const OrdinalType &lda) const
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Array< T > & append(const T &x)
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *.
virtual ~BlockGCRODRIter()
Destructor.
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors...
void iterate()
This method performs block GCRODR iterations until the status test indicates the need to stop or an e...
MultiVecTraits< ScalarType, MV > MVT
Pure virtual base class for defining the status testing capabilities of Belos.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
A pure virtual class for defining the status tests for the Belos iterative solvers.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Class which defines basic traits for the operator type.
BlockGCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
Teuchos::RCP< MV > V
The current Krylov basis.
Traits class which defines basic operations on multivectors.
int curDim
The current dimension of the reduction.
SCT::magnitudeType MagnitudeType
bool isParameter(const std::string &name) const
BlockGCRODRIterOrthoFailure(const std::string &what_arg)
BlockGCRODRIter(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)
BlockGCRODRIter constructor with linear problem, solver utilities, and parameter list of solver optio...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setBlockSize(int blockSize)
Set the blocksize.
BlockGCRODRIterInitFailure(const std::string &what_arg)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
void resetNumIters(int iter=0)
Reset the iteration count.
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
void initialize()
Initialize the solver to an iterate, providing a complete state.
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
BlockGCRODRIterInitFailure is thrown when the BlockGCRODRIter object is unable to generate an initial...
int getRecycledBlocks() const
Set the maximum number of recycled blocks used by the iterative solver.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
bool isInitialized()
States whether the solver has been initialized or not.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Teuchos::SerialDenseVector< int, ScalarType > SDV
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...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos::ScalarTraits< ScalarType > SCT
OperatorTraits< ScalarType, MV, OP > OPT
int getNumIters() const
Get the current iteration count.
OrdinalType stride() const
int sizeUninitialized(OrdinalType length_in)
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated on Mon Jul 15 2024 09:24:23 for Belos by doxygen 1.8.5