Belos  Version of the Day
BelosBlockGCRODRSolMgr.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 
45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP
47 
48 #include "BelosConfigDefs.hpp"
49 #include "BelosTypes.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosSolverManager.hpp"
53 #include "BelosGmresIteration.hpp"
54 #include "BelosBlockGCRODRIter.hpp"
55 #include "BelosBlockGmresIter.hpp"
56 #include "BelosBlockFGmresIter.hpp"
59 #include "BelosStatusTestCombo.hpp"
61 #include "BelosOutputManager.hpp"
62 #include "Teuchos_BLAS.hpp"
63 #include "Teuchos_LAPACK.hpp"
64 #include "Teuchos_as.hpp"
65 #ifdef BELOS_TEUCHOS_TIME_MONITOR
66 #include "Teuchos_TimeMonitor.hpp"
67 #endif // BELOS_TEUCHOS_TIME_MONITOR
68 
69 // MLP: Add links to examples here
70 
71 namespace Belos{
72 
74 
75 
81  public:
82  BlockGCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
83  };
84 
90  public:
91  BlockGCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
92  };
93 
99  public:
100  BlockGCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
101  };
102 
108  public:
109  BlockGCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
110  };
111 
113 
125 
126 template<class ScalarType, class MV, class OP>
127 class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP> {
128 private:
129 
133  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
139 
140 public:
142 
143 
150 
178 
180  virtual ~BlockGCRODRSolMgr() {};
182 
185 
187  std::string description() const;
188 
190 
191 
193 
194 
197  return *problem_;
198  }
199 
202 
205 
207  MagnitudeType achievedTol() const {
208  return achievedTol_;
209  }
210 
212  int getNumIters() const { return numIters_; }
213 
215  bool isLOADetected() const { return loaDetected_; }
216 
218 
220 
221 
224  TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
225  "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
226 
227  // Check state of problem object before proceeding
228  if (! problem->isProblemSet()) {
229  const bool success = problem->setProblem();
230  TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
231  "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the "
232  "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
233  }
234 
235  problem_ = problem;
236  }
237 
240 
242 
244 
245 
252  void reset (const ResetType type) {
253  if ((type & Belos::Problem) && ! problem_.is_null ())
254  problem_->setProblem();
255  else if (type & Belos::RecycleSubspace)
256  keff = 0;
257  }
258 
260 
262 
263 
279  // \returns ::ReturnType specifying:
282  ReturnType solve();
283 
285 
286 private:
287 
288  // Called by all constructors; Contains init instructions common to all constructors
289  void init();
290 
291  // Initialize solver state storage
292  void initializeStateStorage();
293 
294  // Recycling Methods
295  // Appending Function name by:
296  // "Kryl" indicates it is specialized for building a recycle space after an
297  // initial run of Block GMRES which generates an initial unaugmented block Krylov subspace
298  // "AugKryl" indicates it is specialized for building a recycle space from the augmented Krylov subspace
299 
300  // Functions which control the building of a recycle space
301  void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter);
302  void buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
303 
304  // Recycling with Harmonic Ritz Vectors
305  // Computes harmonic eigenpairs of projected matrix created during the priming solve.
306  // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
307 
308  // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension (numBlocks+1)*blockSize x numBlocks.
309  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
310  int getHarmonicVecsKryl (int m, const SDM& HH, SDM& PP);
311 
312  // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+(numBlocks+1)*blockSize x keff+numBlocksm.
313  // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) (numBlocks+1)*blockSize vectors.
314  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
315  int getHarmonicVecsAugKryl (int keff,
316  int m,
317  const SDM& HH,
318  const Teuchos::RCP<const MV>& VV,
319  SDM& PP);
320 
321  // Sort list of n floating-point numbers and return permutation vector
322  void sort (std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm);
323 
324  // Lapack interface
326 
329 
330  //Output Manager
332  Teuchos::RCP<std::ostream> outputStream_;
333 
334  //Status Test
338  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
340 
342  ortho_factory_type orthoFactory_;
343 
350 
353 
364  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
365 
366  //Default Solver Values
367  static const bool adaptiveBlockSize_default_;
368  static const std::string recycleMethod_default_;
369 
370  //Current Solver Values
371  MagnitudeType convTol_, orthoKappa_, achievedTol_;
372  int blockSize_, maxRestarts_, maxIters_, numIters_;
373  int verbosity_, outputStyle_, outputFreq_;
374  bool adaptiveBlockSize_;
375  std::string orthoType_, recycleMethod_;
376  std::string impResScale_, expResScale_;
377  std::string label_;
378 
380  // Solver State Storage
382  //
383  // The number of blocks and recycle blocks (m and k, respectively)
384  int numBlocks_, recycledBlocks_;
385  // Current size of recycled subspace
386  int keff;
387  //
388  // Residual Vector
389  Teuchos::RCP<MV> R_;
390  //
391  // Search Space
392  Teuchos::RCP<MV> V_;
393  //
394  // Recycle subspace and its image under action of the operator
395  Teuchos::RCP<MV> U_, C_;
396  //
397  // Updated recycle Space and its image under action of the operator
398  Teuchos::RCP<MV> U1_, C1_;
399  //
400  // Storage used in constructing recycle space
404  Teuchos::RCP<SDM > PP_;
405  Teuchos::RCP<SDM > HP_;
406  std::vector<ScalarType> tau_;
407  std::vector<ScalarType> work_;
409  std::vector<int> ipiv_;
410 
412  Teuchos::RCP<Teuchos::Time> timerSolve_;
413 
415  bool isSet_;
416 
418  bool loaDetected_;
419 
421  bool builtRecycleSpace_;
422 };
423 
424  //
425  // Set default solver values
426  //
427  template<class ScalarType, class MV, class OP>
428  const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ = true;
429 
430  template<class ScalarType, class MV, class OP>
431  const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ = "harmvecs";
432 
433  //
434  // Method definitions
435  //
436 
437  template<class ScalarType, class MV, class OP>
439  init();
440  }
441 
442  //Basic Constructor
443  template<class ScalarType, class MV, class OP>
447  // Initialize local pointers to null, and initialize local
448  // variables to default values.
449  init ();
450 
451  TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
452  "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
453  "the linear problem argument 'problem' to be nonnull.");
454 
455  problem_ = problem;
456 
457  // Set the parameters using the list that was passed in. If null, we defer initialization until a non-null list is set (by the
458  // client calling setParameters(), or by calling solve() -- in either case, a null parameter list indicates that default parameters should be used).
459  if (! pl.is_null())
460  setParameters (pl);
461  }
462 
463  template<class ScalarType, class MV, class OP>
465  adaptiveBlockSize_ = adaptiveBlockSize_default_;
466  recycleMethod_ = recycleMethod_default_;
467  isSet_ = false;
468  loaDetected_ = false;
469  builtRecycleSpace_ = false;
470  keff = 0;//Effective Size of Recycle Space
471  //The following are all RCP smart pointers to indicated matrices/vectors.
472  //Some MATLAB notation used in comments.
473  R_ = Teuchos::null;//The Block Residual
474  V_ = Teuchos::null;//Block Arnoldi Vectors
475  U_ = Teuchos::null;//Recycle Space
476  C_ = Teuchos::null;//Image of U Under Action of Operator
477  U1_ = Teuchos::null;//Newly Computed Recycle Space
478  C1_ = Teuchos::null;//Image of Newly Computed Recycle Space
479  PP_ = Teuchos::null;//Coordinates of New Recyc. Vectors in Augmented Space
480  HP_ = Teuchos::null;//H_*PP_ or G_*PP_
481  G_ = Teuchos::null;//G_ such that A*[U V(:,1:numBlocks_*blockSize_)] = [C V_]*G_
482  F_ = Teuchos::null;//Upper Triangular portion of QR factorization for HP_
483  H_ = Teuchos::null;//H_ such that A*V(:,1:numBlocks_*blockSize_) = V_*H_ + C_*C_'*V_
484  B_ = Teuchos::null;//B_ = C_'*V_
485 
486  //THIS BLOCK OF CODE IS COMMENTED OUT AND PLACED ELSEWHERE IN THE CODE
487 /*
488  //WE TEMPORARILY INITIALIZE STATUS TESTS HERE FOR TESTING PURPOSES, BUT THEY SHOULD BE
489  //INITIALIZED SOMEWHERE ELSE, LIKE THE setParameters() FUNCTION
490 
491  //INSTANTIATE AND INITIALIZE TEST OBJECTS AS NEEDED
492  if (maxIterTest_.is_null()){
493  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
494  }
495  //maxIterTest_->setMaxIters (maxIters_);
496 
497  //INSTANTIATE THE PRINTER
498  if (printer_.is_null()) {
499  printer_ = Teuchos::rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
500  }
501 
502  if (ortho_.is_null()) // || changedOrthoType) %%%In other codes, this is also triggered if orthogonalization type changed {
503  // Create orthogonalization manager. This requires that the OutputManager (printer_) already be initialized.
504  Teuchos::RCP<const Teuchos::ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
505  ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, orthoParams);
506  }
507 
508  // Convenience typedefs
509  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
510  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
511 
512  if (impConvTest_.is_null()) {
513  impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
514  impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), Belos::TwoNorm);
515  impConvTest_->setShowMaxResNormOnly( true );
516  }
517 
518  if (expConvTest_.is_null()) {
519  expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
520  expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
521  expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), Belos::TwoNorm);
522  expConvTest_->setShowMaxResNormOnly( true );
523  }
524 
525  if (convTest_.is_null()) {
526  convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ, impConvTest_, expConvTest_));
527  }
528 
529  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
530 
531  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
532  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
533  Passed+Failed+Undefined);
534 
535  */
536 
537  }
538 
539  // This method requires the solver manager to return a string that describes itself.
540  template<class ScalarType, class MV, class OP>
542  std::ostringstream oss;
543  oss << "Belos::BlockGCRODRSolMgr<" << SCT::name() << ", ...>";
544  oss << "{";
545  oss << "Ortho Type='"<<orthoType_ ;
546  oss << ", Num Blocks=" <<numBlocks_;
547  oss << ", Block Size=" <<blockSize_;
548  oss << ", Num Recycle Blocks=" << recycledBlocks_;
549  oss << ", Max Restarts=" << maxRestarts_;
550  oss << "}";
551  return oss.str();
552  }
553 
554  template<class ScalarType, class MV, class OP>
558  using Teuchos::parameterList;
559  using Teuchos::RCP;
560  using Teuchos::rcpFromRef;
561 
562  if (defaultParams_.is_null()) {
563  RCP<ParameterList> pl = parameterList ();
564 
565  const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
566  const int maxRestarts = 1000;
567  const int maxIters = 5000;
568  const int blockSize = 2;
569  const int numBlocks = 100;
570  const int numRecycledBlocks = 25;
572  const Belos::OutputType outputStyle = Belos::General;
573  const int outputFreq = 1;
574  RCP<std::ostream> outputStream = rcpFromRef (std::cout);
575  const std::string impResScale ("Norm of Preconditioned Initial Residual");
576  const std::string expResScale ("Norm of Initial Residual");
577  const std::string timerLabel ("Belos");
578  const std::string orthoType ("DGKS");
579  RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
580  //const MagnitudeType orthoKappa = SCT::magnitude (SCT::eps());
581  //
582  // mfh 31 Dec 2011: Negative values of orthoKappa are ignored.
583  // Setting this to a negative value by default ensures that
584  // this parameter is only _not_ ignored if the user
585  // specifically sets a valid value.
586  const MagnitudeType orthoKappa = -SMT::one();
587 
588  // Set all the valid parameters and their default values.
589  pl->set ("Convergence Tolerance", convTol,
590  "The tolerance that the solver needs to achieve in order for "
591  "the linear system(s) to be declared converged. The meaning "
592  "of this tolerance depends on the convergence test details.");
593  pl->set("Maximum Restarts", maxRestarts,
594  "The maximum number of restart cycles (not counting the first) "
595  "allowed for each set of right-hand sides solved.");
596  pl->set("Maximum Iterations", maxIters,
597  "The maximum number of iterations allowed for each set of "
598  "right-hand sides solved.");
599  pl->set("Block Size", blockSize,
600  "How many linear systems to solve at once.");
601  pl->set("Num Blocks", numBlocks,
602  "The maximum number of blocks allowed in the Krylov subspace "
603  "for each set of right-hand sides solved. This includes the "
604  "initial block. Total Krylov basis storage, not counting the "
605  "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
606  pl->set("Num Recycled Blocks", numRecycledBlocks,
607  "The maximum number of vectors in the recycled subspace." );
608  pl->set("Verbosity", verbosity,
609  "What type(s) of solver information should be written "
610  "to the output stream.");
611  pl->set("Output Style", outputStyle,
612  "What style is used for the solver information to write "
613  "to the output stream.");
614  pl->set("Output Frequency", outputFreq,
615  "How often convergence information should be written "
616  "to the output stream.");
617  pl->set("Output Stream", outputStream,
618  "A reference-counted pointer to the output stream where all "
619  "solver output is sent.");
620  pl->set("Implicit Residual Scaling", impResScale,
621  "The type of scaling used in the implicit residual convergence test.");
622  pl->set("Explicit Residual Scaling", expResScale,
623  "The type of scaling used in the explicit residual convergence test.");
624  pl->set("Timer Label", timerLabel,
625  "The string to use as a prefix for the timer labels.");
626  // pl->set("Restart Timers", restartTimers_);
627  pl->set("Orthogonalization", orthoType,
628  "The orthogonalization method to use. Valid options: " +
629  orthoFactory_.validNamesString());
630  pl->set ("Orthogonalization Parameters", *orthoParams,
631  "Sublist of parameters specific to the orthogonalization method to use.");
632  pl->set("Orthogonalization Constant", orthoKappa,
633  "When using DGKS orthogonalization: the \"depTol\" constant, used "
634  "to determine whether another step of classical Gram-Schmidt is "
635  "necessary. Otherwise ignored. Nonpositive values are ignored.");
636  defaultParams_ = pl;
637  }
638  return defaultParams_;
639  }
640 
641  template<class ScalarType, class MV, class OP>
642  void
645  using Teuchos::isParameterType;
646  using Teuchos::getParameter;
647  using Teuchos::null;
649  using Teuchos::parameterList;
650  using Teuchos::RCP;
651  using Teuchos::rcp;
652  using Teuchos::rcp_dynamic_cast;
653  using Teuchos::rcpFromRef;
657 
658  // The default parameter list contains all parameters that BlockGCRODRSolMgr understands, and none that it doesn't
659  RCP<const ParameterList> defaultParams = getValidParameters();
660 
661  if (params.is_null()) {
662  // Users really should supply a non-null input ParameterList,
663  // so that we can fill it in. However, if they really did give
664  // us a null list, be nice and use default parameters. In this
665  // case, we don't have to do any validation.
666  params_ = parameterList (*defaultParams);
667  }
668  else {
669  // Validate the user's parameter list and fill in any missing parameters with default values.
670  //
671  // Currently, validation is quite strict. The following line
672  // will throw an exception for misspelled or extra parameters.
673  // However, it will _not_ throw an exception if there are
674  // missing parameters: these will be filled in with their
675  // default values. There is additional discussion of other
676  // validation strategies in the comments of this function for
677  // Belos::GCRODRSolMgr.
678  params->validateParametersAndSetDefaults (*defaultParams);
679  // No side effects on the solver until after validation passes.
680  params_ = params;
681  }
682 
683  //
684  // Validate and set parameters relating to the Krylov subspace
685  // dimensions and the number of blocks.
686  //
687  // We try to read and validate all these parameters' values
688  // before setting them in the solver. This is because the
689  // validity of each may depend on the others' values. Our goal
690  // is to avoid incorrect side effects, so we don't set the values
691  // in the solver until we know they are right.
692  //
693  {
694  const int maxRestarts = params_->get<int> ("Maximum Restarts");
695  TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
696  "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
697  "must be nonnegative, but you specified a negative value of "
698  << maxRestarts << ".");
699 
700  const int maxIters = params_->get<int> ("Maximum Iterations");
701  TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
702  "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
703  "must be positive, but you specified a nonpositive value of "
704  << maxIters << ".");
705 
706  const int numBlocks = params_->get<int> ("Num Blocks");
707  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
708  "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
709  "positive, but you specified a nonpositive value of " << numBlocks
710  << ".");
711 
712  const int blockSize = params_->get<int> ("Block Size");
713  TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
714  "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
715  "positive, but you specified a nonpositive value of " << blockSize
716  << ".");
717 
718  const int recycledBlocks = params_->get<int> ("Num Recycled Blocks");
719  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
720  "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
721  "be positive, but you specified a nonpositive value of "
722  << recycledBlocks << ".");
723  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
724  std::invalid_argument, "Belos::BlockGCRODRSolMgr: The \"Num Recycled "
725  "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
726  "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
727  << " and \"Num Blocks\" = " << numBlocks << ".");
728 
729  // Now that we've validated the various dimension-related
730  // parameters, we can set them in the solvers.
731  maxRestarts_ = maxRestarts;
732  maxIters_ = maxIters;
733  numBlocks_ = numBlocks;
734  blockSize_ = blockSize;
735  recycledBlocks_ = recycledBlocks;
736  }
737 
738  // Check to see if the timer label changed. If it did, update it in
739  // the parameter list, and create a new timer with that label (if
740  // Belos was compiled with timers enabled).
741  {
742  std::string tempLabel = params_->get<std::string> ("Timer Label");
743  const bool labelChanged = (tempLabel != label_);
744 
745 #ifdef BELOS_TEUCHOS_TIME_MONITOR
746  std::string solveLabel = label_ + ": BlockGCRODRSolMgr total solve time";
747  if (timerSolve_.is_null()) {
748  // We haven't created a timer yet.
749  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
750  } else if (labelChanged) {
751  // We created a timer before with a different label. In that
752  // case, clear the old timer and create a new timer with the
753  // new label. Clearing old timers prevents them from piling
754  // up, since they are stored statically, possibly without
755  // checking for duplicates.
757  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
758  }
759 #endif // BELOS_TEUCHOS_TIME_MONITOR
760  }
761 
762  // Check for a change in verbosity level. Just in case, we allow
763  // the type of this parameter to be either int or MsgType (an
764  // enum).
765  if (params_->isParameter ("Verbosity")) {
766  if (isParameterType<int> (*params_, "Verbosity")) {
767  verbosity_ = params_->get<int> ("Verbosity");
768  }
769  else {
770  verbosity_ = (int) getParameter<MsgType> (*params_, "Verbosity");
771  }
772  }
773 
774  // Check for a change in output style. Just in case, we allow
775  // the type of this parameter to be either int or OutputType (an
776  // enum).
777  if (params_->isParameter ("Output Style")) {
778  if (isParameterType<int> (*params_, "Output Style")) {
779  outputStyle_ = params_->get<int> ("Output Style");
780  }
781  else {
782  outputStyle_ = (int) getParameter<OutputType> (*params_, "Output Style");
783  }
784 
785  // Currently, the only way we can update the output style is to
786  // create a new output test. This is because we must first
787  // instantiate a new StatusTestOutputFactory with the new
788  // style, and then use the factory to make a new output test.
789  // Thus, we set outputTest_ to null for now, so we remember
790  // later to create a new one.
791  outputTest_ = null;
792  }
793 
794  // Get the output stream for the output manager.
795  //
796  // It has been pointed out (mfh 28 Feb 2011 in GCRODRSolMgr code)
797  // that including an RCP<std::ostream> parameter makes it
798  // impossible to serialize the parameter list. If you serialize
799  // the list and read it back in from the serialized
800  // representation, you might not even get a valid
801  // RCP<std::ostream>, let alone the same output stream that you
802  // serialized.
803  //
804  // However, existing Belos users depend on setting the output
805  // stream in the parameter list. We retain this behavior for
806  // backwards compatibility.
807  //
808  // In case the output stream can't be read back in, we default to
809  // stdout (std::cout), just to ensure reasonable behavior.
810  if (params_->isParameter ("Output Stream")) {
811  try {
812  outputStream_ = getParameter<RCP<std::ostream> > (*params_, "Output Stream");
813  }
814  catch (InvalidParameter&) {
815  outputStream_ = rcpFromRef (std::cout);
816  }
817  // We assume that a null output stream indicates that the user
818  // doesn't want to print anything, so we replace it with a
819  // "black hole" stream that prints nothing sent to it. (We
820  // can't use outputStream_ = Teuchos::null, since the output
821  // manager assumes that operator<< works on its output stream.)
822  if (outputStream_.is_null()) {
823  outputStream_ = rcp (new Teuchos::oblackholestream);
824  }
825  }
826 
827  outputFreq_ = params_->get<int> ("Output Frequency");
828  if (verbosity_ & Belos::StatusTestDetails) {
829  // Update parameter in our output status test.
830  if (! outputTest_.is_null()) {
831  outputTest_->setOutputFrequency (outputFreq_);
832  }
833  }
834 
835  // Create output manager if we need to, using the verbosity level
836  // and output stream that we fetched above. We do this here because
837  // instantiating an OrthoManager using OrthoManagerFactory requires
838  // a valid OutputManager.
839  if (printer_.is_null()) {
840  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
841  } else {
842  printer_->setVerbosity (verbosity_);
843  printer_->setOStream (outputStream_);
844  }
845 
846  // Get the orthogonalization manager name ("Orthogonalization").
847  //
848  // Getting default values for the orthogonalization manager
849  // parameters ("Orthogonalization Parameters") requires knowing the
850  // orthogonalization manager name. Save it for later, and also
851  // record whether it's different than before.
852 
853  //bool changedOrthoType = false; // silence warning about not being used
854  if (params_->isParameter ("Orthogonalization")) {
855  const std::string& tempOrthoType =
856  params_->get<std::string> ("Orthogonalization");
857  // Ensure that the specified orthogonalization type is valid.
858  if (! orthoFactory_.isValidName (tempOrthoType)) {
859  std::ostringstream os;
860  os << "Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
861  << tempOrthoType << "\". The following are valid options "
862  << "for the \"Orthogonalization\" name parameter: ";
863  orthoFactory_.printValidNames (os);
864  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
865  }
866  if (tempOrthoType != orthoType_) {
867  //changedOrthoType = true;
868  orthoType_ = tempOrthoType;
869  }
870  }
871 
872  // Get any parameters for the orthogonalization
873  // ("Orthogonalization Parameters"). If not supplied, the
874  // orthogonalization manager factory will supply default values.
875  //
876  // NOTE (mfh 12 Jan 2011, summer 2011, 31 Dec 2011) For backwards
877  // compatibility with other Belos GMRES-type solvers, if params
878  // has an "Orthogonalization Constant" parameter and the DGKS
879  // orthogonalization manager is to be used, the value of this
880  // parameter will override DGKS's "depTol" parameter.
881  //
882  // Users must supply the orthogonalization manager parameters as
883  // a sublist (supplying it as an RCP<ParameterList> would make
884  // the resulting parameter list not serializable).
885  RCP<ParameterList> orthoParams = sublist (params, "Orthogonalization Parameters", true);
886  TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
887  "Failed to get orthogonalization parameters. "
888  "Please report this bug to the Belos developers.");
889 
890  // Check if the desired orthogonalization method changed, or if
891  // the orthogonalization manager has not yet been instantiated.
892  // If either is the case, instantiate a new MatOrthoManager
893  // subclass instance corresponding to the desired
894  // orthogonalization method. We've already fetched the
895  // orthogonalization method name (orthoType_) and its parameters
896  // (orthoParams) above.
897  //
898  // As suggested (by mfh 12 Jan 2011 in Belos::GCRODRSolMgr) In
899  // order to ensure parameter changes get propagated to the
900  // orthomanager we simply reinstantiate the OrthoManager every
901  // time, whether or not the orthogonalization method name or
902  // parameters have changed. This is not efficient for some
903  // orthogonalization managers that keep a lot of internal state.
904  // A more general way to fix this inefficiency would be to
905  //
906  // 1. make each orthogonalization implement
907  // Teuchos::ParameterListAcceptor, and
908  //
909  // 2. make each orthogonalization implement a way to set the
910  // OutputManager instance and timer label.
911  //
912  // Create orthogonalization manager. This requires that the
913  // OutputManager (printer_) already be initialized.
914  ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
915  label_, orthoParams);
916 
917  // The DGKS orthogonalization accepts a "Orthogonalization
918  // Constant" parameter (also called kappa in the code, but not in
919  // the parameter list). If its value is provided in the given
920  // parameter list and its value is positive, use it. Ignore
921  // negative values.
922  //
923  // The "Orthogonalization Constant" parameter is retained only
924  // for backwards compatibility.
925 
926  //bool gotValidOrthoKappa = false; // silence warning about not being used
927  if (params_->isParameter ("Orthogonalization Constant")) {
928  const MagnitudeType orthoKappa =
929  params_->get<MagnitudeType> ("Orthogonalization Constant");
930  if (orthoKappa > 0) {
931  orthoKappa_ = orthoKappa;
932  //gotValidOrthoKappa = true;
933  // Only DGKS currently accepts this parameter.
934  if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
935  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
936  // This cast should always succeed; it's a bug otherwise.
937  // (If the cast fails, then orthoType_ doesn't correspond
938  // to the OrthoManager subclass instance that we think we
939  // have, so we initialized the wrong subclass somehow.)
940  rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
941  }
942  }
943  }
944 
945  // Convergence
946  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
947  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
948 
949  // Check for convergence tolerance
950  convTol_ = params_->get<MagnitudeType> ("Convergence Tolerance");
951  if (! impConvTest_.is_null()) {
952  impConvTest_->setTolerance (convTol_);
953  }
954  if (! expConvTest_.is_null()) {
955  expConvTest_->setTolerance (convTol_);
956  }
957 
958  // Check for a change in scaling, if so we need to build new residual tests.
959  if (params_->isParameter ("Implicit Residual Scaling")) {
960  std::string tempImpResScale =
961  getParameter<std::string> (*params_, "Implicit Residual Scaling");
962 
963  // Only update the scaling if it's different.
964  if (impResScale_ != tempImpResScale) {
965  ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
966  impResScale_ = tempImpResScale;
967 
968  if (! impConvTest_.is_null()) {
969  try {
970  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
971  }
972  catch (StatusTestError&) {
973  // Delete the convergence test so it gets constructed again.
974  impConvTest_ = null;
975  convTest_ = null;
976  }
977  }
978  }
979  }
980 
981  if (params_->isParameter("Explicit Residual Scaling")) {
982  std::string tempExpResScale =
983  getParameter<std::string> (*params_, "Explicit Residual Scaling");
984 
985  // Only update the scaling if it's different.
986  if (expResScale_ != tempExpResScale) {
987  ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
988  expResScale_ = tempExpResScale;
989 
990  if (! expConvTest_.is_null()) {
991  try {
992  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
993  }
994  catch (StatusTestError&) {
995  // Delete the convergence test so it gets constructed again.
996  expConvTest_ = null;
997  convTest_ = null;
998  }
999  }
1000  }
1001  }
1002 
1003  //
1004  // Create iteration stopping criteria ("status tests") if we need
1005  // to, by combining three different stopping criteria.
1006  //
1007  // First, construct maximum-number-of-iterations stopping criterion.
1008  if (maxIterTest_.is_null()) {
1009  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
1010  } else {
1011  maxIterTest_->setMaxIters (maxIters_);
1012  }
1013 
1014  // Implicit residual test, using the native residual to determine if
1015  // convergence was achieved.
1016  if (impConvTest_.is_null()) {
1017  impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1018  impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
1019  Belos::TwoNorm);
1020  }
1021 
1022  // Explicit residual test once the native residual is below the tolerance
1023  if (expConvTest_.is_null()) {
1024  expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1025  expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
1026  expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
1027  Belos::TwoNorm);
1028  }
1029  // Convergence test first tests the implicit residual, then the
1030  // explicit residual if the implicit residual test passes.
1031  if (convTest_.is_null()) {
1032  convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1033  impConvTest_,
1034  expConvTest_));
1035  }
1036  // Construct the complete iteration stopping criterion:
1037  //
1038  // "Stop iterating if the maximum number of iterations has been
1039  // reached, or if the convergence test passes."
1040  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
1041  maxIterTest_, convTest_));
1042  // Create the status test output class.
1043  // This class manages and formats the output from the status test.
1044  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
1045  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
1047 
1048  // Set the solver string for the output test
1049  std::string solverDesc = "Block GCRODR ";
1050  outputTest_->setSolverDesc (solverDesc);
1051 
1052  // Inform the solver manager that the current parameters were set.
1053  isSet_ = true;
1054  }
1055 
1056  // initializeStateStorage.
1057  template<class ScalarType, class MV, class OP>
1058  void
1060  {
1061 
1062  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1063 
1064  // Check if there is any multivector to clone from.
1065  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1066 
1067  //The Dimension of the Krylov Subspace
1068  int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1069 
1070  //Number of columns in [U_ V_(:,1:KrylSpaDim)]
1071  int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;// + 1 is for possible extra recycle vector
1072 
1073  //Number of columns in [C_ V_]
1074  int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1075 
1076  //TEMPORARY SKELETON DEFINITION OF THIS FUNCTION TO GET THINGS WORKING
1077  //NOT EVERYTHING IS INITIALIZE CORRECTLY YET.
1078 
1079  //INITIALIZE RECYCLE SPACE VARIABLES HERE
1080 
1081  //WE DO NOT ALLOCATE V HERE IN THIS SOLVERMANAGER. If no recycle space exists, then block_gmres_iter
1082  //will allocated V for us. If a recycle space already exists, then we will allocate V after updating the
1083  //recycle space for the current problem.
1084  // If the block Krylov subspace has not been initialized before, generate it using the RHS from lp_.
1085  /*if (V_ == Teuchos::null) {
1086  V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1087  }
1088  else{
1089  // Generate V_ by cloning itself ONLY if more space is needed.
1090  if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1091  Teuchos::RCP<const MV> tmp = V_;
1092  V_ = MVT::Clone( *tmp, numBlocks_+1 );
1093  }
1094  }*/
1095 
1096  //INTITIALIZE SPACE FOR THE NEWLY COMPUTED RECYCLE SPACE VARIABLES HERE
1097 
1098  if (U_ == Teuchos::null) {
1099  U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1100  }
1101  else {
1102  // Generate U_ by cloning itself ONLY if more space is needed.
1103  if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1104  Teuchos::RCP<const MV> tmp = U_;
1105  U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1106  }
1107  }
1108 
1109  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1110  if (C_ == Teuchos::null) {
1111  C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1112  }
1113  else {
1114  // Generate C_ by cloning itself ONLY if more space is needed.
1115  if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1116  Teuchos::RCP<const MV> tmp = C_;
1117  C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1118  }
1119  }
1120 
1121  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1122  if (U1_ == Teuchos::null) {
1123  U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1124  }
1125  else {
1126  // Generate U1_ by cloning itself ONLY if more space is needed.
1127  if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1128  Teuchos::RCP<const MV> tmp = U1_;
1129  U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1130  }
1131  }
1132 
1133  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1134  if (C1_ == Teuchos::null) {
1135  C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1136  }
1137  else {
1138  // Generate C1_ by cloning itself ONLY if more space is needed.
1139  if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1140  Teuchos::RCP<const MV> tmp = C1_;
1141  C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1142  }
1143  }
1144 
1145  // Generate R_ only if it doesn't exist
1146  if (R_ == Teuchos::null){
1147  R_ = MVT::Clone( *rhsMV, blockSize_ );
1148  }
1149 
1150  //INITIALIZE SOME WORK VARIABLES
1151 
1152  // Generate G_ only if it doesn't exist, otherwise resize it.
1153  if (G_ == Teuchos::null){
1154  G_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1155  }
1156  else{
1157  if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1158  {
1159  G_->reshape( augSpaImgDim, augSpaDim );
1160  }
1161  G_->putScalar(zero);
1162  }
1163 
1164  // Generate H_ only if it doesn't exist by pointing it to a view of G_.
1165  if (H_ == Teuchos::null){
1166  H_ = Teuchos::rcp (new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1167  }
1168 
1169  // Generate F_ only if it doesn't exist, otherwise resize it.
1170  if (F_ == Teuchos::null){
1171  F_ = Teuchos::rcp( new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1172  }
1173  else {
1174  if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1175  F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1176  }
1177  }
1178  F_->putScalar(zero);
1179 
1180  // Generate PP_ only if it doesn't exist, otherwise resize it.
1181  if (PP_ == Teuchos::null){
1182  PP_ = Teuchos::rcp( new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1183  }
1184  else{
1185  if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1186  PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1187  }
1188  }
1189 
1190  // Generate HP_ only if it doesn't exist, otherwise resize it.
1191  if (HP_ == Teuchos::null)
1192  HP_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1193  else{
1194  if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1195  HP_->reshape( augSpaImgDim, augSpaDim );
1196  }
1197  }
1198 
1199  // Size of tau_ will change during computation, so just be sure it starts with appropriate size
1200  tau_.resize(recycledBlocks_+1);
1201 
1202  // Size of work_ will change during computation, so just be sure it starts with appropriate size
1203  work_.resize(recycledBlocks_+1);
1204 
1205  // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
1206  ipiv_.resize(recycledBlocks_+1);
1207 
1208  }
1209 
1210 template<class ScalarType, class MV, class OP>
1211 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
1212 
1213  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1214  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1215 
1216  int p = block_gmres_iter->getState().curDim; //Dimension of the Krylov space generated
1217  std::vector<int> index(keff);//we use this to index certain columns of U, C, and V to
1218  //get views into pieces of these matrices.
1219 
1220  //GET CORRECT PIECE OF MATRIX H APPROPRIATE TO SIZE OF KRYLOV SUBSPACE
1221  SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1222  if(recycledBlocks_ >= p + blockSize_){//keep whole block Krylov subspace
1223  //IF THIS HAS HAPPENED, THIS MEANS WE CONVERGED DURING THIS CYCLE
1224  //THEREFORE, WE DO NOT CONSTRUCT C = A*U;
1225  index.resize(p);
1226  for (int ii=0; ii < p; ++ii) index[ii] = ii;
1227  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1228  MVT::SetBlock(*V_, index, *Utmp);
1229  keff = p;
1230  }
1231  else{ //use a subspace selection method to get recycle space
1232  int info = 0;
1233  Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1234  if(recycleMethod_ == "harmvecs"){
1235  keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1236  printer_->stream(Debug) << "keff = " << keff << std::endl;
1237  }
1238 // Hereafter, only keff columns of PP are needed
1239 PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, keff ) );
1240 // Now get views into C, U, V
1241 index.resize(keff);
1242 for (int ii=0; ii<keff; ++ii) index[ii] = ii;
1243 Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1244 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1245 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1246 index.resize(p);
1247 for (int ii=0; ii < p; ++ii) index[ii] = ii;
1248 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1249 
1250 // Form U (the subspace to recycle)
1251 // U = newstate.V(:,1:p) * PP;
1252 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1253 
1254 // Step #1: Form HP = H*P
1255 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1256 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1257 // Step #1.5: Perform workspace size query for QR factorization of HP
1258 int lwork = -1;
1259 tau_.resize(keff);
1260 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1261 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1262 
1263 // Step #2: Compute QR factorization of HP
1264  lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1265 work_.resize(lwork);
1266 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1267 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1268 
1269 // Step #3: Explicitly construct Q and R factors
1270 // The upper triangular part of HP is copied into R and HP becomes Q.
1271 SDM Rtmp( Teuchos::View, *F_, keff, keff );
1272 for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1273 //lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1274 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1275 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1276  // Now we have [Q,R] = qr(H*P)
1277 
1278  // Now compute C = V(:,1:p+blockSize_) * Q
1279  index.resize( p+blockSize_ );
1280  for (int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1281  Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+blockSize_ vectors now; needed p above)
1282  MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1283 
1284  // Finally, compute U = U*R^{-1}.
1285  // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
1286  // backsolve capabilities don't exist in the Belos::MultiVec class
1287 
1288 // Step #1: First, compute LU factorization of R
1289 ipiv_.resize(Rtmp.numRows());
1290 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1291 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1292 // Step #2: Form inv(R)
1293 lwork = Rtmp.numRows();
1294 work_.resize(lwork);
1295 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1296 //TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1297 // Step #3: Let U = U * R^{-1}
1298 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1299 
1300 } //end else from if(recycledBlocks_ >= p + 1)
1301 return;
1302 } // end buildRecycleSpaceKryl defnition
1303 
1304 template<class ScalarType, class MV, class OP>
1305 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
1306  const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1307  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1308 
1309  std::vector<MagnitudeType> d(keff);
1310  std::vector<ScalarType> dscalar(keff);
1311  std::vector<int> index(numBlocks_+1);
1312 
1313  // Get the state
1314  BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1315  int p = oldState.curDim;
1316 
1317  // insufficient new information to update recycle space
1318  if (p<1) return;
1319 
1320  if(recycledBlocks_ >= keff + p){ //we add new Krylov vectors to existing recycle space
1321  // If this has happened, we have converged during this cycle. Therefore, do not construct C = A*C;
1322  index.resize(p);
1323  for (int ii=0; ii < p; ++ii) { index[ii] = keff+ii; } //get a view after current reycle vectors
1324  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1325  for (int ii=0; ii < p; ++ii) { index[ii] =ii; }
1326  MVT::SetBlock(*V_, index, *Utmp);
1327  keff += p;
1328  }
1329 
1330  // Take the norm of the recycled vectors
1331  {
1332  index.resize(keff);
1333  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1334  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1335  d.resize(keff);
1336  dscalar.resize(keff);
1337  MVT::MvNorm( *Utmp, d );
1338  for (int i=0; i<keff; ++i) {
1339  d[i] = one / d[i];
1340  dscalar[i] = (ScalarType)d[i];
1341  }
1342  MVT::MvScale( *Utmp, dscalar );
1343  }
1344 
1345  // Get view into current "full" upper Hessnburg matrix
1346  // note that p describes the dimension of the iter+1 block Krylov space so we have to adjust how we use it to index Gtmp
1347  Teuchos::RCP<SDM> Gtmp = Teuchos::rcp( new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1348 
1349  // Insert D into the leading keff x keff block of G
1350  for (int i=0; i<keff; ++i)
1351  (*Gtmp)(i,i) = d[i];
1352 
1353  // Compute the harmoic Ritz pairs for the generalized eigenproblem
1354  // getHarmonicVecsKryl assumes PP has recycledBlocks_+1 columns available
1355  // See previous block of comments for why we subtract p-blockSize_
1356  int keff_new;
1357  {
1358  SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1359  keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1360  }
1361 
1362  // Code to form new U, C
1363  // U = [U V(:,1:p)] * P; (in two steps)
1364 
1365  // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
1366  Teuchos::RCP<MV> U1tmp;
1367  {
1368  index.resize( keff );
1369  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1370  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1371  index.resize( keff_new );
1372  for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1373  U1tmp = MVT::CloneViewNonConst( *U1_, index );
1374  SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1375  MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1376  }
1377 
1378  // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
1379  {
1380  index.resize(p-blockSize_);
1381  for (int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1382  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1383  SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1384  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1385  }
1386 
1387  // Form GP = G*P
1388  SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1389  {
1390  SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1391  HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1392  }
1393 
1394  // Workspace size query for QR factorization HP= QF (the worksize will be placed in work_[0])
1395  int info = 0, lwork = -1;
1396  tau_.resize(keff_new);
1397  lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1398  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1399 
1400  lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1401  work_.resize(lwork);
1402  lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1403  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1404 
1405  // Explicitly construct Q and F factors
1406  // NOTE: The upper triangular part of HP is copied into F and HP becomes Q.
1407  SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1408  for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1409  //lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1410  lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1411  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1412 
1413  // Form orthonormalized C and adjust U accordingly so that C = A*U
1414  // C = [C V] * Q;
1415 
1416  // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
1417  {
1418  Teuchos::RCP<MV> C1tmp;
1419  {
1420  index.resize(keff);
1421  for (int i=0; i < keff; i++) { index[i] = i; }
1422  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1423  index.resize(keff_new);
1424  for (int i=0; i < keff_new; i++) { index[i] = i; }
1425  C1tmp = MVT::CloneViewNonConst( *C1_, index );
1426  SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1427  MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1428  }
1429  // Now compute C += V(:,1:p+1) * Q
1430  {
1431  index.resize( p );
1432  for (int i=0; i < p; ++i) { index[i] = i; }
1433  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1434  SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1435  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1436  }
1437  }
1438 
1439  // C_ = C1_; (via a swap)
1440  std::swap(C_, C1_);
1441 
1442  // Finally, compute U_ = U_*R^{-1}
1443  // First, compute LU factorization of R
1444  ipiv_.resize(Ftmp.numRows());
1445  lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1446  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1447 
1448  // Now, form inv(R)
1449  lwork = Ftmp.numRows();
1450  work_.resize(lwork);
1451  lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1452  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1453 
1454  {
1455  index.resize(keff_new);
1456  for (int i=0; i < keff_new; i++) { index[i] = i; }
1457  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1458  MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1459  }
1460 
1461  // Set the current number of recycled blocks and subspace
1462  // dimension with the Block GCRO-DR iteration.
1463  if (keff != keff_new) {
1464  keff = keff_new;
1465  block_gcrodr_iter->setSize( keff, numBlocks_ );
1466  // Important to zero this out before next cyle
1467  SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1468  b1.putScalar(zero);
1469  }
1470 
1471 } //end buildRecycleSpaceAugKryl definition
1472 
1473 template<class ScalarType, class MV, class OP>
1474 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(int keff, int m, const SDM& GG, const Teuchos::RCP<const MV>& VV, SDM& PP){
1475  int i, j;
1476  int m2 = GG.numCols();
1477  bool xtraVec = false;
1478  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1479  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1480  std::vector<int> index;
1481 
1482  // Real and imaginary eigenvalue components
1483  std::vector<MagnitudeType> wr(m2), wi(m2);
1484 
1485  // Magnitude of harmonic Ritz values
1486  std::vector<MagnitudeType> w(m2);
1487 
1488  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1489  SDM vr(m2,m2,false);
1490 
1491  // Sorted order of harmonic Ritz values
1492  std::vector<int> iperm(m2);
1493 
1494  // Set flag indicating recycle space has been generated this solve
1495  builtRecycleSpace_ = true;
1496 
1497  // Form matrices for generalized eigenproblem
1498 
1499  // B = G' * G; Don't zero out matrix when constructing
1500  SDM B(m2,m2,false);
1501  B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
1502 
1503  // A_tmp = | C'*U 0 |
1504  // | V_{m+1}'*U I |
1505  SDM A_tmp( keff+m+blockSize_, keff+m );
1506 
1507 
1508  // A_tmp(1:keff,1:keff) = C' * U;
1509  index.resize(keff);
1510  for (int i=0; i<keff; ++i) { index[i] = i; }
1511  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1512  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1513  SDM A11( Teuchos::View, A_tmp, keff, keff );
1514  MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1515 
1516  // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
1517  SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1518  index.resize(m+blockSize_);
1519  for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1520  Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1521  MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1522 
1523  // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
1524  for( i=keff; i<keff+m; i++)
1525  A_tmp(i,i) = one;
1526 
1527  // A = G' * A_tmp;
1528  SDM A( m2, A_tmp.numCols() );
1529  A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1530 
1531  // Compute k smallest harmonic Ritz pairs
1532  // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
1533  // ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
1534  // IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
1535  // RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
1536  // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
1537  char balanc='P', jobvl='N', jobvr='V', sense='N';
1538  int ld = A.numRows();
1539  int lwork = 6*ld;
1540  int ldvl = ld, ldvr = ld;
1541  int info = 0,ilo = 0,ihi = 0;
1542  MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1543  ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
1544  std::vector<ScalarType> beta(ld);
1545  std::vector<ScalarType> work(lwork);
1546  std::vector<MagnitudeType> rwork(lwork);
1547  std::vector<MagnitudeType> lscale(ld), rscale(ld);
1548  std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1549  std::vector<int> iwork(ld+6);
1550  int *bwork = 0; // If sense == 'N', bwork is never referenced
1551  //lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1552  // &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1553  // &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
1554  lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1555  &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1556  &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1557  &iwork[0], bwork, &info);
1558  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1559 
1560  // Construct magnitude of each harmonic Ritz value
1561  // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
1562  for( i=0; i<ld; i++ ) // Construct magnitude of each harmonic Ritz value
1563  w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1564 
1565  this->sort(w,ld,iperm);
1566 
1567  bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1568 
1569  // Select recycledBlocks_ smallest eigenvectors
1570  for( i=0; i<recycledBlocks_; i++ )
1571  for( j=0; j<ld; j++ )
1572  PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1573 
1574  if(scalarTypeIsComplex==false) {
1575 
1576  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1577  if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1578  int countImag = 0;
1579  for ( i=ld-recycledBlocks_; i<ld; i++ )
1580  if (wi[iperm[i]] != 0.0) countImag++;
1581  // Check to see if this count is even or odd:
1582  if (countImag % 2) xtraVec = true;
1583  }
1584 
1585  if (xtraVec) { // we need to store one more vector
1586  if (wi[iperm[ld-recycledBlocks_]] > 0.0) { // I picked the "real" component
1587  for( j=0; j<ld; j++ ) // so get the "imag" component
1588  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1589  }
1590  else { // I picked the "imag" component
1591  for( j=0; j<ld; j++ ) // so get the "real" component
1592  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1593  }
1594  }
1595 
1596  }
1597 
1598  // Return whether we needed to store an additional vector
1599  if (xtraVec)
1600  return recycledBlocks_+1;
1601  else
1602  return recycledBlocks_;
1603 
1604 } //end getHarmonicVecsAugKryl definition
1605 
1606 template<class ScalarType, class MV, class OP>
1607 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(int m, const SDM& HH, SDM& PP){
1608  bool xtraVec = false;
1609  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1610  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1611 
1612  // Real and imaginary eigenvalue components
1613  std::vector<MagnitudeType> wr(m), wi(m);
1614 
1615  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1616  SDM vr(m,m,false);
1617 
1618  // Magnitude of harmonic Ritz values
1619  std::vector<MagnitudeType> w(m);
1620 
1621  // Sorted order of harmonic Ritz values, also used for DGEEV
1622  std::vector<int> iperm(m);
1623 
1624  // Size of workspace and workspace for DGEEV
1625  int lwork = 4*m;
1626  std::vector<ScalarType> work(lwork);
1627  std::vector<MagnitudeType> rwork(lwork);
1628 
1629  // Output info
1630  int info = 0;
1631 
1632  // Set flag indicating recycle space has been generated this solve
1633  builtRecycleSpace_ = true;
1634 
1635  // Solve linear system: H_m^{-H}*E_m where E_m is the last blockSize_ columns of the identity matrix
1636  SDM HHt( HH, Teuchos::TRANS );
1637  Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp( new SDM( m, blockSize_));
1638 
1639  //Initialize harmRitzMatrix as E_m
1640  for(int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1641 
1642  //compute harmRitzMatrix <- H_m^{-H}*E_m
1643  lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1644 
1645  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1646  // Compute H_m + H_m^{-H}*E_m*H_lbl^{H}*H_lbl
1647  // H_lbl is bottom-right block of H_, which is a blockSize_ x blockSize_ matrix
1648 
1649  Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1651 
1652  { //So that HTemp will fall out of scope
1653 
1654  // HH_lbl_t <- H_lbl_t*H_lbl
1655  Teuchos::RCP<SDM> Htemp = Teuchos::null;
1656  Htemp = Teuchos::rcp(new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1657  Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1658  H_lbl_t.assign(*Htemp);
1659  // harmRitzMatrix <- harmRitzMatrix*HH_lbl_t
1660  Htemp = Teuchos::rcp(new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1661  Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1662  harmRitzMatrix -> assign(*Htemp);
1663 
1664  // We need to add harmRitzMatrix to the last blockSize_ columns of HH and store the result harmRitzMatrix
1665  int harmColIndex, HHColIndex;
1666  Htemp = Teuchos::rcp(new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1667  for(int i = 0; i<blockSize_; i++){
1668  harmColIndex = harmRitzMatrix -> numCols() - i -1;
1669  HHColIndex = m-i-1;
1670  for(int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1671  }
1672  harmRitzMatrix = Htemp;
1673  }
1674 
1675  // Revise to do query for optimal workspace first
1676  // Create simple storage for the left eigenvectors, which we don't care about.
1677 
1678  const int ldvl = m;
1679  ScalarType* vl = 0;
1680  //lapack.GEEV('N', 'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
1681  // vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
1682  lapack.GEEV('N', 'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
1683  vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1684 
1685  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1686 
1687  // Construct magnitude of each harmonic Ritz value
1688  for( int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1689 
1690  this->sort(w, m, iperm);
1691 
1692  bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1693 
1694  // Select recycledBlocks_ smallest eigenvectors
1695  for( int i=0; i<recycledBlocks_; ++i )
1696  for(int j=0; j<m; j++ )
1697  PP(j,i) = vr(j,iperm[i]);
1698 
1699  if(scalarTypeIsComplex==false) {
1700 
1701  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1702  if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1703  int countImag = 0;
1704  for (int i=0; i<recycledBlocks_; ++i )
1705  if (wi[iperm[i]] != 0.0) countImag++;
1706  // Check to see if this count is even or odd:
1707  if (countImag % 2) xtraVec = true;
1708  }
1709 
1710  if (xtraVec) { // we need to store one more vector
1711  if (wi[iperm[recycledBlocks_-1]] > 0.0) { // I picked the "real" component
1712  for(int j=0; j<m; ++j ) // so get the "imag" component
1713  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1714  }
1715  else{ // I picked the "imag" component
1716  for(int j=0; j<m; ++j ) // so get the "real" component
1717  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1718  }
1719  }
1720 
1721  }
1722 
1723  // Return whether we needed to store an additional vector
1724  if (xtraVec) {
1725  printer_->stream(Debug) << "Recycled " << recycledBlocks_+1 << " vectors" << std::endl;
1726  return recycledBlocks_+1;
1727  }
1728  else {
1729  printer_->stream(Debug) << "Recycled " << recycledBlocks_ << " vectors" << std::endl;
1730  return recycledBlocks_;
1731  }
1732 } //end getHarmonicVecsKryl
1733 
1734 // This method sorts list of n floating-point numbers and return permutation vector
1735 template<class ScalarType, class MV, class OP>
1736 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm) {
1737  int l, r, j, i, flag;
1738  int RR2;
1739  MagnitudeType dRR, dK;
1740 
1741  // Initialize the permutation vector.
1742  for(j=0;j<n;j++)
1743  iperm[j] = j;
1744 
1745  if (n <= 1) return;
1746 
1747  l = n / 2 + 1;
1748  r = n - 1;
1749  l = l - 1;
1750  dRR = dlist[l - 1];
1751  dK = dlist[l - 1];
1752 
1753  RR2 = iperm[l - 1];
1754  while (r != 0) {
1755  j = l;
1756  flag = 1;
1757  while (flag == 1) {
1758  i = j;
1759  j = j + j;
1760  if (j > r + 1)
1761  flag = 0;
1762  else {
1763  if (j < r + 1)
1764  if (dlist[j] > dlist[j - 1]) j = j + 1;
1765  if (dlist[j - 1] > dK) {
1766  dlist[i - 1] = dlist[j - 1];
1767  iperm[i - 1] = iperm[j - 1];
1768  }
1769  else {
1770  flag = 0;
1771  }
1772  }
1773  }
1774  dlist[i - 1] = dRR;
1775  iperm[i - 1] = RR2;
1776  if (l == 1) {
1777  dRR = dlist [r];
1778  RR2 = iperm[r];
1779  dK = dlist[r];
1780  dlist[r] = dlist[0];
1781  iperm[r] = iperm[0];
1782  r = r - 1;
1783  }
1784  else {
1785  l = l - 1;
1786  dRR = dlist[l - 1];
1787  RR2 = iperm[l - 1];
1788  dK = dlist[l - 1];
1789  }
1790  }
1791  dlist[0] = dRR;
1792  iperm[0] = RR2;
1793 } //end sort() definition
1794 
1795 template<class ScalarType, class MV, class OP>
1797  using Teuchos::RCP;
1798  using Teuchos::rcp;
1799  using Teuchos::rcp_const_cast;
1800 
1801  // MLP: NEED TO ADD CHECK IF PARAMETERS ARE SET LATER
1802 
1803  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1804  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1805  std::vector<int> index(numBlocks_+1);
1806 
1807  // MLP: EXCEPTION TESTING SHOULD GO HERE
1808 
1809  //THE FOLLOWING BLOCK OF CODE IS TO INFORM THE PROBLEM HOW MANY RIGHT HAND
1810  //SIDES WILL BE SOLVED. IF THE CHOSEN BLOCK SIZE IS THE SAME AS THE NUMBER
1811  //OF RIGHT HAND SIDES IN THE PROBLEM THEN WE PASS THE PROBLEM
1812  //INDICES FOR ALL RIGHT HAND SIDES. IF BLOCK SIZE IS GREATER THAN
1813  //THE NUMBER OF RIGHT HAND SIDES AND THE USER HAS ADAPTIVE BLOCK SIZING
1814  //TURNED OFF, THEN THE PROBLEM WILL GENERATE RANDOM RIGHT HAND SIDES SO WE HAVE AS MANY
1815  //RIGHT HAND SIDES AS BLOCK SIZE INDICATES. THIS MAY NEED TO BE CHANGED
1816  //LATER TO ACCOMADATE SMARTER CHOOSING OF FICTITIOUS RIGHT HAND SIDES.
1817  //THIS CODE IS DIRECTLY LIFTED FROM BelosBlockGmresSolMgr.hpp.
1818 
1819  // Create indices for the linear systems to be solved.
1820  int startPtr = 0;
1821  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1822  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1823 
1824  //currIdx holds indices to all right-hand sides to be solved
1825  //or -1's to indicate that random right-hand sides should be generated
1826  std::vector<int> currIdx;
1827 
1828  if ( adaptiveBlockSize_ ) {
1829  blockSize_ = numCurrRHS;
1830  currIdx.resize( numCurrRHS );
1831  for (int i=0; i<numCurrRHS; ++i)
1832  currIdx[i] = startPtr+i;
1833  }
1834  else {
1835  currIdx.resize( blockSize_ );
1836  for (int i=0; i<numCurrRHS; ++i)
1837  currIdx[i] = startPtr+i;
1838  for (int i=numCurrRHS; i<blockSize_; ++i)
1839  currIdx[i] = -1;
1840  }
1841 
1842  // Inform the linear problem of the linear systems to solve/generate.
1843  problem_->setLSIndex( currIdx );
1844 
1845  //ADD ERROR CHECKING TO MAKE SURE SIZE OF BLOCK KRYLOV SUBSPACE NOT LARGER THAN dim
1846  //ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1847 
1848  // reset loss of orthogonality flag
1849  loaDetected_ = false;
1850 
1851  // Assume convergence is achieved, then let any failed convergence set this to false.
1852  bool isConverged = true;
1853 
1854  // Initialize storage for all state variables
1855  initializeStateStorage();
1856 
1857  // Parameter list MLP: WE WILL NEED TO ASSIGN SOME PARAMETER VALUES LATER
1858  Teuchos::ParameterList plist;
1859 
1860  while (numRHS2Solve > 0){ //This loops through each block of RHS's being solved
1861  // **************************************************************************************************************************
1862  // Begin initial solver preparations. Either update U,C for new operator or generate an initial space using a cycle of GMRES
1863  // **************************************************************************************************************************
1864  //int prime_iterations; // silence warning about not being used
1865  if(keff > 0){//If there is already a subspace to recycle, then use it
1866  // Update U, C for the new operator
1867 
1868 // TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,BlockGCRODRSolMgrRecyclingFailure,
1869 // "Belos::BlockGCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1870  printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
1871 
1872  // Compute image of U_ under the new operator
1873  index.resize(keff);
1874  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1875  RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1876  RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1877  problem_->apply( *Utmp, *Ctmp );
1878 
1879  RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1880 
1881  // Orthogonalize this block
1882  // Get a matrix to hold the orthonormalization coefficients.
1883  SDM Ftmp( Teuchos::View, *F_, keff, keff );
1884  int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,false));
1885  // Throw an error if we could not orthogonalize this block
1886  TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,BlockGCRODRSolMgrOrthoFailure,"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1887 
1888  // U_ = U_*F^{-1}
1889  // First, compute LU factorization of R
1890  int info = 0;
1891  ipiv_.resize(Ftmp.numRows());
1892  lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1893  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1894 
1895  // Now, form inv(F)
1896  int lwork = Ftmp.numRows();
1897  work_.resize(lwork);
1898  lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1899  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1900 
1901  // U_ = U1_; (via a swap)
1902  MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1903  std::swap(U_, U1_);
1904 
1905  // Must reinitialize after swap
1906  index.resize(keff);
1907  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1908  Ctmp = MVT::CloneViewNonConst( *C_, index );
1909  Utmp = MVT::CloneView( *U_, index );
1910 
1911  // Compute C_'*R_
1912  SDM Ctr(keff,blockSize_);
1913  problem_->computeCurrPrecResVec( &*R_ );
1914  MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1915 
1916  // Update solution ( x += U_*C_'*R_ )
1917  RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1918  MVT::MvInit( *update, 0.0 );
1919  MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1920  problem_->updateSolution( update, true );
1921 
1922  // Update residual norm ( r -= C_*C_'*R_ )
1923  MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1924 
1925  // We recycled space from previous call
1926  //prime_iterations = 0;
1927 
1928  // Since we have a previous recycle space, we do not need block_gmres_iter
1929  // and we must allocate V ourselves. See comments in initialize() routine for
1930  // further explanation.
1931  if (V_ == Teuchos::null) {
1932  // Check if there is a multivector to clone from.
1933  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1934  V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1935  }
1936  else{
1937  // Generate V_ by cloning itself ONLY if more space is needed.
1938  if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1939  Teuchos::RCP<const MV> tmp = V_;
1940  V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1941  }
1942  }
1943  } //end if(keff > 0)
1944  else{ //if there was no subspace to recycle, then build one
1945  printer_->stream(Debug) << " No recycled subspace available for RHS index " << std::endl << std::endl;
1946 
1947  Teuchos::ParameterList primeList;
1948  //we set this as numBlocks-1 because that is the Krylov dimension we want
1949  //so that the size of the Hessenberg created matches that of what we were expecting.
1950  primeList.set("Num Blocks",numBlocks_-1);
1951  primeList.set("Block Size",blockSize_);
1952  primeList.set("Recycled Blocks",0);
1953  primeList.set("Keep Hessenberg",true);
1954  primeList.set("Initialize Hessenberg",true);
1955 
1956  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1957  if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {//if user has selected a total subspace dimension larger than system dimension
1958  ptrdiff_t tmpNumBlocks = 0;
1959  if (blockSize_ == 1)
1960  tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
1961  else{
1962  tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
1963  printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
1964  << std::endl << "The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1965  primeList.set("Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1966  }
1967  }
1968  else{
1969  primeList.set("Num Blocks",numBlocks_-1);
1970  }
1971  //Create Block GMRES iteration object to perform one cycle of GMRES
1973  block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
1974 
1975  // MLP: ADD LOGIC TO DEAL WITH USER ASKING TO GENERATE A LARGER SPACE THAN dim AS IN HEIDI'S BlockGmresSolMgr CODE (DIDN'T WE ALREADY DO THIS SOMEWHERE?)
1976  block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1977 
1978  //Create the first block in the current BLOCK Krylov basis (using residual)
1979  Teuchos::RCP<MV> V_0;
1980  if (currIdx[blockSize_-1] == -1) {//if augmented with random RHS's
1981  V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1982  problem_->computeCurrPrecResVec( &*V_0 );
1983  }
1984  else {
1985  V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1986  }
1987 
1988  // Get a matrix to hold the orthonormalization coefficients.
1989  Teuchos::RCP<SDM > z_0 = Teuchos::rcp( new SDM( blockSize_, blockSize_ ) );
1990 
1991  // Orthonormalize the new V_0
1992  int rank = ortho_->normalize( *V_0, z_0 );
1993  (void) rank; // silence compiler warning for unused variable
1994  // MLP: ADD EXCEPTION IF INITIAL BLOCK IS RANK DEFFICIENT
1995 
1996  // Set the new state and initialize the iteration.
1998  newstate.V = V_0;
1999  newstate.z = z_0;
2000  newstate.curDim = 0;
2001  block_gmres_iter->initializeGmres(newstate);
2002 
2003  bool primeConverged = false;
2004 
2005  try {
2006  printer_->stream(Debug) << " Preparing to Iterate!!!!" << std::endl << std::endl;
2007  block_gmres_iter->iterate();
2008 
2009  // *********************
2010  // check for convergence
2011  // *********************
2012  if ( convTest_->getStatus() == Passed ) {
2013  printer_->stream(Debug) << "We converged during the prime the pump stage" << std::endl << std::endl;
2014  primeConverged = !(expConvTest_->getLOADetected());
2015  if ( expConvTest_->getLOADetected() ) {
2016  // we don't have convergence
2017  loaDetected_ = true;
2018  printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2019  }
2020  }
2021  // *******************************************
2022  // check for maximum iterations of block GMRES
2023  // *******************************************
2024  else if( maxIterTest_->getStatus() == Passed ) {
2025  // we don't have convergence
2026  primeConverged = false;
2027  }
2028  // ****************************************************************
2029  // We need to recycle and continue, print a message indicating this
2030  // ****************************************************************
2031  else{ //debug statements
2032  printer_->stream(Debug) << " We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2033  }
2034  } //end try
2035  catch (const GmresIterationOrthoFailure &e) {
2036  // If the block size is not one, it's not considered a lucky breakdown.
2037  if (blockSize_ != 1) {
2038  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2039  << block_gmres_iter->getNumIters() << std::endl
2040  << e.what() << std::endl;
2041  if (convTest_->getStatus() != Passed)
2042  primeConverged = false;
2043  }
2044  else {
2045  // If the block size is one, try to recover the most recent least-squares solution
2046  block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2047  // Check to see if the most recent least-squares solution yielded convergence.
2048  sTest_->checkStatus( &*block_gmres_iter );
2049  if (convTest_->getStatus() != Passed)
2050  isConverged = false;
2051  }
2052  } // end catch (const GmresIterationOrthoFailure &e)
2053  catch (const std::exception &e) {
2054  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2055  << block_gmres_iter->getNumIters() << std::endl
2056  << e.what() << std::endl;
2057  throw;
2058  }
2059 
2060  // Record number of iterations in generating initial recycle spacec
2061  //prime_iterations = block_gmres_iter->getNumIters();//instantiated here because it is not needed outside of else{} scope; we'll see if this is true or not
2062 
2063  // Update the linear problem.
2064  RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2065  problem_->updateSolution( update, true );
2066 
2067  // Update the block residual
2068 
2069  problem_->computeCurrPrecResVec( &*R_ );
2070 
2071  // Get the state.
2072  newstate = block_gmres_iter->getState();
2073  int p = newstate.curDim;
2074  (void) p; // silence compiler warning for unused variable
2075  H_->assign(*(newstate.H));//IS THIS A GOOD IDEA? I DID IT SO MEMBERS WOULD HAVE ACCESS TO H.
2076  // Get new Krylov vectors, store in V_
2077 
2078  // Since the block_gmres_iter returns the state
2079  // with const RCP's we need to cast newstate.V as
2080  // a non const RCP. This is okay because block_gmres_iter
2081  // is about to fall out of scope, so we are simply
2082  // rescuing *newstate.V from being destroyed so we can
2083  // use it for future block recycled GMRES cycles
2084  V_ = rcp_const_cast<MV>(newstate.V);
2085  newstate.V.release();
2086  //COMPUTE NEW RECYCLE SPACE SOMEHOW
2087  buildRecycleSpaceKryl(keff, block_gmres_iter);
2088  printer_->stream(Debug) << "Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2089 
2090  // Return to outer loop if the priming solve
2091  // converged, set the next linear system.
2092  if (primeConverged) {
2093  /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2094  // Inform the linear problem that we are
2095  // finished with this block linear system.
2096  problem_->setCurrLS();
2097 
2098  // Update indices for the linear systems to be solved.
2099  // KMS: Fix the numRHS2Solve; Copy from Heidi's BlockGmres code
2100  numRHS2Solve -= 1;
2101  printer_->stream(Debug) << numRHS2Solve << " Right hand sides left to solve" << std::endl;
2102  if ( numRHS2Solve > 0 ) {
2103  currIdx[0]++;
2104  // Set the next indices
2105  problem_->setLSIndex( currIdx );
2106  }
2107  else {
2108  currIdx.resize( numRHS2Solve );
2109  }
2110  *****************************************************************************/
2111 
2112  // Inform the linear problem that we are finished with this block linear system.
2113  problem_->setCurrLS();
2114 
2115  // Update indices for the linear systems to be solved.
2116  startPtr += numCurrRHS;
2117  numRHS2Solve -= numCurrRHS;
2118  if ( numRHS2Solve > 0 ) {
2119  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2120  if ( adaptiveBlockSize_ ) {
2121  blockSize_ = numCurrRHS;
2122  currIdx.resize( numCurrRHS );
2123  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2124  }
2125  else {
2126  currIdx.resize( blockSize_ );
2127  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2128  for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2129  }
2130  // Set the next indices.
2131  problem_->setLSIndex( currIdx );
2132  }
2133  else {
2134  currIdx.resize( numRHS2Solve );
2135  }
2136  continue;//Begin solving the next block of RHS's
2137  } //end if (primeConverged)
2138  } //end else [if(keff > 0)]
2139 
2140  // *****************************************
2141  // Initial subspace update/construction done
2142  // Start cycles of recycled GMRES
2143  // *****************************************
2144  Teuchos::ParameterList blockgcrodrList;
2145  blockgcrodrList.set("Num Blocks",numBlocks_);
2146  blockgcrodrList.set("Block Size",blockSize_);
2147  blockgcrodrList.set("Recycled Blocks",keff);
2148 
2150  block_gcrodr_iter = Teuchos::rcp( new BlockGCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,blockgcrodrList) );
2152 
2153  index.resize( blockSize_ );
2154  for(int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2155  Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2156 
2157  // MLP: MVT::SetBlock(*R_,index,*V0);
2158  MVT::Assign(*R_,*V0);
2159 
2160  index.resize(keff);//resize to get appropriate recycle space vectors
2161  for(int i=0; i < keff; i++){ index[i] = i;};
2162  B_ = rcp(new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2163  H_ = rcp(new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2164 
2165  newstate.V = V_;
2166  newstate.B= B_;
2167  newstate.U = MVT::CloneViewNonConst(*U_, index);
2168  newstate.C = MVT::CloneViewNonConst(*C_, index);
2169  newstate.H = H_;
2170  newstate.curDim = blockSize_;
2171  block_gcrodr_iter -> initialize(newstate);
2172 
2173  int numRestarts = 0;
2174 
2175  while(1){//Each execution of this loop is a cycle of block GCRODR
2176  try{
2177  block_gcrodr_iter -> iterate();
2178 
2179  // ***********************
2180  // Check Convergence First
2181  // ***********************
2182  if( convTest_->getStatus() == Passed ) {
2183  // we have convergence
2184  break;//from while(1)
2185  } //end if converged
2186 
2187  // ***********************************
2188  // Check if maximum iterations reached
2189  // ***********************************
2190  else if(maxIterTest_->getStatus() == Passed ){
2191  // no convergence; hit maxit
2192  isConverged = false;
2193  break; // from while(1)
2194  } //end elseif reached maxit
2195 
2196  // **********************************************
2197  // Check if subspace full; do we need to restart?
2198  // **********************************************
2199  else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2200 
2201  // Update recycled space even if we have reached max number of restarts
2202 
2203  // Update linear problem
2204  Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2205  problem_->updateSolution(update, true);
2206  buildRecycleSpaceAugKryl(block_gcrodr_iter);
2207 
2208  printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2209  // NOTE: If we have hit the maximum number of restarts, then we will quit
2210  if(numRestarts >= maxRestarts_) {
2211  isConverged = false;
2212  break; //from while(1)
2213  } //end if max restarts
2214 
2215  numRestarts++;
2216 
2217  printer_ -> stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
2218 
2219  // Create the restart vector (first block in the current Krylov basis)
2220  problem_->computeCurrPrecResVec( &*R_ );
2221  index.resize( blockSize_ );
2222  for (int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2223  Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2224  MVT::SetBlock(*R_,index,*V0);
2225 
2226  // Set the new state and initialize the solver.
2228  index.resize( numBlocks_*blockSize_ );
2229  for (int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2230  restartState.V = MVT::CloneViewNonConst( *V_, index );
2231  index.resize( keff );
2232  for (int ii=0; ii<keff; ++ii) index[ii] = ii;
2233  restartState.U = MVT::CloneViewNonConst( *U_, index );
2234  restartState.C = MVT::CloneViewNonConst( *C_, index );
2235  B_ = rcp(new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2236  H_ = rcp(new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2237  restartState.B = B_;
2238  restartState.H = H_;
2239  restartState.curDim = blockSize_;
2240  block_gcrodr_iter->initialize(restartState);
2241 
2242  } //end else if need to restart
2243 
2244  // ****************************************************************
2245  // We returned from iterate(), but none of our status tests passed.
2246  // Something is wrong, and it is probably our fault.
2247  // ****************************************************************
2248  else {
2249  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2250  } //end else (no status test passed)
2251 
2252  }//end try
2253  catch(const BlockGCRODRIterOrthoFailure &e){ //there was an exception
2254  // Try to recover the most recent least-squares solution
2255  block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2256  // Check to see if the most recent least-squares solution yielded convergence.
2257  sTest_->checkStatus( &*block_gcrodr_iter );
2258  if (convTest_->getStatus() != Passed) isConverged = false;
2259  break;
2260  } // end catch orthogonalization failure
2261  catch(const std::exception &e){
2262  printer_->stream(Errors) << "Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
2263  << block_gcrodr_iter->getNumIters() << std::endl
2264  << e.what() << std::endl;
2265  throw;
2266  } //end catch standard exception
2267  } //end while(1)
2268 
2269  // Compute the current solution.
2270  // Update the linear problem.
2271  Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2272  problem_->updateSolution( update, true );
2273 
2274  /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2275  // Inform the linear problem that we are finished with this block linear system.
2276  problem_->setCurrLS();
2277 
2278  // Update indices for the linear systems to be solved.
2279  numRHS2Solve -= 1;
2280  if ( numRHS2Solve > 0 ) {
2281  currIdx[0]++;
2282  // Set the next indices.
2283  problem_->setLSIndex( currIdx );
2284  }
2285  else {
2286  currIdx.resize( numRHS2Solve );
2287  }
2288  ******************************************************************************/
2289 
2290  // Inform the linear problem that we are finished with this block linear system.
2291  problem_->setCurrLS();
2292 
2293  // Update indices for the linear systems to be solved.
2294  startPtr += numCurrRHS;
2295  numRHS2Solve -= numCurrRHS;
2296  if ( numRHS2Solve > 0 ) {
2297  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2298  if ( adaptiveBlockSize_ ) {
2299  blockSize_ = numCurrRHS;
2300  currIdx.resize( numCurrRHS );
2301  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2302  }
2303  else {
2304  currIdx.resize( blockSize_ );
2305  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2306  for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2307  }
2308  // Set the next indices.
2309  problem_->setLSIndex( currIdx );
2310  }
2311  else {
2312  currIdx.resize( numRHS2Solve );
2313  }
2314 
2315  // If we didn't build a recycle space this solve but ran at least k iterations, force build of new recycle space
2316  if (!builtRecycleSpace_) {
2317  buildRecycleSpaceAugKryl(block_gcrodr_iter);
2318  printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2319  }
2320  }//end while (numRHS2Solve > 0)
2321 
2322  // print final summary
2323  sTest_->print( printer_->stream(FinalSummary) );
2324 
2325  // print timing information
2326  #ifdef BELOS_TEUCHOS_TIME_MONITOR
2327  // Calling summarize() can be expensive, so don't call unless the user wants to print out timing details.
2328  // summarize() will do all the work even if it's passed a "black hole" output stream.
2329  if (verbosity_ & TimingDetails) Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
2330  #endif
2331  // get iteration information for this solve
2332  numIters_ = maxIterTest_->getNumIters();
2333 
2334  // get residual information for this solve
2335  const std::vector<MagnitudeType>* pTestValues = NULL;
2336  pTestValues = impConvTest_->getTestValue();
2337  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2338  "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2339  "getTestValue() method returned NULL. Please report this bug to the "
2340  "Belos developers.");
2341  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2342  "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2343  "getTestValue() method returned a vector of length zero. Please report "
2344  "this bug to the Belos developers.");
2345  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
2346  // achieved tolerances for all vectors in the current solve(), or
2347  // just for the vectors from the last deflation?
2348  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
2349 
2350  if (!isConverged) return Unconverged; // return from BlockGCRODRSolMgr::solve()
2351  return Converged; // return from BlockGCRODRSolMgr::solve()
2352 } //end solve()
2353 
2354 } //End Belos Namespace
2355 
2356 #endif /* BELOS_BLOCK_GCRODR_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:87
Collection of types and exceptions used within the Belos solvers.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
ScalarType * values() const
virtual ~BlockGCRODRSolMgr()
Destructor.
static void clearCounter(const std::string &name)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
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)
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
Thrown when the linear problem was not set up correctly.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
static magnitudeType real(T a)
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set the parameters the solver should use to solve the linear problem.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Structure to contain pointers to GmresIteration state variables.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
Thrown if any problem occurs in using or creating the recycle subspace.
A factory class for generating StatusTestOutput objects.
bool is_null() const
Teuchos::RCP< MV > V
The current Krylov basis.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
int curDim
The current dimension of the reduction.
Traits class which defines basic operations on multivectors.
std::string description() const
A description of the Block GCRODR solver manager.
Belos::StatusTest for logically combining several status tests.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
OrdinalType numRows() const
int curDim
The current dimension of the reduction.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
BlockGCRODRSolMgr()
Default constructor.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Thrown when the solution manager was unable to orthogonalize the basis vectors.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Belos concrete class for performing the block GMRES iteration.
Ptr< T > release()
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
static magnitudeType magnitude(T a)
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
Class which defines basic traits for the operator type.
OrdinalType stride() const
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
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 ...
OrdinalType numCols() const
Belos header file which uses auto-configuration information to include necessary C++ headers...
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
ReturnType solve()
Solve the current linear problem.
Thrown when an LAPACK call fails.
Belos concrete class for performing the block, flexible GMRES iteration.
OutputType
Style of output used to display status test information.
Definition: BelosTypes.hpp:140
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.

Generated for Belos by doxygen 1.8.14