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 //
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_GCRODR_ITER_HPP
43 #define BELOS_BLOCK_GCRODR_ITER_HPP
44 
45 
50 #include "BelosConfigDefs.hpp"
51 #include "BelosTypes.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 
67 // MLP
68 #include <unistd.h>
69 
84 namespace Belos{
85 
87 
88 
94  template <class ScalarType, class MV>
101  int curDim;
102 
105 
108 
115 
119 
120  BlockGCRODRIterState() : curDim(0), V(Teuchos::null),
121  U(Teuchos::null), C(Teuchos::null),
122  H(Teuchos::null), B(Teuchos::null)
123  {}
124 
125  };
126 
128 
129 
130 
132 
133 
147  public:
148  BlockGCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
149  };
150 
159  public:
160  BlockGCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
161  };
162 
164 
165 
166  template<class ScalarType, class MV, class OP>
167  class BlockGCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
168  public:
169 
170  //
171  //Convenience typedefs
172  //
179 
181 
182 
192  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
195  Teuchos::ParameterList &params );
196 
198  virtual ~BlockGCRODRIter() {};
200 
202 
203 
225  void iterate();
226 
249  void initialize() {
251  initialize(empty);
252  }
253 
258 
268  state.curDim = curDim_;
269  state.V = V_;
270  state.U = U_;
271  state.C = C_;
272  state.H = H_;
273  state.B = B_;
274  return state;
275  }
277 
279 
280 
282  bool isInitialized(){ return initialized_;};
283 
285  int getNumIters() const { return iter_; };
286 
288  void resetNumIters( int iter = 0 ) { iter_ = iter; };
289 
292  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
293 
295 
301 
302 
304 
305 
307 
308 
309 
311  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; };
312 
314  int getNumBlocks() const { return numBlocks_; }
315 
317  int getBlockSize() const { return blockSize_; };
318 
320  int getCurSubspaceDim() const {
321  if (!initialized_) return 0;
322  return curDim_;
323  };
324 
326  int getMaxSubspaceDim() const { return numBlocks_*blockSize_; };
327 
329  int getRecycledBlocks() const { return recycledBlocks_; };
330 
332 
333 
335 
336 
337  void updateLSQR( int dim = -1);
338 
340  void setBlockSize(int blockSize){ blockSize_ = blockSize; }
341 
343  void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
344 
346  void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
347 
349  void setSize( int recycledBlocks, int numBlocks ) {
350  // only call resize if size changed
351  if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
352  recycledBlocks_ = recycledBlocks;
353  numBlocks_ = numBlocks;
354  cs_.sizeUninitialized( numBlocks_ );
355  sn_.sizeUninitialized( numBlocks_ );
356  Z_.shapeUninitialized( numBlocks_*blockSize_, blockSize_ );
357  }
358  }
359 
361 
362  private:
363 
364  //
365  // Internal methods
366  //
367 
368 
369 
370  //Classes inputed through constructor that define the linear problem to be solved
371  //
376 
377  //
378  //Algorithmic Parameters
379  //
380 
381  //numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
382  //blockSize_ is the number of columns in each block Krylov vector.
383  int numBlocks_, blockSize_;
384 
385  //boolean vector indicating which right-hand sides we care about
386  //when we are testing residual norms. THIS IS NOT IMPLEMENTED. RIGHT NOW JUST
387  //SELECTS ALL RIGHT HANDS SIDES FOR NORM COMPUTATION.
388  std::vector<bool> trueRHSIndices_;
389 
390  // recycledBlocks_ is the size of the allocated space for the recycled subspace, in vectors.
391  int recycledBlocks_;
392 
393  //Storage for QR factorization of the least squares system if using plane rotations.
394  SDV sn_;
396 
397  //Storage for QR factorization of the least squares system if using Householder reflections
398  //Per block Krylov vector, we actually construct a 2*blockSize_ by 2*blockSize_ matrix which
399  //is the product of all Householder transformations for that block. This has been shown to yield
400  //speed ups without losing accuracy because we can apply previous Householder transformations
401  //with BLAS3 operations.
402  std::vector< SDM >House_;
403  SDV beta_;
404 
405  //
406  //Current Solver State
407  //
408  //initialized_ specifies that the basis vectors have been initialized and the iterate() routine
409  //is capable of running; _initialize is controlled by the initialize() member method
410  //For the implications of the state of initialized_, please see documentation for initialize()
411  bool initialized_;
412 
413  // Current subspace dimension, number of iterations performed, and number of iterations performed in this cycle.
414  int curDim_, iter_, lclIter_;
415 
416  //
417  // Recycled Krylov Method Storage
418  //
419 
420 
422  Teuchos::RCP<MV> V_;
423 
425  Teuchos::RCP<MV> U_, C_;
426 
427 
429 
430 
435 
440 
448 
450  SDM Z_;
451 
453 
454  // File stream variables to use Mike Parks' Matlab output codes.
455  std::ofstream ofs;
456  char filename[30];
457 
458  };//End BlockGCRODRIter Class Definition
459 
461  //Constructor.
462  template<class ScalarType, class MV, class OP>
464  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
467  Teuchos::ParameterList &params ):lp_(problem),
468  om_(printer), stest_(tester), ortho_(ortho) {
469  numBlocks_ = 0;
470  blockSize_ = 0;
471  recycledBlocks_ = 0;
472  initialized_ = false;
473  curDim_ = 0;
474  iter_ = 0;
475  lclIter_ = 0;
476  V_ = Teuchos::null;
477  U_ = Teuchos::null;
478  C_ = Teuchos::null;
479  H_ = Teuchos::null;
480  B_ = Teuchos::null;
481  R_ = Teuchos::null;
482  // Get the maximum number of blocks allowed for this Krylov subspace
483  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::BlockGCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
484  int nb = Teuchos::getParameter<int>(params, "Num Blocks");
485 
486  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
487  int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
488 
489  TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
490  TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::BlockGCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
491 
492 
493  int bs = Teuchos::getParameter<int>(params, "Block Size");
494 
495  TEUCHOS_TEST_FOR_EXCEPTION(bs <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() the block size was passed a non-postitive argument.");
496 
497 
498  numBlocks_ = nb;
499  recycledBlocks_ = rb;
500  blockSize_ = bs;
501 
502  //INITIALIZE ENTRIES OF trueRHSIndices_ TO CHECK EVERY NORM FOR NOW. LATER, THE USER
503  //SHOULD SPECIFY WHICH RIGHT HAND SIDES ARE IMPORTANT FOR CONVERGENCE TESTING
504  trueRHSIndices_.resize(blockSize_);
505  int i;
506  for(i=0; i<blockSize_; i++){
507  trueRHSIndices_[i] = true;
508  }
509 
510  //THIS MAKES SPACE FOR GIVENS ROTATIONS BUT IN REALITY WE NEED TO DO TESTING ON BLOCK SIZE
511  //AND CHOOSE BETWEEN GIVENS ROTATIONS AND HOUSEHOLDER TRANSFORMATIONS.
512  cs_.sizeUninitialized( numBlocks_+1 );
513  sn_.sizeUninitialized( numBlocks_+1 );
514  Z_.shapeUninitialized( (numBlocks_+1)*blockSize_,blockSize_ );
515 
516  House_.resize(numBlocks_);
517 
518  for(i=0; i<numBlocks_;i++){
519  House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
520  }
521  }//End Constructor Definition
522 
524  // Iterate until the status test informs us we should stop.
525  template <class ScalarType, class MV, class OP>
527  TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, BlockGCRODRIterInitFailure,"Belos::BlockGCRODRIter::iterate(): GCRODRIter class not initialized." );
528 
529 // MLP
530 //sleep(1);
531 //std::cout << "Calling setSize" << std::endl;
532  // Force call to setsize to ensure internal storage is correct dimension
533  setSize( recycledBlocks_, numBlocks_ );
534 
535  Teuchos::RCP<MV> Vnext;
537  std::vector<int> curind(blockSize_);
538 
539  // z_ must be zeroed out in order to compute Givens rotations correctly
540  Z_.putScalar(0.0);
541 
542  // Orthonormalize the new V_0
543  for(int i = 0; i<blockSize_; i++){curind[i] = i;};
544 
545 // MLP
546 //sleep(1);
547 //std::cout << "Calling normalize" << std::endl;
548  Vnext = MVT::CloneViewNonConst(*V_,curind);
549  //Orthonormalize Initial Columns
550  //Store orthogonalization coefficients in Z0
551  Teuchos::RCP<SDM > Z0 =
552  Teuchos::rcp( new SDM(blockSize_,blockSize_) );
553  int rank = ortho_->normalize(*Vnext,Z0);
554 
555 // MLP
556 //sleep(1);
557 //std::cout << "Assigning Z" << std::endl;
558  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank at the initial step.");
559  // Copy Z0 into the leading blockSize_ by blockSize_ block of Z_
560  Teuchos::RCP<SDM > Z_block = Teuchos::rcp( new SDM(Teuchos::View, Z_, blockSize_,blockSize_) );
561  Z_block->assign(*Z0);
562 
563  std::vector<int> prevind(blockSize_*(numBlocks_ + 1));
564 
566  // iterate until the status test tells us to stop.
567  //
568  // also break if the basis is full
569  //
570  while( (stest_->checkStatus(this) != Passed) && (curDim_+blockSize_-1) < (numBlocks_*blockSize_)) {
571  lclIter_++;
572  iter_++;
573 //KMS
574 //std::cout << "Iter=" << iter_ << std::endl << "lclIter=" << lclIter_ << std::endl;
575 
576  int HFirstCol = curDim_-blockSize_;//First column of H we need view of
577  int HLastCol = HFirstCol + blockSize_-1 ;//last column of H we need a view of
578  int HLastOrthRow = HLastCol;//Last row of H we will put orthog coefficients in
579  int HFirstNormRow = HLastOrthRow + 1;//First row of H where normalization matrix goes
580 //KMS
581 //std::cout << "curDim_ = " << curDim_ << ", HFirstCol = " << HFirstCol << ", HLastCol = " << HLastCol <<", HLastOrthRow = " << HLastOrthRow << ", HFirstNormRow = " << HFirstNormRow << std::endl;
582  // Get next basis indices
583  for(int i = 0; i< blockSize_; i++){
584  curind[i] = curDim_ + i;
585  }
586  Vnext = MVT::CloneViewNonConst(*V_,curind);
587 
588  //Get a view of the previous block Krylov vector.
589  //This is used for orthogonalization and for computing V^H K H
590  // Get next basis indices
591  for(int i = 0; i< blockSize_; i++){
592  curind[blockSize_ - 1 - i] = curDim_ - i - 1;
593  }
594  Vprev = MVT::CloneView(*V_,curind);
595  // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
596  lp_->apply(*Vprev,*Vnext);
597  Vprev = Teuchos::null;
598 
599  //First, remove the recycled subspace (C) from Vnext and put coefficients in B.
600 
601  //Get a view of the matrix B and put the pointer into an array
602  //Put a pointer to C in another array
604 
606  subB = Teuchos::rcp( new SDM ( Teuchos::View,*B_,recycledBlocks_,blockSize_,0, HFirstCol ) );
607 
609  AsubB.append( subB );
610  // Project out the recycled subspace.
611  ortho_->project( *Vnext, AsubB, C );
612  //Now, remove block Krylov Subspace from Vnext and store coefficients in H_ and R_
613 
614  // Get a view of all the previous vectors
615  prevind.resize(curDim_);
616  for (int i=0; i<curDim_; i++) { prevind[i] = i; }
617  Vprev = MVT::CloneView(*V_,prevind);
618  Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
619 
620  // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
621  Teuchos::RCP<SDM> subH = Teuchos::rcp( new SDM ( Teuchos::View,*H_,curDim_,blockSize_,0,HFirstCol ) );
623  AsubH.append( subH );
624  // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
625  Teuchos::RCP<SDM > subR = Teuchos::rcp( new SDM ( Teuchos::View,*H_,blockSize_,blockSize_,HFirstNormRow,HFirstCol ) );
626  // Project out the previous Krylov vectors and normalize the next vector.
627  int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
628 
629  // Copy over the coefficients to R just in case we run into an error.
630  SDM subR2( Teuchos::View,*R_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
631  SDM subH2( Teuchos::View,*H_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
632  subR2.assign(subH2);
633 
634  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank.");
635 
636  // Update the QR factorization of the upper Hessenberg matrix
637  updateLSQR();
638 
639  // Update basis dim
640  curDim_ = curDim_ + blockSize_;
641 
642 
643 
644  }//end while(stest_->checkStatus(this) ~= Passed && curDim_+1 <= numBlocks_*blockSize_)
645 
646  }//end iterate() defintition
647 
649  //Initialize this iteration object.
650  template <class ScalarType, class MV, class OP>
652  if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
653  curDim_ = newstate.curDim;
654  V_ = newstate.V;
655  U_ = newstate.U;
656  C_ = newstate.C;
657  H_ = newstate.H;
658  B_ = newstate.B;
659  lclIter_ = 0;//resets the count of local iterations for the new cycle
660  R_ = Teuchos::rcp(new SDM(H_->numRows(), H_->numCols() )); //R_ should look like H but point to separate memory
661 
662  //All Householder product matrices start out as identity matrices.
663  //We construct an identity from which to copy.
664  SDM Identity(2*blockSize_, 2*blockSize_);
665  for(int i=0;i<2*blockSize_; i++){
666  Identity[i][i] = 1;
667  }
668  for(int i=0; i<numBlocks_;i++){
669  House_[i].assign(Identity);
670  }
671  }
672  else {
673  TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have V initialized.");
674  TEUCHOS_TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have H initialized.");
675  }
676  // the solver is initialized
677  initialized_ = true;
678  }//end initialize() defintition
679 
681  //Get the native residuals stored in this iteration.
682  //This function will only compute the native residuals for
683  //right-hand sides we are interested in, as dictated by
684  //std::vector<int> trueRHSIndices_ (THIS IS NOT YET IMPLEMENTED. JUST GETS ALL RESIDUALS)
685  //A norm of -1 is entered for all residuals about which we do not care.
686  template <class ScalarType, class MV, class OP>
688  BlockGCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const
689  {
690  //
691  // NOTE: Make sure the incoming std::vector is the correct size!
692  //
693  if (norms != NULL) {
694  if (static_cast<int> (norms->size()) < blockSize_) {
695  norms->resize( blockSize_ );
696  }
698  for (int j=0; j<blockSize_; j++) {
699  if(trueRHSIndices_[j]){
700  (*norms)[j] = blas.NRM2( blockSize_, &Z_(curDim_-blockSize_+j, j), 1);
701  }
702  else{
703  (*norms)[j] = -1;
704  }
705  }
706  return Teuchos::null;
707  } else { // norms is NULL
708  // FIXME If norms is NULL, return residual vectors.
709  return Teuchos::null;
710  }
711  }//end getNativeResiduals() definition
712 
714  //Get the current update from this subspace.
715  template <class ScalarType, class MV, class OP>
717  //
718  // If this is the first iteration of the Arnoldi factorization,
719  // there is no update, so return Teuchos::null.
720  //
721  Teuchos::RCP<MV> currentUpdate = Teuchos::null;
722 //KMS if(curDim_==0) {
723  if(curDim_<=blockSize_) {
724  return currentUpdate;
725  }
726  else{
727  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
728  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
730  currentUpdate = MVT::Clone( *V_, blockSize_ );
731  //
732  // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
733  //
734  SDM Y( Teuchos::Copy, Z_, curDim_-blockSize_, blockSize_ );
735  Teuchos::RCP<SDM> Rtmp = Teuchos::rcp(new SDM(Teuchos::View, *R_, curDim_, curDim_-blockSize_));
736 //KMS
737 //sleep(1);
738 //std::cout << "Before TRSM" << std::endl;
739 //sleep(1);
740 //std::cout << "The size of Rtmp is " << Rtmp -> numRows() << " by " << Rtmp -> numCols() << std::endl;
741 //std::cout << "The size of Y is " << Y.numRows() << " by " << Y.numCols() << std::endl;
742 //std::cout << "blockSize_ = " << blockSize_ << std::endl;
743 //std::cout << "curDim_ = " << curDim_ << std::endl;
744 //std::cout << "curDim_ - blockSize_ = " << curDim_ - blockSize_ << std::endl;
745  //
746  // Solve the least squares problem.
747  // Observe that in calling TRSM, we use the value
748  // curDim_ -blockSize_. This is because curDim_ has
749  // already been incremented but the proper size of R is still
750  // based on the previous value.
751  //
753  Teuchos::NON_UNIT_DIAG, curDim_-blockSize_, blockSize_, one,
754  Rtmp->values(), Rtmp->stride(), Y.values(), Y.stride() );
755 //KMS
756 //sleep(1);
757 //std::cout << "After TRSM" << std::endl;
758  //
759  // Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
760  //
761  std::vector<int> index(curDim_-blockSize_);
762  for ( int i=0; i<curDim_-blockSize_; i++ ) index[i] = i;
763  Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
764  MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
765 
766 
767 
768  //
769  // Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
770  //
771  if (U_ != Teuchos::null) {
772  SDM z(recycledBlocks_,blockSize_);
773  SDM subB( Teuchos::View, *B_, recycledBlocks_, curDim_-blockSize_ );
774  z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, Y, zero );
775 
776  //std::cout << (*U_).MyLength() << " " << (*U_).NumVectors() << " " << subB.numRows() << " " << subB.numCols() << " " << Y.numRows() << " " << Y.numCols()<< " " << curDim_ << " " << recycledBlocks_;
777  MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
778  }
779  }
780 
781 
782 
783  return currentUpdate;
784  }//end getCurrentUpdate() definition
785 
786  template<class ScalarType, class MV, class OP>
788 
789  int i;
790  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
791  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
792 
793  using Teuchos::rcp;
794 
795  // Get correct dimension based on input "dim"
796  // Remember that ortho failures result in an exit before updateLSQR() is called.
797  // Therefore, it is possible that dim == curDim_.
798  int curDim = curDim_;
799  if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
800  curDim = dim;
801  }
802 
804 
805  if(blockSize_ == 1){//if only one right-hand side then use Givens rotations
806  //
807  // Apply previous rotations and compute new rotations to reduce upper-Hessenberg
808  // system to upper-triangular form.
809  //
810  // QR factorization of Least-Squares system with Givens rotations
811  //
812  for (i=0; i<curDim-1; i++) {
813  //
814  // Apply previous Givens rotations to new column of Hessenberg matrix
815  //
816  blas.ROT( 1, &(*R_)(i,curDim-1), 1, &(*R_)(i+1, curDim-1), 1, &cs_[i], &sn_[i] );
817  }
818  //
819  // Calculate new Givens rotation
820  //
821  blas.ROTG( &(*R_)(curDim-1,curDim-1), &(*R_)(curDim,curDim-1), &cs_[curDim-1], &sn_[curDim-1] );
822  (*R_)(curDim,curDim-1) = zero;
823  //
824  // Update RHS w/ new transformation
825  //
826  blas.ROT( 1, &Z_(curDim-1,0), 1, &Z_(curDim,0), 1, &cs_[curDim-1], &sn_[curDim-1] );
827  }
828  else{//if multiple right-hand sides then use Householder transormations
829  //
830  //apply previous reflections and compute new reflections to reduce upper-Hessenberg
831  //system to upper-triagular form.
832 
833  //In Matlab, applying the reflection to a matrix M would look like
834  // M_refl = M - 2*v_refl*(v_refl'*M)/(norm(v_refl)^2)
835 
836  //In order to take advantage of BLAS while applying reflections to a matrix, we
837  //perform it in a two step process
838  //1. workvec = M'*v_refl {using BLAS.GEMV()}
839  //2. M_refl = M_refl - 2*v_refl*workvec'/(norm(v_refl)^2) {using BLAS.GER()}
840 
841  Teuchos::RCP< SDM > workmatrix = Teuchos::null;//matrix of column vectors onto which we apply the reflection
842  Teuchos::RCP< SDV > workvec = Teuchos::null;//where we store the result of the first step of the 2-step reflection process
843  Teuchos::RCP<SDV> v_refl = Teuchos::null;//the reflection vector
844  int R_colStart = curDim_-blockSize_;
845  Teuchos::RCP< SDM >Rblock = Teuchos::null;
846 
847  //
848  //Apply previous reflections
849  //
850  for(i=0; i<lclIter_-1; i++){
851  int R_rowStart = i*blockSize_;
852  //get a view of the part of R_ effected by these reflections.
853  Teuchos::RCP< SDM > RblockCopy = rcp(new SDM (Teuchos::Copy, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
854  Teuchos::RCP< SDM > RblockView = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
855  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());
856 
857  }
858 
859 
860  //Get a view of last 2*blockSize entries of entire block to
861  //generate new reflections.
862  Rblock = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, curDim_-blockSize_, curDim_-blockSize_));
863 
864  //Calculate and apply the new reflections
865  for(i=0; i<blockSize_; i++){
866  //
867  //Calculating Reflection
868  //
869  int curcol = (lclIter_ - 1)*blockSize_ + i;//current column of R_
870  int lclCurcol = i;//current column of Rblock
871  ScalarType signDiag = (*R_)(curcol,curcol) / Teuchos::ScalarTraits<ScalarType>::magnitude((*R_)(curcol,curcol));
872 
873  // Norm of the vector to be reflected.
874  // BLAS returns a ScalarType, but it really should be a magnitude.
875  ScalarType nvs = blas.NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
876  ScalarType alpha = -signDiag*nvs;
877 
878  //norm of reflection vector which is just the vector being reflected
879  //i.e. v = R_(curcol:curcol+blockSize_,curcol))
880  //v_refl = v - alpha*e1
881  //norm(v_refl) = norm(v) + alpha^2 - 2*v*alpha
882  //store in v_refl
883 
884  // Beware, nvs should have a zero imaginary part (since
885  // it is a norm of a vector), but it may not due to rounding
886  // error.
887  //nvs = nvs + alpha*alpha - 2*(*R_)(curcol,curcol)*alpha;
888  //(*R_)(curcol,curcol) -= alpha;
889 
890  //Copy relevant values of the current column of R_ into the reflection vector
891  //Modify first entry
892  //Take norm of reflection vector
893  //Square the norm
894  v_refl = rcp(new SDV(Teuchos::Copy, &((*R_)(curcol,curcol)), blockSize_ + 1 ));
895  (*v_refl)[0] -= alpha;
896  nvs = blas.NRM2(blockSize_+1,v_refl -> values(),1);
897  nvs *= nvs;
898 
899  //
900  //Apply new reflector to:
901  //1. To subsequent columns of R_ in the current block
902  //2. To House[iter_] to store product of reflections for this column
903  //3. To the least-squares right-hand side.
904  //4. The current column
905  //
906  //
907 
908 
909  //
910  //1.
911  //
912  if(i < blockSize_ - 1){//only do this when there are subsquent columns in the block to apply to
913  workvec = Teuchos::rcp(new SDV(blockSize_ - i -1));
914  //workvec = Teuchos::rcp(new SDV(2*blockSize_));
915  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, *Rblock, blockSize_+1, blockSize_ - i -1, lclCurcol, lclCurcol +1 ) );
916  blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix->values(), workmatrix->stride(), v_refl->values(), 1, zero, workvec->values(), 1);
917  blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl->values(),1,workvec->values(),1,workmatrix->values(),workmatrix->stride());
918  }
919 
920 
921  //
922  //2.
923  //
924  workvec = Teuchos::rcp(new SDV(2*blockSize_));
925  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, House_[lclIter_ -1], blockSize_+1, 2*blockSize_, i, 0 ) );
926  blas.GEMV(Teuchos::TRANS,workmatrix->numRows(),workmatrix->numCols(),one,workmatrix->values(),workmatrix->stride(), v_refl->values(), 1,zero,workvec->values(),1);
927  blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl -> values(),1,workvec->values(),1,workmatrix->values(),(*workmatrix).stride());
928 
929  //
930  //3.
931  //
932  workvec = Teuchos::rcp(new SDV(blockSize_));
933  workmatrix = Teuchos::rcp(new SDM (Teuchos::View, Z_, blockSize_+1, blockSize_, curcol, 0 ) );
934  blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix-> values(), workmatrix->stride(), v_refl -> values(), 1, zero, workvec->values(), 1);
935  blas.GER((*workmatrix).numRows(),(*workmatrix).numCols(), -2.0*one/nvs,v_refl -> values(), 1,&((*workvec)[0]),1,(*workmatrix)[0],(*workmatrix).stride());
936 
937  //
938  //4.
939  //
940  (*R_)[curcol][curcol] = alpha;
941  for(int ii=1; ii<= blockSize_; ii++){
942  (*R_)[curcol][curcol+ii] = 0;
943  }
944  }
945 
946  }
947 
948  } // end updateLSQR()
949 
950 
951 }//End Belos Namespace
952 
953 #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:60
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 Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5