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