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