Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosGCRODRSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_GCRODR_SOLMGR_HPP
11 #define BELOS_GCRODR_SOLMGR_HPP
12 
16 
17 #include "BelosConfigDefs.hpp"
19 #include "BelosSolverManager.hpp"
20 #include "BelosLinearProblem.hpp"
21 #include "BelosTypes.hpp"
22 
23 #include "BelosGCRODRIter.hpp"
24 #include "BelosBlockFGmresIter.hpp"
27 #include "BelosStatusTestCombo.hpp"
29 #include "BelosOutputManager.hpp"
30 #include "Teuchos_BLAS.hpp" // includes Teuchos_ConfigDefs.hpp
31 #include "Teuchos_LAPACK.hpp"
32 #include "Teuchos_as.hpp"
33 
34 #ifdef BELOS_TEUCHOS_TIME_MONITOR
35 # include "Teuchos_TimeMonitor.hpp"
36 #endif // BELOS_TEUCHOS_TIME_MONITOR
37 
38 #if defined(HAVE_TEUCHOSCORE_CXX11)
39 # include <type_traits>
40 # if defined(HAVE_TEUCHOS_COMPLEX)
41 #include "Kokkos_Complex.hpp"
42 # endif
43 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
44 
55 namespace Belos {
56 
58 
59 
67  public:
68  GCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
69  };
70 
78  public:
79  GCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
80  };
81 
89  public:
90  GCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
91  };
92 
100  public:
101  GCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
102  };
103 
105 
130  template<class ScalarType, class MV, class OP,
131  const bool lapackSupportsScalarType =
133  class GCRODRSolMgr :
134  public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
135  {
136  static const bool requiresLapack =
138  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
139  requiresLapack> base_type;
140 
141  public:
143  base_type ()
144  {}
147  base_type ()
148  {}
149  virtual ~GCRODRSolMgr () {}
150  };
151 
155  template<class ScalarType, class MV, class OP>
156  class GCRODRSolMgr<ScalarType, MV, OP, true> :
157  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
158  {
159 
160 #if defined(HAVE_TEUCHOSCORE_CXX11)
161 # if defined(HAVE_TEUCHOS_COMPLEX)
162  #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
163  static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
164  std::is_same<ScalarType, std::complex<double> >::value ||
165  std::is_same<ScalarType, Kokkos::complex<double> >::value ||
166  std::is_same<ScalarType, float>::value ||
167  std::is_same<ScalarType, double>::value ||
168  std::is_same<ScalarType, long double>::value,
169  "Belos::GCRODRSolMgr: ScalarType must be one of the four "
170  "types (S,D,C,Z) supported by LAPACK or long double (largely not impl'd).");
171  #else
172  static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
173  std::is_same<ScalarType, std::complex<double> >::value ||
174  std::is_same<ScalarType, Kokkos::complex<double> >::value ||
175  std::is_same<ScalarType, float>::value ||
176  std::is_same<ScalarType, double>::value,
177  "Belos::GCRODRSolMgr: ScalarType must be one of the four "
178  "types (S,D,C,Z) supported by LAPACK.");
179  #endif
180 # else
181  #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
182  static_assert (std::is_same<ScalarType, float>::value ||
183  std::is_same<ScalarType, double>::value ||
184  std::is_same<ScalarType, long double>::value,
185  "Belos::GCRODRSolMgr: ScalarType must be float, double or long double. "
186  "Complex arithmetic support is currently disabled. To "
187  "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
188  #else
189  static_assert (std::is_same<ScalarType, float>::value ||
190  std::is_same<ScalarType, double>::value,
191  "Belos::GCRODRSolMgr: ScalarType must be float or double. "
192  "Complex arithmetic support is currently disabled. To "
193  "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
194  #endif
195 # endif // defined(HAVE_TEUCHOS_COMPLEX)
196 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
197 
198  private:
202  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
205 
206  public:
208 
209 
215  GCRODRSolMgr();
216 
271 
273  virtual ~GCRODRSolMgr() {};
274 
278  }
280 
282 
283 
286  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
287  return *problem_;
288  }
289 
292  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
293 
297  return params_;
298  }
299 
306  return Teuchos::tuple(timerSolve_);
307  }
308 
314  MagnitudeType achievedTol() const override {
315  return achievedTol_;
316  }
317 
319  int getNumIters() const override {
320  return numIters_;
321  }
322 
325  bool isLOADetected() const override { return false; }
326 
328 
330 
331 
333  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override {
334  problem_ = problem;
335  }
336 
338  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
339 
341 
343 
344 
348  void reset( const ResetType type ) override {
349  if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) {
350  bool set = problem_->setProblem();
351  if (!set)
352  throw "Could not set problem.";
353  }
354  else if (type & Belos::RecycleSubspace) {
355  keff = 0;
356  }
357  }
359 
361 
362 
389  ReturnType solve() override;
390 
392 
394 
396  std::string description() const override;
397 
399 
400  private:
401 
402  // Called by all constructors; Contains init instructions common to all constructors
403  void init();
404 
405  // Initialize solver state storage
406  void initializeStateStorage();
407 
408  // Compute updated recycle space given existing recycle space and newly generated Krylov space
409  void buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
410 
411  // Computes harmonic eigenpairs of projected matrix created during the priming solve.
412  // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
413  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
414  // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
415  int getHarmonicVecs1(int m,
418 
419  // Computes harmonic eigenpairs of projected matrix created during one cycle.
420  // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
421  // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
422  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
423  // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
424  int getHarmonicVecs2(int keff, int m,
426  const Teuchos::RCP<const MV>& VV,
428 
429  // Sort list of n floating-point numbers and return permutation vector
430  void sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm);
431 
432  // Lapack interface
434 
435  // Linear problem.
437 
438  // Output manager.
440  Teuchos::RCP<std::ostream> outputStream_;
441 
442  // Status test.
446  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
448 
453 
454  // Current parameter list.
456 
457  // Default solver values.
458  static constexpr double orthoKappa_default_ = 0.0;
459  static constexpr int maxRestarts_default_ = 100;
460  static constexpr int maxIters_default_ = 1000;
461  static constexpr int numBlocks_default_ = 50;
462  static constexpr int blockSize_default_ = 1;
463  static constexpr int recycledBlocks_default_ = 5;
464  static constexpr int verbosity_default_ = Belos::Errors;
465  static constexpr int outputStyle_default_ = Belos::General;
466  static constexpr int outputFreq_default_ = -1;
467  static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
468  static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
469  static constexpr const char * label_default_ = "Belos";
470  static constexpr const char * orthoType_default_ = "ICGS";
471 
472  // Current solver values.
473  MagnitudeType convTol_, orthoKappa_, achievedTol_;
474  int maxRestarts_, maxIters_, numIters_;
475  int verbosity_, outputStyle_, outputFreq_;
476  std::string orthoType_;
477  std::string impResScale_, expResScale_;
478 
480  // Solver State Storage
482  //
483  // The number of blocks and recycle blocks (m and k, respectively)
484  int numBlocks_, recycledBlocks_;
485  // Current size of recycled subspace
486  int keff;
487  //
488  // Residual vector
489  Teuchos::RCP<MV> r_;
490  //
491  // Search space
492  Teuchos::RCP<MV> V_;
493  //
494  // Recycled subspace and its image
495  Teuchos::RCP<MV> U_, C_;
496  //
497  // Updated recycle space and its image
498  Teuchos::RCP<MV> U1_, C1_;
499  //
500  // Storage used in constructing harmonic Ritz values/vectors
506  std::vector<ScalarType> tau_;
507  std::vector<ScalarType> work_;
509  std::vector<int> ipiv_;
511 
512  // Timers.
513  std::string label_;
514  Teuchos::RCP<Teuchos::Time> timerSolve_;
515 
516  // Internal state variables.
517  bool isSet_;
518 
519  // Have we generated or regenerated a recycle space yet this solve?
520  bool builtRecycleSpace_;
521  };
522 
523 
524 // Empty Constructor
525 template<class ScalarType, class MV, class OP>
527  achievedTol_(0.0),
528  numIters_(0)
529 {
530  init ();
531 }
532 
533 
534 // Basic Constructor
535 template<class ScalarType, class MV, class OP>
539  achievedTol_(0.0),
540  numIters_(0)
541 {
542  // Initialize local pointers to null, and initialize local variables
543  // to default values.
544  init ();
545 
547  problem == Teuchos::null, std::invalid_argument,
548  "Belos::GCRODRSolMgr constructor: The solver manager's "
549  "constructor needs the linear problem argument 'problem' "
550  "to be non-null.");
551  problem_ = problem;
552 
553  // Set the parameters using the list that was passed in. If null,
554  // we defer initialization until a non-null list is set (by the
555  // client calling setParameters(), or by calling solve() -- in
556  // either case, a null parameter list indicates that default
557  // parameters should be used).
558  if (! pl.is_null ()) {
559  setParameters (pl);
560  }
561 }
562 
563 // Common instructions executed in all constructors
564 template<class ScalarType, class MV, class OP>
566  outputStream_ = Teuchos::rcpFromRef(std::cout);
568  orthoKappa_ = orthoKappa_default_;
569  maxRestarts_ = maxRestarts_default_;
570  maxIters_ = maxIters_default_;
571  numBlocks_ = numBlocks_default_;
572  recycledBlocks_ = recycledBlocks_default_;
573  verbosity_ = verbosity_default_;
574  outputStyle_ = outputStyle_default_;
575  outputFreq_ = outputFreq_default_;
576  orthoType_ = orthoType_default_;
577  impResScale_ = impResScale_default_;
578  expResScale_ = expResScale_default_;
579  label_ = label_default_;
580  isSet_ = false;
581  builtRecycleSpace_ = false;
582  keff = 0;
583  r_ = Teuchos::null;
584  V_ = Teuchos::null;
585  U_ = Teuchos::null;
586  C_ = Teuchos::null;
587  U1_ = Teuchos::null;
588  C1_ = Teuchos::null;
589  PP_ = Teuchos::null;
590  HP_ = Teuchos::null;
591  H2_ = Teuchos::null;
592  R_ = Teuchos::null;
593  H_ = Teuchos::null;
594  B_ = Teuchos::null;
595 }
596 
597 template<class ScalarType, class MV, class OP>
598 void
601 {
602  using Teuchos::isParameterType;
603  using Teuchos::getParameter;
604  using Teuchos::null;
606  using Teuchos::parameterList;
607  using Teuchos::RCP;
608  using Teuchos::rcp;
609  using Teuchos::rcp_dynamic_cast;
610  using Teuchos::rcpFromRef;
614 
615  // The default parameter list contains all parameters that
616  // GCRODRSolMgr understands, and none that it doesn't understand.
617  RCP<const ParameterList> defaultParams = getValidParameters();
618 
619  // Create the internal parameter list if one doesn't already exist.
620  //
621  // (mfh 28 Feb 2011, 10 Mar 2011) At the time this code was written,
622  // ParameterList did not have validators or validateParameters().
623  // This is why the code below carefully validates the parameters one
624  // by one and fills in defaults. This code could be made a lot
625  // shorter by using validators. To do so, we would have to define
626  // appropriate validators for all the parameters. (This would more
627  // or less just move all that validation code out of this routine
628  // into to getValidParameters().)
629  //
630  // For an analogous reason, GCRODRSolMgr defines default parameter
631  // values as class data, as well as in the default ParameterList.
632  // This redundancy could be removed by defining the default
633  // parameter values only in the default ParameterList (which
634  // documents each parameter as well -- handy!).
635  if (params_.is_null()) {
636  params_ = parameterList (*defaultParams);
637  } else {
638  // A common case for setParameters() is for it to be called at the
639  // beginning of the solve() routine. This follows the Belos
640  // pattern of delaying initialization until the last possible
641  // moment (when the user asks Belos to perform the solve). In
642  // this common case, we save ourselves a deep copy of the input
643  // parameter list.
644  if (params_ != params) {
645  // Make a deep copy of the input parameter list. This allows
646  // the caller to modify or change params later, without
647  // affecting the behavior of this solver. This solver will then
648  // only change its internal parameters if setParameters() is
649  // called again.
650  params_ = parameterList (*params);
651  }
652 
653  // Fill in any missing parameters and their default values. Also,
654  // throw an exception if the parameter list has any misspelled or
655  // "extra" parameters. If you don't like this behavior, you'll
656  // want to replace the line of code below with your desired
657  // validation scheme. Note that Teuchos currently only implements
658  // two options:
659  //
660  // 1. validateParameters() requires that params_ has all the
661  // parameters that the default list has, and none that it
662  // doesn't have.
663  //
664  // 2. validateParametersAndSetDefaults() fills in missing
665  // parameters in params_ using the default list, but requires
666  // that any parameter provided in params_ is also in the
667  // default list.
668  //
669  // Here is an easy way to ignore any "extra" or misspelled
670  // parameters: Make a deep copy of the default list, fill in any
671  // "missing" parameters from the _input_ list, and then validate
672  // the input list using the deep copy of the default list. We
673  // show this in code:
674  //
675  // RCP<ParameterList> defaultCopy = parameterList (*getValidParameters ());
676  // defaultCopy->validateParametersAndSetDefaults (params);
677  // params->validateParametersAndSetDefaults (defaultCopy);
678  //
679  // This method is not entirely robust, because the input list may
680  // have incorrect validators set for existing parameters in the
681  // default list. This would then cause "validation" of the
682  // default list to throw an exception. As a result, we've chosen
683  // for now to be intolerant of misspellings and "extra" parameters
684  // in the input list.
685  params_->validateParametersAndSetDefaults (*defaultParams);
686  }
687 
688  // Check for maximum number of restarts.
689  if (params->isParameter ("Maximum Restarts")) {
690  maxRestarts_ = params->get("Maximum Restarts", maxRestarts_default_);
691 
692  // Update parameter in our list.
693  params_->set ("Maximum Restarts", maxRestarts_);
694  }
695 
696  // Check for maximum number of iterations
697  if (params->isParameter ("Maximum Iterations")) {
698  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
699 
700  // Update parameter in our list and in status test.
701  params_->set ("Maximum Iterations", maxIters_);
702  if (! maxIterTest_.is_null())
703  maxIterTest_->setMaxIters (maxIters_);
704  }
705 
706  // Check for the maximum number of blocks.
707  if (params->isParameter ("Num Blocks")) {
708  numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
709  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
710  "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
711  "be strictly positive, but you specified a value of "
712  << numBlocks_ << ".");
713  // Update parameter in our list.
714  params_->set ("Num Blocks", numBlocks_);
715  }
716 
717  // Check for the maximum number of blocks.
718  if (params->isParameter ("Num Recycled Blocks")) {
719  recycledBlocks_ = params->get ("Num Recycled Blocks",
720  recycledBlocks_default_);
721  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
722  "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
723  "parameter must be strictly positive, but you specified "
724  "a value of " << recycledBlocks_ << ".");
725  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
726  "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
727  "parameter must be less than the \"Num Blocks\" "
728  "parameter, but you specified \"Num Recycled Blocks\" "
729  "= " << recycledBlocks_ << " and \"Num Blocks\" = "
730  << numBlocks_ << ".");
731  // Update parameter in our list.
732  params_->set("Num Recycled Blocks", recycledBlocks_);
733  }
734 
735  // Check to see if the timer label changed. If it did, update it in
736  // the parameter list, and create a new timer with that label (if
737  // Belos was compiled with timers enabled).
738  if (params->isParameter ("Timer Label")) {
739  std::string tempLabel = params->get ("Timer Label", label_default_);
740 
741  // Update parameter in our list and solver timer
742  if (tempLabel != label_) {
743  label_ = tempLabel;
744  params_->set ("Timer Label", label_);
745  std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
746 #ifdef BELOS_TEUCHOS_TIME_MONITOR
747  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
748 #endif
749  if (ortho_ != Teuchos::null) {
750  ortho_->setLabel( label_ );
751  }
752  }
753  }
754 
755  // Check for a change in verbosity level
756  if (params->isParameter ("Verbosity")) {
757  if (isParameterType<int> (*params, "Verbosity")) {
758  verbosity_ = params->get ("Verbosity", verbosity_default_);
759  } else {
760  verbosity_ = (int) getParameter<Belos::MsgType> (*params, "Verbosity");
761  }
762  // Update parameter in our list.
763  params_->set ("Verbosity", verbosity_);
764  // If the output manager (printer_) is null, then we will
765  // instantiate it later with the correct verbosity.
766  if (! printer_.is_null())
767  printer_->setVerbosity (verbosity_);
768  }
769 
770  // Check for a change in output style
771  if (params->isParameter ("Output Style")) {
772  if (isParameterType<int> (*params, "Output Style")) {
773  outputStyle_ = params->get ("Output Style", outputStyle_default_);
774  } else {
775  outputStyle_ = (int) getParameter<OutputType> (*params, "Output Style");
776  }
777 
778  // Update parameter in our list.
779  params_->set ("Output Style", outputStyle_);
780  // We will (re)instantiate the output status test afresh below.
781  outputTest_ = null;
782  }
783 
784  // Get the output stream for the output manager.
785  //
786  // While storing the output stream in the parameter list (either as
787  // an RCP or as a nonconst reference) is convenient and safe for
788  // programming, it makes it impossible to serialize the parameter
789  // list, read it back in from the serialized representation, and get
790  // the same output stream as before. This is because output streams
791  // may be arbitrary constructed objects.
792  //
793  // In case the user tried reading in the parameter list from a
794  // serialized representation and the output stream can't be read
795  // back in, we set the output stream to point to std::cout. This
796  // ensures reasonable behavior.
797  if (params->isParameter ("Output Stream")) {
798  try {
799  outputStream_ = getParameter<RCP<std::ostream> > (*params, "Output Stream");
800  } catch (InvalidParameter&) {
801  outputStream_ = rcpFromRef (std::cout);
802  }
803  // We assume that a null output stream indicates that the user
804  // doesn't want to print anything, so we replace it with a "black
805  // hole" stream that prints nothing sent to it. (We can't use a
806  // null output stream, since the output manager always sends
807  // things it wants to print to the output stream.)
808  if (outputStream_.is_null()) {
809  outputStream_ = rcp (new Teuchos::oblackholestream);
810  }
811  // Update parameter in our list.
812  params_->set ("Output Stream", outputStream_);
813  // If the output manager (printer_) is null, then we will
814  // instantiate it later with the correct output stream.
815  if (! printer_.is_null()) {
816  printer_->setOStream (outputStream_);
817  }
818  }
819 
820  // frequency level
821  if (verbosity_ & Belos::StatusTestDetails) {
822  if (params->isParameter ("Output Frequency")) {
823  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
824  }
825 
826  // Update parameter in out list and output status test.
827  params_->set("Output Frequency", outputFreq_);
828  if (! outputTest_.is_null())
829  outputTest_->setOutputFrequency (outputFreq_);
830  }
831 
832  // Create output manager if we need to, using the verbosity level
833  // and output stream that we fetched above. We do this here because
834  // instantiating an OrthoManager using OrthoManagerFactory requires
835  // a valid OutputManager.
836  if (printer_.is_null()) {
837  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
838  }
839 
840  // Get the orthogonalization manager name ("Orthogonalization").
841  //
842  // Getting default values for the orthogonalization manager
843  // parameters ("Orthogonalization Parameters") requires knowing the
844  // orthogonalization manager name. Save it for later, and also
845  // record whether it's different than before.
847  bool changedOrthoType = false;
848  if (params->isParameter ("Orthogonalization")) {
849  const std::string& tempOrthoType =
850  params->get ("Orthogonalization", orthoType_default_);
851  // Ensure that the specified orthogonalization type is valid.
852  if (! factory.isValidName (tempOrthoType)) {
853  std::ostringstream os;
854  os << "Belos::GCRODRSolMgr: Invalid orthogonalization name \""
855  << tempOrthoType << "\". The following are valid options "
856  << "for the \"Orthogonalization\" name parameter: ";
857  factory.printValidNames (os);
858  throw std::invalid_argument (os.str());
859  }
860  if (tempOrthoType != orthoType_) {
861  changedOrthoType = true;
862  orthoType_ = tempOrthoType;
863  // Update parameter in our list.
864  params_->set ("Orthogonalization", orthoType_);
865  }
866  }
867 
868  // Get any parameters for the orthogonalization ("Orthogonalization
869  // Parameters"). If not supplied, the orthogonalization manager
870  // factory will supply default values.
871  //
872  // NOTE (mfh 12 Jan 2011) For the sake of backwards compatibility,
873  // if params has an "Orthogonalization Constant" parameter and the
874  // DGKS orthogonalization manager is to be used, the value of this
875  // parameter will override DGKS's "depTol" parameter.
876  //
877  // Users must supply the orthogonalization manager parameters as a
878  // sublist (supplying it as an RCP<ParameterList> would make the
879  // resulting parameter list not serializable).
880  RCP<ParameterList> orthoParams;
881  { // The nonmember function sublist() returns an RCP<ParameterList>,
882  // which is what we want here.
883  using Teuchos::sublist;
884  // Abbreviation to avoid typos.
885  const std::string paramName ("Orthogonalization Parameters");
886 
887  try {
888  orthoParams = sublist (params_, paramName, true);
889  } catch (InvalidParameter&) {
890  // We didn't get the parameter list from params, so get a
891  // default parameter list from the OrthoManagerFactory. Modify
892  // params_ so that it has the default parameter list, and set
893  // orthoParams to ensure it's a sublist of params_ (and not just
894  // a copy of one).
895  params_->set (paramName, factory.getDefaultParameters (orthoType_));
896  orthoParams = sublist (params_, paramName, true);
897  }
898  }
899  TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
900  "Failed to get orthogonalization parameters. "
901  "Please report this bug to the Belos developers.");
902 
903  // Instantiate a new MatOrthoManager subclass instance if necessary.
904  // If not necessary, then tell the existing instance about the new
905  // parameters.
906  if (ortho_.is_null() || changedOrthoType) {
907  // We definitely need to make a new MatOrthoManager, since either
908  // we haven't made one yet, or we've changed orthogonalization
909  // methods. Creating the orthogonalization manager requires that
910  // the OutputManager (printer_) already be initialized.
911  ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_,
912  label_, orthoParams);
913  } else {
914  // If the MatOrthoManager implements the ParameterListAcceptor
915  // mix-in interface, we can propagate changes to its parameters
916  // without reinstantiating the MatOrthoManager.
917  //
918  // We recommend that all MatOrthoManager subclasses implement
919  // Teuchos::ParameterListAcceptor, but do not (yet) require this.
920  typedef Teuchos::ParameterListAcceptor PLA;
921  RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
922  if (pla.is_null()) {
923  // Oops, it's not a ParameterListAcceptor. We have to
924  // reinstantiate the MatOrthoManager in order to pass in the
925  // possibly new parameters.
926  ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_,
927  label_, orthoParams);
928  } else {
929  pla->setParameterList (orthoParams);
930  }
931  }
932 
933  // The DGKS orthogonalization accepts a "Orthogonalization Constant"
934  // parameter (also called kappa in the code, but not in the
935  // parameter list). If its value is provided in the given parameter
936  // list, and its value is positive, use it. Ignore negative values.
937  //
938  // NOTE (mfh 12 Jan 2011) This overrides the "depTol" parameter that
939  // may have been specified in "Orthogonalization Parameters". We
940  // retain this behavior for backwards compatibility.
941  if (params->isParameter ("Orthogonalization Constant")) {
942  MagnitudeType orthoKappa = orthoKappa_default_;
943  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
944  orthoKappa = params->get ("Orthogonalization Constant", orthoKappa);
945  }
946  else {
947  orthoKappa = params->get ("Orthogonalization Constant", orthoKappa_default_);
948  }
949 
950  if (orthoKappa > 0) {
951  orthoKappa_ = orthoKappa;
952  // Update parameter in our list.
953  params_->set("Orthogonalization Constant", orthoKappa_);
954  // Only DGKS currently accepts this parameter.
955  if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
956  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
957  // This cast should always succeed; it's a bug
958  // otherwise. (If the cast fails, then orthoType_
959  // doesn't correspond to the OrthoManager subclass
960  // instance that we think we have, so we initialized the
961  // wrong subclass somehow.)
962  rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
963  }
964  }
965  }
966 
967  // Convergence
968  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
969  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
970 
971  // Check for convergence tolerance
972  if (params->isParameter("Convergence Tolerance")) {
973  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
974  convTol_ = params->get ("Convergence Tolerance",
975  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
976  }
977  else {
978  convTol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
979  }
980 
981  // Update parameter in our list and residual tests.
982  params_->set ("Convergence Tolerance", convTol_);
983  if (! impConvTest_.is_null())
984  impConvTest_->setTolerance (convTol_);
985  if (! expConvTest_.is_null())
986  expConvTest_->setTolerance (convTol_);
987  }
988 
989  // Check for a change in scaling, if so we need to build new residual tests.
990  if (params->isParameter ("Implicit Residual Scaling")) {
991  std::string tempImpResScale =
992  getParameter<std::string> (*params, "Implicit Residual Scaling");
993 
994  // Only update the scaling if it's different.
995  if (impResScale_ != tempImpResScale) {
996  ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
997  impResScale_ = tempImpResScale;
998 
999  // Update parameter in our list and residual tests
1000  params_->set("Implicit Residual Scaling", impResScale_);
1001  // NOTE (mfh 28 Feb 2011) StatusTestImpResNorm only lets you
1002  // call defineScaleForm() once. The code below attempts to call
1003  // defineScaleForm(); if the scale form has already been
1004  // defined, it constructs a new StatusTestImpResNorm instance.
1005  // StatusTestImpResNorm really should not expose the
1006  // defineScaleForm() method, since it's serving an
1007  // initialization purpose; all initialization should happen in
1008  // the constructor whenever possible. In that case, the code
1009  // below could be simplified into a single (re)instantiation.
1010  if (! impConvTest_.is_null()) {
1011  try {
1012  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
1013  }
1014  catch (StatusTestError&) {
1015  // Delete the convergence test so it gets constructed again.
1016  impConvTest_ = null;
1017  convTest_ = null;
1018  }
1019  }
1020  }
1021  }
1022 
1023  if (params->isParameter("Explicit Residual Scaling")) {
1024  std::string tempExpResScale =
1025  getParameter<std::string> (*params, "Explicit Residual Scaling");
1026 
1027  // Only update the scaling if it's different.
1028  if (expResScale_ != tempExpResScale) {
1029  ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
1030  expResScale_ = tempExpResScale;
1031 
1032  // Update parameter in our list and residual tests
1033  params_->set("Explicit Residual Scaling", expResScale_);
1034  // NOTE (mfh 28 Feb 2011) See note above on the (re)construction
1035  // of StatusTestImpResNorm.
1036  if (! expConvTest_.is_null()) {
1037  try {
1038  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
1039  }
1040  catch (StatusTestError&) {
1041  // Delete the convergence test so it gets constructed again.
1042  expConvTest_ = null;
1043  convTest_ = null;
1044  }
1045  }
1046  }
1047  }
1048  //
1049  // Create iteration stopping criteria ("status tests") if we need
1050  // to, by combining three different stopping criteria.
1051  //
1052  // First, construct maximum-number-of-iterations stopping criterion.
1053  if (maxIterTest_.is_null())
1054  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
1055 
1056  // Implicit residual test, using the native residual to determine if
1057  // convergence was achieved.
1058  if (impConvTest_.is_null()) {
1059  impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1060  impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
1061  Belos::TwoNorm);
1062  }
1063 
1064  // Explicit residual test once the native residual is below the tolerance
1065  if (expConvTest_.is_null()) {
1066  expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1067  expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
1068  expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
1069  Belos::TwoNorm);
1070  }
1071  // Convergence test first tests the implicit residual, then the
1072  // explicit residual if the implicit residual test passes.
1073  if (convTest_.is_null()) {
1074  convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1075  impConvTest_,
1076  expConvTest_));
1077  }
1078  // Construct the complete iteration stopping criterion:
1079  //
1080  // "Stop iterating if the maximum number of iterations has been
1081  // reached, or if the convergence test passes."
1082  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
1083  maxIterTest_,
1084  convTest_));
1085  // Create the status test output class.
1086  // This class manages and formats the output from the status test.
1087  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
1088  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
1090 
1091  // Set the solver string for the output test
1092  std::string solverDesc = " GCRODR ";
1093  outputTest_->setSolverDesc( solverDesc );
1094 
1095  // Create the timer if we need to.
1096  if (timerSolve_.is_null()) {
1097  std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
1098 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1099  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1100 #endif
1101  }
1102 
1103  // Inform the solver manager that the current parameters were set.
1104  isSet_ = true;
1105 }
1106 
1107 
1108 template<class ScalarType, class MV, class OP>
1111 {
1112  using Teuchos::ParameterList;
1113  using Teuchos::parameterList;
1114  using Teuchos::RCP;
1115 
1116  static RCP<const ParameterList> validPL;
1117  if (is_null(validPL)) {
1118  RCP<ParameterList> pl = parameterList ();
1119 
1120  // Set all the valid parameters and their default values.
1121  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
1122  "The relative residual tolerance that needs to be achieved by the\n"
1123  "iterative solver in order for the linear system to be declared converged.");
1124  pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
1125  "The maximum number of cycles allowed for each\n"
1126  "set of RHS solved.");
1127  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
1128  "The maximum number of iterations allowed for each\n"
1129  "set of RHS solved.");
1130  // mfh 25 Oct 2010: "Block Size" must be 1 because GCRODR is
1131  // currently not a block method: i.e., it does not work on
1132  // multiple right-hand sides at once.
1133  pl->set("Block Size", static_cast<int>(blockSize_default_),
1134  "Block Size Parameter -- currently must be 1 for GCRODR");
1135  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
1136  "The maximum number of vectors allowed in the Krylov subspace\n"
1137  "for each set of RHS solved.");
1138  pl->set("Num Recycled Blocks", static_cast<int>(recycledBlocks_default_),
1139  "The maximum number of vectors in the recycled subspace." );
1140  pl->set("Verbosity", static_cast<int>(verbosity_default_),
1141  "What type(s) of solver information should be outputted\n"
1142  "to the output stream.");
1143  pl->set("Output Style", static_cast<int>(outputStyle_default_),
1144  "What style is used for the solver information outputted\n"
1145  "to the output stream.");
1146  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
1147  "How often convergence information should be outputted\n"
1148  "to the output stream.");
1149  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
1150  "A reference-counted pointer to the output stream where all\n"
1151  "solver output is sent.");
1152  pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1153  "The type of scaling used in the implicit residual convergence test.");
1154  pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1155  "The type of scaling used in the explicit residual convergence test.");
1156  pl->set("Timer Label", static_cast<const char *>(label_default_),
1157  "The string to use as a prefix for the timer labels.");
1158  {
1160  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
1161  "The type of orthogonalization to use. Valid options: " +
1162  factory.validNamesString());
1163  RCP<const ParameterList> orthoParams =
1164  factory.getDefaultParameters (orthoType_default_);
1165  pl->set ("Orthogonalization Parameters", *orthoParams,
1166  "Parameters specific to the type of orthogonalization used.");
1167  }
1168  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(orthoKappa_default_),
1169  "When using DGKS orthogonalization: the \"depTol\" constant, used "
1170  "to determine whether another step of classical Gram-Schmidt is "
1171  "necessary. Otherwise ignored.");
1172  validPL = pl;
1173  }
1174  return validPL;
1175 }
1176 
1177 // initializeStateStorage
1178 template<class ScalarType, class MV, class OP>
1180 
1181  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1182 
1183  // Check if there is any multivector to clone from.
1184  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1185  if (rhsMV == Teuchos::null) {
1186  // Nothing to do
1187  return;
1188  }
1189  else {
1190 
1191  // Initialize the state storage
1192  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1193  "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1194 
1195  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1196  if (U_ == Teuchos::null) {
1197  U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1198  }
1199  else {
1200  // Generate U_ by cloning itself ONLY if more space is needed.
1201  if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1202  Teuchos::RCP<const MV> tmp = U_;
1203  U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1204  }
1205  }
1206 
1207  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1208  if (C_ == Teuchos::null) {
1209  C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1210  }
1211  else {
1212  // Generate C_ by cloning itself ONLY if more space is needed.
1213  if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1214  Teuchos::RCP<const MV> tmp = C_;
1215  C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1216  }
1217  }
1218 
1219  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1220  if (V_ == Teuchos::null) {
1221  V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1222  }
1223  else {
1224  // Generate V_ by cloning itself ONLY if more space is needed.
1225  if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1226  Teuchos::RCP<const MV> tmp = V_;
1227  V_ = MVT::Clone( *tmp, numBlocks_+1 );
1228  }
1229  }
1230 
1231  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1232  if (U1_ == Teuchos::null) {
1233  U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1234  }
1235  else {
1236  // Generate U1_ by cloning itself ONLY if more space is needed.
1237  if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1238  Teuchos::RCP<const MV> tmp = U1_;
1239  U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1240  }
1241  }
1242 
1243  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1244  if (C1_ == Teuchos::null) {
1245  C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1246  }
1247  else {
1248  // Generate C1_ by cloning itself ONLY if more space is needed.
1249  if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1250  Teuchos::RCP<const MV> tmp = C1_;
1251  C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1252  }
1253  }
1254 
1255  // Generate r_ only if it doesn't exist
1256  if (r_ == Teuchos::null)
1257  r_ = MVT::Clone( *rhsMV, 1 );
1258 
1259  // Size of tau_ will change during computation, so just be sure it starts with appropriate size
1260  tau_.resize(recycledBlocks_+1);
1261 
1262  // Size of work_ will change during computation, so just be sure it starts with appropriate size
1263  work_.resize(recycledBlocks_+1);
1264 
1265  // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
1266  ipiv_.resize(recycledBlocks_+1);
1267 
1268  // Generate H2_ only if it doesn't exist, otherwise resize it.
1269  if (H2_ == Teuchos::null)
1270  H2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1271  else {
1272  if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1273  H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1274  }
1275  H2_->putScalar(zero);
1276 
1277  // Generate R_ only if it doesn't exist, otherwise resize it.
1278  if (R_ == Teuchos::null)
1279  R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1280  else {
1281  if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1282  R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1283  }
1284  R_->putScalar(zero);
1285 
1286  // Generate PP_ only if it doesn't exist, otherwise resize it.
1287  if (PP_ == Teuchos::null)
1288  PP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1289  else {
1290  if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1291  PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1292  }
1293 
1294  // Generate HP_ only if it doesn't exist, otherwise resize it.
1295  if (HP_ == Teuchos::null)
1296  HP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1297  else {
1298  if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1299  HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1300  }
1301 
1302  } // end else
1303 }
1304 
1305 
1306 // solve()
1307 template<class ScalarType, class MV, class OP>
1309  using Teuchos::RCP;
1310  using Teuchos::rcp;
1311 
1312  // Set the current parameters if they were not set before.
1313  // NOTE: This may occur if the user generated the solver manager with the default constructor and
1314  // then didn't set any parameters using setParameters().
1315  if (!isSet_) { setParameters( params_ ); }
1316 
1317  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1318  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1319  std::vector<int> index(numBlocks_+1);
1320 
1321  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1322 
1323  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1324 
1325  // Create indices for the linear systems to be solved.
1326  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1327  std::vector<int> currIdx(1);
1328  currIdx[0] = 0;
1329 
1330  // Inform the linear problem of the current linear system to solve.
1331  problem_->setLSIndex( currIdx );
1332 
1333  // Check the number of blocks and change them is necessary.
1334  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1335  if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1336  numBlocks_ = Teuchos::as<int>(dim);
1337  printer_->stream(Warnings) <<
1338  "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1339  " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1340  params_->set("Num Blocks", numBlocks_);
1341  }
1342 
1343  // Assume convergence is achieved, then let any failed convergence set this to false.
1344  bool isConverged = true;
1345 
1346  // Initialize storage for all state variables
1347  initializeStateStorage();
1348 
1350  // Parameter list
1351  Teuchos::ParameterList plist;
1352 
1353  plist.set("Num Blocks",numBlocks_);
1354  plist.set("Recycled Blocks",recycledBlocks_);
1355 
1357  // GCRODR solver
1358 
1359  RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1360  gcrodr_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
1361  // Number of iterations required to generate initial recycle space (if needed)
1362  int prime_iterations = 0;
1363 
1364  // Enter solve() iterations
1365  {
1366 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1367  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1368 #endif
1369 
1370  while ( numRHS2Solve > 0 ) {
1371 
1372  // Set flag indicating recycle space has not been generated this solve
1373  builtRecycleSpace_ = false;
1374 
1375  // Reset the status test.
1376  outputTest_->reset();
1377 
1379  // Initialize recycled subspace for GCRODR
1380 
1381  // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
1382  if (keff > 0) {
1384  "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1385 
1386  printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
1387  // Compute image of U_ under the new operator
1388  index.resize(keff);
1389  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1390  RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1391  RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1392  problem_->apply( *Utmp, *Ctmp );
1393 
1394  RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1395 
1396  // Orthogonalize this block
1397  // Get a matrix to hold the orthonormalization coefficients.
1399  int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,false));
1400  // Throw an error if we could not orthogonalize this block
1401  TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1402 
1403  // U_ = U_*R^{-1}
1404  // First, compute LU factorization of R
1405  int info = 0;
1406  ipiv_.resize(Rtmp.numRows());
1407  lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1408  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1409 
1410  // Now, form inv(R)
1411  int lwork = Rtmp.numRows();
1412  work_.resize(lwork);
1413  lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1414  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1415 
1416  // U_ = U1_; (via a swap)
1417  MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1418  std::swap(U_, U1_);
1419 
1420  // Must reinitialize after swap
1421  index.resize(keff);
1422  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1423  Ctmp = MVT::CloneViewNonConst( *C_, index );
1424  Utmp = MVT::CloneView( *U_, index );
1425 
1426  // Compute C_'*r_
1428  problem_->computeCurrPrecResVec( &*r_ );
1429  MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1430 
1431  // Update solution ( x += U_*C_'*r_ )
1432  RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1433  MVT::MvInit( *update, 0.0 );
1434  MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1435  problem_->updateSolution( update, true );
1436 
1437  // Update residual norm ( r -= C_*C_'*r_ )
1438  MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1439 
1440  // We recycled space from previous call
1441  prime_iterations = 0;
1442 
1443  }
1444  else {
1445 
1446  // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
1447  printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1448 
1449  Teuchos::ParameterList primeList;
1450 
1451  // Tell the block solver that the block size is one.
1452  primeList.set("Num Blocks",numBlocks_);
1453  primeList.set("Recycled Blocks",0);
1454 
1455  // Create GCRODR iterator object to perform one cycle of GMRES.
1456  RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1457  gcrodr_prime_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
1458 
1459  // Create the first block in the current Krylov basis (residual).
1460  problem_->computeCurrPrecResVec( &*r_ );
1461  index.resize( 1 ); index[0] = 0;
1462  RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1463  MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
1464 
1465  // Set the new state and initialize the solver.
1467  index.resize( numBlocks_+1 );
1468  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1469  newstate.V = MVT::CloneViewNonConst( *V_, index );
1470  newstate.U = Teuchos::null;
1471  newstate.C = Teuchos::null;
1472  newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1473  newstate.B = Teuchos::null;
1474  newstate.curDim = 0;
1475  gcrodr_prime_iter->initialize(newstate);
1476 
1477  // Perform one cycle of GMRES
1478  bool primeConverged = false;
1479  try {
1480  gcrodr_prime_iter->iterate();
1481 
1482  // Check convergence first
1483  if ( convTest_->getStatus() == Passed ) {
1484  // we have convergence
1485  primeConverged = true;
1486  }
1487  }
1488  catch (const GCRODRIterOrthoFailure &e) {
1489  // Try to recover the most recent least-squares solution
1490  gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1491 
1492  // Check to see if the most recent least-squares solution yielded convergence.
1493  sTest_->checkStatus( &*gcrodr_prime_iter );
1494  if (convTest_->getStatus() == Passed)
1495  primeConverged = true;
1496  }
1497  catch (const StatusTestNaNError& e) {
1498  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
1499  achievedTol_ = MT::one();
1500  Teuchos::RCP<MV> X = problem_->getLHS();
1501  MVT::MvInit( *X, SCT::zero() );
1502  printer_->stream(Warnings) << "Belos::GCRODRSolMgr::solve(): Warning! NaN has been detected!"
1503  << std::endl;
1504  return Unconverged;
1505  }
1506  catch (const std::exception &e) {
1507  printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration "
1508  << gcrodr_prime_iter->getNumIters() << std::endl
1509  << e.what() << std::endl;
1510  throw;
1511  }
1512  // Record number of iterations in generating initial recycle spacec
1513  prime_iterations = gcrodr_prime_iter->getNumIters();
1514 
1515  // Update the linear problem.
1516  RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1517  problem_->updateSolution( update, true );
1518 
1519  // Get the state.
1520  newstate = gcrodr_prime_iter->getState();
1521  int p = newstate.curDim;
1522 
1523  // Compute harmonic Ritz vectors
1524  // NOTE: The storage for the harmonic Ritz vectors (PP) is made one column larger
1525  // just in case we split a complex conjugate pair.
1526  // NOTE: Generate a recycled subspace only if we have enough vectors. If we converged
1527  // too early, move on to the next linear system and try to generate a subspace again.
1528  if (recycledBlocks_ < p+1) {
1529  int info = 0;
1530  RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1531  // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
1532  keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
1533  // Hereafter, only keff columns of PP are needed
1534  PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1535  // Now get views into C, U, V
1536  index.resize(keff);
1537  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1538  RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1539  RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1540  RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1541  index.resize(p);
1542  for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
1543  RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1544 
1545  // Form U (the subspace to recycle)
1546  // U = newstate.V(:,1:p) * PP;
1547  MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1548 
1549  // Form orthonormalized C and adjust U so that C = A*U
1550 
1551  // First, compute [Q, R] = qr(H*P);
1552 
1553  // Step #1: Form HP = H*P
1554  Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1556  HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1557 
1558  // Step #1.5: Perform workspace size query for QR
1559  // factorization of HP. On input, lwork must be -1.
1560  // _GEQRF will put the workspace size in work_[0].
1561  int lwork = -1;
1562  tau_.resize (keff);
1563  lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1564  HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1566  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve:"
1567  " LAPACK's _GEQRF failed to compute a workspace size.");
1568 
1569  // Step #2: Compute QR factorization of HP
1570  //
1571  // NOTE (mfh 17 Apr 2014) LAPACK promises that the value of
1572  // work_[0] after the workspace query will fit in int. This
1573  // justifies the cast. We call real() first because
1574  // static_cast from std::complex to int doesn't work.
1575  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1576  work_.resize (lwork); // Allocate workspace for the QR factorization
1577  lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1578  HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1580  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve:"
1581  " LAPACK's _GEQRF failed to compute a QR factorization.");
1582 
1583  // Step #3: Explicitly construct Q and R factors
1584  // NOTE: The upper triangular part of HP is copied into R and HP becomes Q.
1586  for (int ii = 0; ii < keff; ++ii) {
1587  for (int jj = ii; jj < keff; ++jj) {
1588  Rtmp(ii,jj) = HPtmp(ii,jj);
1589  }
1590  }
1591  // NOTE (mfh 17 Apr 2014): Teuchos::LAPACK's wrapper for
1592  // UNGQR dispatches to the correct Scalar-specific routine.
1593  // It calls {S,D}ORGQR if Scalar is real, and {C,Z}UNGQR if
1594  // Scalar is complex.
1595  lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1596  HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1597  lwork, &info);
1599  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1600  "LAPACK's _UNGQR failed to construct the Q factor.");
1601 
1602  // Now we have [Q,R] = qr(H*P)
1603 
1604  // Now compute C = V(:,1:p+1) * Q
1605  index.resize (p + 1);
1606  for (int ii = 0; ii < (p+1); ++ii) {
1607  index[ii] = ii;
1608  }
1609  Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
1610  MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1611 
1612  // Finally, compute U = U*R^{-1}.
1613  // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
1614  // backsolve capabilities don't exist in the Belos::MultiVec class
1615 
1616  // Step #1: First, compute LU factorization of R
1617  ipiv_.resize(Rtmp.numRows());
1618  lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1620  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1621  "LAPACK's _GETRF failed to compute an LU factorization.");
1622 
1623  // FIXME (mfh 17 Apr 2014) We have to compute the explicit
1624  // inverse of R here because Belos::MultiVecTraits doesn't
1625  // have a triangular solve (where the triangular matrix is
1626  // globally replicated and the "right-hand side" is the
1627  // distributed MultiVector).
1628 
1629  // Step #2: Form inv(R)
1630  lwork = Rtmp.numRows();
1631  work_.resize(lwork);
1632  lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1634  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1635  "LAPACK's _GETRI failed to invert triangular matrix.");
1636 
1637  // Step #3: Let U = U * R^{-1}
1638  MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1639 
1640  printer_->stream(Debug)
1641  << " Generated recycled subspace using RHS index " << currIdx[0]
1642  << " of dimension " << keff << std::endl << std::endl;
1643 
1644  } // if (recycledBlocks_ < p+1)
1645 
1646  // Return to outer loop if the priming solve converged, set the next linear system.
1647  if (primeConverged) {
1648  // Inform the linear problem that we are finished with this block linear system.
1649  problem_->setCurrLS();
1650 
1651  // Update indices for the linear systems to be solved.
1652  numRHS2Solve -= 1;
1653  if (numRHS2Solve > 0) {
1654  currIdx[0]++;
1655  problem_->setLSIndex (currIdx); // Set the next indices
1656  }
1657  else {
1658  currIdx.resize (numRHS2Solve);
1659  }
1660 
1661  continue;
1662  }
1663  } // if (keff > 0) ...
1664 
1665  // Prepare for the Gmres iterations with the recycled subspace.
1666 
1667  // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
1668  gcrodr_iter->setSize( keff, numBlocks_ );
1669 
1670  // Reset the number of iterations.
1671  gcrodr_iter->resetNumIters(prime_iterations);
1672 
1673  // Reset the number of calls that the status test output knows about.
1674  outputTest_->resetNumCalls();
1675 
1676  // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
1677  problem_->computeCurrPrecResVec( &*r_ );
1678  index.resize( 1 ); index[0] = 0;
1679  RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1680  MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
1681 
1682  // Set the new state and initialize the solver.
1684  index.resize( numBlocks_+1 );
1685  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1686  newstate.V = MVT::CloneViewNonConst( *V_, index );
1687  index.resize( keff );
1688  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1689  newstate.C = MVT::CloneViewNonConst( *C_, index );
1690  newstate.U = MVT::CloneViewNonConst( *U_, index );
1691  newstate.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1692  newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1693  newstate.curDim = 0;
1694  gcrodr_iter->initialize(newstate);
1695 
1696  // variables needed for inner loop
1697  int numRestarts = 0;
1698  while(1) {
1699 
1700  // tell gcrodr_iter to iterate
1701  try {
1702  gcrodr_iter->iterate();
1703 
1705  //
1706  // check convergence first
1707  //
1709  if ( convTest_->getStatus() == Passed ) {
1710  // we have convergence
1711  break; // break from while(1){gcrodr_iter->iterate()}
1712  }
1714  //
1715  // check for maximum iterations
1716  //
1718  else if ( maxIterTest_->getStatus() == Passed ) {
1719  // we don't have convergence
1720  isConverged = false;
1721  break; // break from while(1){gcrodr_iter->iterate()}
1722  }
1724  //
1725  // check for restarting, i.e. the subspace is full
1726  //
1728  else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1729 
1730  // Update the recycled subspace even if we have hit the maximum number of restarts.
1731 
1732  // Update the linear problem.
1733  RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1734  problem_->updateSolution( update, true );
1735 
1736  buildRecycleSpace2(gcrodr_iter);
1737 
1738  printer_->stream(Debug)
1739  << " Generated new recycled subspace using RHS index "
1740  << currIdx[0] << " of dimension " << keff << std::endl
1741  << std::endl;
1742 
1743  // NOTE: If we have hit the maximum number of restarts then quit
1744  if (numRestarts >= maxRestarts_) {
1745  isConverged = false;
1746  break; // break from while(1){gcrodr_iter->iterate()}
1747  }
1748  numRestarts++;
1749 
1750  printer_->stream(Debug)
1751  << " Performing restart number " << numRestarts << " of "
1752  << maxRestarts_ << std::endl << std::endl;
1753 
1754  // Create the restart vector (first block in the current Krylov basis)
1755  problem_->computeCurrPrecResVec( &*r_ );
1756  index.resize( 1 ); index[0] = 0;
1757  RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1758  MVT::SetBlock(*r_,index,*v00); // V(:,0) = r
1759 
1760  // Set the new state and initialize the solver.
1761  GCRODRIterState<ScalarType,MV> restartState;
1762  index.resize( numBlocks_+1 );
1763  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1764  restartState.V = MVT::CloneViewNonConst( *V_, index );
1765  index.resize( keff );
1766  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1767  restartState.U = MVT::CloneViewNonConst( *U_, index );
1768  restartState.C = MVT::CloneViewNonConst( *C_, index );
1769  restartState.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1770  restartState.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1771  restartState.curDim = 0;
1772  gcrodr_iter->initialize(restartState);
1773 
1774 
1775  } // end of restarting
1776 
1778  //
1779  // we returned from iterate(), but none of our status tests Passed.
1780  // something is wrong, and it is probably our fault.
1781  //
1783 
1784  else {
1786  true, std::logic_error, "Belos::GCRODRSolMgr::solve: "
1787  "Invalid return from GCRODRIter::iterate().");
1788  }
1789  }
1790  catch (const GCRODRIterOrthoFailure &e) {
1791  // Try to recover the most recent least-squares solution
1792  gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1793 
1794  // Check to see if the most recent least-squares solution yielded convergence.
1795  sTest_->checkStatus( &*gcrodr_iter );
1796  if (convTest_->getStatus() != Passed)
1797  isConverged = false;
1798  break;
1799  }
1800  catch (const std::exception& e) {
1801  printer_->stream(Errors)
1802  << "Error! Caught exception in GCRODRIter::iterate() at iteration "
1803  << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1804  throw;
1805  }
1806  }
1807 
1808  // Compute the current solution.
1809  // Update the linear problem.
1810  RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1811  problem_->updateSolution( update, true );
1812 
1813  // Inform the linear problem that we are finished with this block linear system.
1814  problem_->setCurrLS();
1815 
1816  // If we didn't build a recycle space this solve but ran at least k iterations,
1817  // force build of new recycle space
1818 
1819  if (!builtRecycleSpace_) {
1820  buildRecycleSpace2(gcrodr_iter);
1821  printer_->stream(Debug)
1822  << " Generated new recycled subspace using RHS index " << currIdx[0]
1823  << " of dimension " << keff << std::endl << std::endl;
1824  }
1825 
1826  // Update indices for the linear systems to be solved.
1827  numRHS2Solve -= 1;
1828  if (numRHS2Solve > 0) {
1829  currIdx[0]++;
1830  problem_->setLSIndex (currIdx); // Set the next indices
1831  }
1832  else {
1833  currIdx.resize (numRHS2Solve);
1834  }
1835  } // while (numRHS2Solve > 0)
1836  }
1837 
1838  sTest_->print (printer_->stream (FinalSummary)); // print final summary
1839 
1840  // print timing information
1841 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1842  // Calling summarize() can be expensive, so don't call unless the
1843  // user wants to print out timing details. summarize() will do all
1844  // the work even if it's passed a "black hole" output stream.
1845  if (verbosity_ & TimingDetails)
1846  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1847 #endif // BELOS_TEUCHOS_TIME_MONITOR
1848 
1849  // get iteration information for this solve
1850  numIters_ = maxIterTest_->getNumIters ();
1851 
1852  // Save the convergence test value ("achieved tolerance") for this
1853  // solve. This solver (unlike BlockGmresSolMgr) always has two
1854  // residual norm status tests: an explicit and an implicit test.
1855  // The master convergence test convTest_ is a SEQ combo of the
1856  // implicit resp. explicit tests. If the implicit test never
1857  // passes, then the explicit test won't ever be executed. This
1858  // manifests as expConvTest_->getTestValue()->size() < 1. We deal
1859  // with this case by using the values returned by
1860  // impConvTest_->getTestValue().
1861  {
1862  const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1863  if (pTestValues == NULL || pTestValues->size() < 1) {
1864  pTestValues = impConvTest_->getTestValue();
1865  }
1866  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1867  "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1868  "method returned NULL. Please report this bug to the Belos developers.");
1869  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1870  "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1871  "method returned a vector of length zero. Please report this bug to the "
1872  "Belos developers.");
1873 
1874  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1875  // achieved tolerances for all vectors in the current solve(), or
1876  // just for the vectors from the last deflation?
1877  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1878  }
1879 
1880  return isConverged ? Converged : Unconverged; // return from solve()
1881 }
1882 
1883 // Given existing recycle space and Krylov space, build new recycle space
1884 template<class ScalarType, class MV, class OP>
1886 
1887  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1888  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1889 
1890  std::vector<MagnitudeType> d(keff);
1891  std::vector<ScalarType> dscalar(keff);
1892  std::vector<int> index(numBlocks_+1);
1893 
1894  // Get the state
1895  GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
1896  int p = oldState.curDim;
1897 
1898  // insufficient new information to update recycle space
1899  if (p<1) return;
1900 
1901  // Take the norm of the recycled vectors
1902  {
1903  index.resize(keff);
1904  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1905  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1906  d.resize(keff);
1907  dscalar.resize(keff);
1908  MVT::MvNorm( *Utmp, d );
1909  for (int i=0; i<keff; ++i) {
1910  d[i] = one / d[i];
1911  dscalar[i] = (ScalarType)d[i];
1912  }
1913  MVT::MvScale( *Utmp, dscalar );
1914  }
1915 
1916  // Get view into current "full" upper Hessnburg matrix
1918 
1919  // Insert D into the leading keff x keff block of H2
1920  for (int i=0; i<keff; ++i) {
1921  (*H2tmp)(i,i) = d[i];
1922  }
1923 
1924  // Compute the harmoic Ritz pairs for the generalized eigenproblem
1925  // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
1926  int keff_new;
1927  {
1928  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1929  keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, PPtmp );
1930  }
1931 
1932  // Code to form new U, C
1933  // U = [U V(:,1:p)] * P; (in two steps)
1934 
1935  // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
1936  Teuchos::RCP<MV> U1tmp;
1937  {
1938  index.resize( keff );
1939  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1940  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1941  index.resize( keff_new );
1942  for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1943  U1tmp = MVT::CloneViewNonConst( *U1_, index );
1944  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1945  MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1946  }
1947 
1948  // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
1949  {
1950  index.resize(p);
1951  for (int ii=0; ii < p; ii++) { index[ii] = ii; }
1952  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1953  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1954  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1955  }
1956 
1957  // Form HP = H*P
1958  Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1959  {
1960  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1961  HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1962  }
1963 
1964  // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
1965  int info = 0, lwork = -1;
1966  tau_.resize (keff_new);
1967  lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1968  HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1970  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1971  "LAPACK's _GEQRF failed to compute a workspace size.");
1972 
1973  // NOTE (mfh 18 Apr 2014) LAPACK promises that the value of work_[0]
1974  // after the workspace query will fit in int. This justifies the
1975  // cast. We call real() first because static_cast from std::complex
1976  // to int doesn't work.
1977  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1978  work_.resize (lwork); // Allocate workspace for the QR factorization
1979  lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1980  HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1982  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1983  "LAPACK's _GEQRF failed to compute a QR factorization.");
1984 
1985  // Explicitly construct Q and R factors
1986  // NOTE: The upper triangular part of HP is copied into R and HP becomes Q.
1987  Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
1988  for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
1989 
1990  // NOTE (mfh 18 Apr 2014): Teuchos::LAPACK's wrapper for UNGQR
1991  // dispatches to the correct Scalar-specific routine. It calls
1992  // {S,D}ORGQR if Scalar is real, and {C,Z}UNGQR if Scalar is
1993  // complex.
1994  lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1995  HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1996  lwork, &info);
1998  info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1999  "LAPACK's _UNGQR failed to construct the Q factor.");
2000 
2001  // Form orthonormalized C and adjust U accordingly so that C = A*U
2002  // C = [C V] * Q;
2003 
2004  // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
2005  {
2006  Teuchos::RCP<MV> C1tmp;
2007  {
2008  index.resize(keff);
2009  for (int i=0; i < keff; i++) { index[i] = i; }
2010  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2011  index.resize(keff_new);
2012  for (int i=0; i < keff_new; i++) { index[i] = i; }
2013  C1tmp = MVT::CloneViewNonConst( *C1_, index );
2014  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2015  MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2016  }
2017  // Now compute C += V(:,1:p+1) * Q
2018  {
2019  index.resize( p+1 );
2020  for (int i=0; i < p+1; ++i) { index[i] = i; }
2021  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2022  Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2023  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2024  }
2025  }
2026 
2027  // C_ = C1_; (via a swap)
2028  std::swap(C_, C1_);
2029 
2030  // Finally, compute U_ = U_*R^{-1}
2031  // First, compute LU factorization of R
2032  ipiv_.resize(Rtmp.numRows());
2033  lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2034  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2035 
2036  // Now, form inv(R)
2037  lwork = Rtmp.numRows();
2038  work_.resize(lwork);
2039  lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2040  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2041 
2042  {
2043  index.resize(keff_new);
2044  for (int i=0; i < keff_new; i++) { index[i] = i; }
2045  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2046  MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2047  }
2048 
2049  // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
2050  if (keff != keff_new) {
2051  keff = keff_new;
2052  gcrodr_iter->setSize( keff, numBlocks_ );
2053  // Important to zero this out before next cyle
2054  Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2055  b1.putScalar(zero);
2056  }
2057 
2058 }
2059 
2060 
2061 // Compute the harmonic eigenpairs of the projected, dense system.
2062 template<class ScalarType, class MV, class OP>
2063 int GCRODRSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs1(int m,
2066  int i, j;
2067  bool xtraVec = false;
2068  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2069 
2070  // Real and imaginary eigenvalue components
2071  std::vector<MagnitudeType> wr(m), wi(m);
2072 
2073  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
2075 
2076  // Magnitude of harmonic Ritz values
2077  std::vector<MagnitudeType> w(m);
2078 
2079  // Sorted order of harmonic Ritz values, also used for GEEV
2080  std::vector<int> iperm(m);
2081 
2082  // Output info
2083  int info = 0;
2084 
2085  // Set flag indicating recycle space has been generated this solve
2086  builtRecycleSpace_ = true;
2087 
2088  // Solve linear system: H_m^{-H}*e_m
2091  e_m[m-1] = one;
2092  lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2093  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2094 
2095  // Compute H_m + d*H_m^{-H}*e_m*e_m^H
2096  ScalarType d = HH(m, m-1) * HH(m, m-1);
2098  for( i=0; i<m; ++i )
2099  harmHH(i, m-1) += d * e_m[i];
2100 
2101  // Revise to do query for optimal workspace first
2102  // Create simple storage for the left eigenvectors, which we don't care about.
2103  const int ldvl = 1;
2104  ScalarType* vl = 0;
2105 
2106  // Size of workspace and workspace for GEEV
2107  int lwork = -1;
2108  std::vector<ScalarType> work(1);
2109  std::vector<MagnitudeType> rwork(2*m);
2110 
2111  // First query GEEV for the optimal workspace size
2112  lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2113  vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2114 
2115  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
2116  work.resize( lwork );
2117 
2118  lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2119  vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2120  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2121 
2122  // Construct magnitude of each harmonic Ritz value
2123  for( i=0; i<m; ++i )
2124  w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2125 
2126  // Construct magnitude of each harmonic Ritz value
2127  this->sort(w, m, iperm);
2128 
2129  const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2130 
2131  // Select recycledBlocks_ smallest eigenvectors
2132  for( i=0; i<recycledBlocks_; ++i ) {
2133  for( j=0; j<m; j++ ) {
2134  PP(j,i) = vr(j,iperm[i]);
2135  }
2136  }
2137 
2138  if(!scalarTypeIsComplex) {
2139 
2140  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
2141  if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2142  int countImag = 0;
2143  for ( i=0; i<recycledBlocks_; ++i ) {
2144  if (wi[iperm[i]] != 0.0)
2145  countImag++;
2146  }
2147  // Check to see if this count is even or odd:
2148  if (countImag % 2)
2149  xtraVec = true;
2150  }
2151 
2152  if (xtraVec) { // we need to store one more vector
2153  if (wi[iperm[recycledBlocks_-1]] > 0.0) { // I picked the "real" component
2154  for( j=0; j<m; ++j ) { // so get the "imag" component
2155  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2156  }
2157  }
2158  else { // I picked the "imag" component
2159  for( j=0; j<m; ++j ) { // so get the "real" component
2160  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2161  }
2162  }
2163  }
2164 
2165  }
2166 
2167  // Return whether we needed to store an additional vector
2168  if (xtraVec) {
2169  return recycledBlocks_+1;
2170  }
2171  else {
2172  return recycledBlocks_;
2173  }
2174 
2175 }
2176 
2177 // Compute the harmonic eigenpairs of the projected, dense system.
2178 template<class ScalarType, class MV, class OP>
2179 int GCRODRSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs2(int keffloc, int m,
2181  const Teuchos::RCP<const MV>& VV,
2183  int i, j;
2184  int m2 = HH.numCols();
2185  bool xtraVec = false;
2186  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2187  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2188  std::vector<int> index;
2189 
2190  // Real and imaginary eigenvalue components
2191  std::vector<MagnitudeType> wr(m2), wi(m2);
2192 
2193  // Magnitude of harmonic Ritz values
2194  std::vector<MagnitudeType> w(m2);
2195 
2196  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
2198 
2199  // Sorted order of harmonic Ritz values
2200  std::vector<int> iperm(m2);
2201 
2202  // Set flag indicating recycle space has been generated this solve
2203  builtRecycleSpace_ = true;
2204 
2205  // Form matrices for generalized eigenproblem
2206 
2207  // B = H2' * H2; Don't zero out matrix when constructing
2209  B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
2210 
2211  // A_tmp = | C'*U 0 |
2212  // | V_{m+1}'*U I |
2213  Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2214 
2215  // A_tmp(1:keffloc,1:keffloc) = C' * U;
2216  index.resize(keffloc);
2217  for (i=0; i<keffloc; ++i) { index[i] = i; }
2218  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2219  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2220  Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2221  MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2222 
2223  // A_tmp(keffloc+1:m-k+keffloc+1,1:keffloc) = V' * U;
2224  Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2225  index.resize(m+1);
2226  for (i=0; i < m+1; i++) { index[i] = i; }
2227  Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2228  MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2229 
2230  // A_tmp(keffloc+1:m-k+keffloc,keffloc+1:m-k+keffloc) = eye(m-k);
2231  for( i=keffloc; i<keffloc+m; i++ ) {
2232  A_tmp(i,i) = one;
2233  }
2234 
2235  // A = H2' * A_tmp;
2236  Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2237  A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2238 
2239  // Compute k smallest harmonic Ritz pairs
2240  // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
2241  // ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
2242  // IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
2243  // RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
2244  // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
2245  char balanc='P', jobvl='N', jobvr='V', sense='N';
2246  int ld = A.numRows();
2247  int lwork = 6*ld;
2248  int ldvl = ld, ldvr = ld;
2249  int info = 0,ilo = 0,ihi = 0;
2250  MagnitudeType abnrm = 0.0, bbnrm = 0.0;
2251  ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
2252  std::vector<ScalarType> beta(ld);
2253  std::vector<ScalarType> work(lwork);
2254  std::vector<MagnitudeType> rwork(lwork);
2255  std::vector<MagnitudeType> lscale(ld), rscale(ld);
2256  std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2257  std::vector<int> iwork(ld+6);
2258  int *bwork = 0; // If sense == 'N', bwork is never referenced
2259  //lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2260  // &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2261  // &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
2262  lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2263  &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2264  &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2265  &iwork[0], bwork, &info);
2266  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2267 
2268  // Construct magnitude of each harmonic Ritz value
2269  // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
2270  for( i=0; i<ld; i++ ) {
2271  w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2273  }
2274 
2275  // Construct magnitude of each harmonic Ritz value
2276  this->sort(w,ld,iperm);
2277 
2278  const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2279 
2280  // Select recycledBlocks_ smallest eigenvectors
2281  for( i=0; i<recycledBlocks_; i++ ) {
2282  for( j=0; j<ld; j++ ) {
2283  PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2284  }
2285  }
2286 
2287  if(!scalarTypeIsComplex) {
2288 
2289  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
2290  if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2291  int countImag = 0;
2292  for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2293  if (wi[iperm[i]] != 0.0)
2294  countImag++;
2295  }
2296  // Check to see if this count is even or odd:
2297  if (countImag % 2)
2298  xtraVec = true;
2299  }
2300 
2301  if (xtraVec) { // we need to store one more vector
2302  if (wi[iperm[ld-recycledBlocks_]] > 0.0) { // I picked the "real" component
2303  for( j=0; j<ld; j++ ) { // so get the "imag" component
2304  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2305  }
2306  }
2307  else { // I picked the "imag" component
2308  for( j=0; j<ld; j++ ) { // so get the "real" component
2309  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2310  }
2311  }
2312  }
2313 
2314  }
2315 
2316  // Return whether we needed to store an additional vector
2317  if (xtraVec) {
2318  return recycledBlocks_+1;
2319  }
2320  else {
2321  return recycledBlocks_;
2322  }
2323 
2324 }
2325 
2326 
2327 // This method sorts list of n floating-point numbers and return permutation vector
2328 template<class ScalarType, class MV, class OP>
2329 void GCRODRSolMgr<ScalarType,MV,OP,true>::sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm) {
2330  int l, r, j, i, flag;
2331  int RR2;
2332  MagnitudeType dRR, dK;
2333 
2334  // Initialize the permutation vector.
2335  for(j=0;j<n;j++)
2336  iperm[j] = j;
2337 
2338  if (n <= 1) return;
2339 
2340  l = n / 2 + 1;
2341  r = n - 1;
2342  l = l - 1;
2343  dRR = dlist[l - 1];
2344  dK = dlist[l - 1];
2345 
2346  RR2 = iperm[l - 1];
2347  while (r != 0) {
2348  j = l;
2349  flag = 1;
2350 
2351  while (flag == 1) {
2352  i = j;
2353  j = j + j;
2354 
2355  if (j > r + 1)
2356  flag = 0;
2357  else {
2358  if (j < r + 1)
2359  if (dlist[j] > dlist[j - 1]) j = j + 1;
2360 
2361  if (dlist[j - 1] > dK) {
2362  dlist[i - 1] = dlist[j - 1];
2363  iperm[i - 1] = iperm[j - 1];
2364  }
2365  else {
2366  flag = 0;
2367  }
2368  }
2369  }
2370  dlist[i - 1] = dRR;
2371  iperm[i - 1] = RR2;
2372 
2373  if (l == 1) {
2374  dRR = dlist [r];
2375  RR2 = iperm[r];
2376  dK = dlist[r];
2377  dlist[r] = dlist[0];
2378  iperm[r] = iperm[0];
2379  r = r - 1;
2380  }
2381  else {
2382  l = l - 1;
2383  dRR = dlist[l - 1];
2384  RR2 = iperm[l - 1];
2385  dK = dlist[l - 1];
2386  }
2387  }
2388  dlist[0] = dRR;
2389  iperm[0] = RR2;
2390 }
2391 
2392 
2393 template<class ScalarType, class MV, class OP>
2395  std::ostringstream out;
2396  out << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
2397  out << "{";
2398  out << "Ortho Type: \"" << orthoType_ << "\"";
2399  out << ", Num Blocks: " <<numBlocks_;
2400  out << ", Num Recycle Blocks: " << recycledBlocks_;
2401  out << ", Max Restarts: " << maxRestarts_;
2402  out << "}";
2403  return out.str ();
2404 }
2405 
2406 } // namespace Belos
2407 
2408 #endif /* BELOS_GCRODR_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:74
ScalarType * values() const
Collection of types and exceptions used within the Belos solvers.
Partial specialization for ScalarType types for which Teuchos::LAPACK has a valid implementation...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
GCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:88
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
T & get(const std::string &name, T def_value)
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call...
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
static T squareroot(T x)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
An implementation of StatusTestResNorm using a family of residual norms.
GCRODRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
GCRODRSolMgrRecyclingFailure(const std::string &what_arg)
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
Belos::StatusTest class for specifying a maximum number of iterations.
int curDim
The current dimension of the reduction.
static std::string name()
A factory class for generating StatusTestOutput objects.
std::string validNamesString() const
List (as a string) of recognized MatOrthoManager names.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
A Belos::StatusTest class for specifying a maximum number of iterations.
GCRODRSolMgrOrthoFailure(const std::string &what_arg)
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:174
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::RCP< MV > V
The current Krylov basis.
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)
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
A linear system to solve, and its associated information.
GCRODRSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorm...
Class which describes the linear problem to be solved by the iterative solver.
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:123
Belos concrete class for performing the GCRO-DR iteration.
Structure to contain pointers to GCRODRIter state variables.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
GCRODRSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspa...
static magnitudeType magnitude(T a)
OrdinalType numCols() const
GCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
GCRODRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
bool isType(const std::string &name) const
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
A class for extending the status testing capabilities of Belos via logical combinations.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Class which defines basic traits for the operator type.
Teuchos::RCP< MV > C
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
OrdinalType stride() const
OrdinalType numRows() const
Belos concrete class for performing the block, flexible GMRES iteration.
bool is_null() const
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.

Generated on Wed Jul 30 2025 09:23:52 for Belos by doxygen 1.8.5