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

Generated on Fri Jun 5 2020 10:20:37 for Belos by doxygen 1.8.5