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

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11