Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBlockGCRODRSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP
47 
48 #include "BelosConfigDefs.hpp"
49 #include "BelosTypes.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosSolverManager.hpp"
53 #include "BelosGmresIteration.hpp"
54 #include "BelosBlockGCRODRIter.hpp"
55 #include "BelosBlockGmresIter.hpp"
56 #include "BelosBlockFGmresIter.hpp"
59 #include "BelosStatusTestCombo.hpp"
61 #include "BelosOutputManager.hpp"
62 #include "Teuchos_LAPACK.hpp"
63 #include "Teuchos_as.hpp"
64 #ifdef BELOS_TEUCHOS_TIME_MONITOR
65 #include "Teuchos_TimeMonitor.hpp"
66 #endif // BELOS_TEUCHOS_TIME_MONITOR
67 
68 // MLP: Add links to examples here
69 
70 namespace Belos{
71 
73 
74 
80  public:
81  BlockGCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
82  };
83 
89  public:
90  BlockGCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
91  };
92 
98  public:
99  BlockGCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
100  };
101 
107  public:
108  BlockGCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
109  };
110 
112 
124 
125 template<class ScalarType, class MV, class OP>
126 class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP> {
127 private:
128 
138 
139 public:
141 
142 
149 
177 
179  virtual ~BlockGCRODRSolMgr() {};
181 
184 
186  std::string description() const;
187 
189 
190 
192 
193 
196  return *problem_;
197  }
198 
201 
204 
207  return achievedTol_;
208  }
209 
211  int getNumIters() const { return numIters_; }
212 
214  bool isLOADetected() const { return loaDetected_; }
215 
217 
219 
220 
223  TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
224  "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
225 
226  // Check state of problem object before proceeding
227  if (! problem->isProblemSet()) {
228  const bool success = problem->setProblem();
229  TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
230  "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the "
231  "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
232  }
233 
234  problem_ = problem;
235  }
236 
239 
241 
243 
244 
251  void reset (const ResetType type) {
252  if ((type & Belos::Problem) && ! problem_.is_null ())
253  problem_->setProblem();
254  else if (type & Belos::RecycleSubspace)
255  keff = 0;
256  }
257 
259 
261 
262 
278  // \returns ::ReturnType specifying:
281  ReturnType solve();
282 
284 
285 private:
286 
287  // Called by all constructors; Contains init instructions common to all constructors
288  void init();
289 
290  // Initialize solver state storage
291  void initializeStateStorage();
292 
293  // Recycling Methods
294  // Appending Function name by:
295  // "Kryl" indicates it is specialized for building a recycle space after an
296  // initial run of Block GMRES which generates an initial unaugmented block Krylov subspace
297  // "AugKryl" indicates it is specialized for building a recycle space from the augmented Krylov subspace
298 
299  // Functions which control the building of a recycle space
302 
303  // Recycling with Harmonic Ritz Vectors
304  // Computes harmonic eigenpairs of projected matrix created during the priming solve.
305  // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
306 
307  // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension (numBlocks+1)*blockSize x numBlocks.
308  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
309  int getHarmonicVecsKryl (int m, const SDM& HH, SDM& PP);
310 
311  // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+(numBlocks+1)*blockSize x keff+numBlocksm.
312  // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) (numBlocks+1)*blockSize vectors.
313  // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
314  int getHarmonicVecsAugKryl (int keff,
315  int m,
316  const SDM& HH,
317  const Teuchos::RCP<const MV>& VV,
318  SDM& PP);
319 
320  // Sort list of n floating-point numbers and return permutation vector
321  void sort (std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm);
322 
323  // Lapack interface
325 
328 
329  //Output Manager
332 
333  //Status Test
339 
342 
349 
352 
364 
365  //Default Solver Values
366  static const bool adaptiveBlockSize_default_;
367  static const std::string recycleMethod_default_;
368 
369  //Current Solver Values
376  std::string label_;
377 
379  // Solver State Storage
381  //
382  // The number of blocks and recycle blocks (m and k, respectively)
384  // Current size of recycled subspace
385  int keff;
386  //
387  // Residual Vector
389  //
390  // Search Space
392  //
393  // Recycle subspace and its image under action of the operator
395  //
396  // Updated recycle Space and its image under action of the operator
398  //
399  // Storage used in constructing recycle space
405  std::vector<ScalarType> tau_;
406  std::vector<ScalarType> work_;
408  std::vector<int> ipiv_;
409 
412 
414  bool isSet_;
415 
418 
421 };
422 
423  //
424  // Set default solver values
425  //
426  template<class ScalarType, class MV, class OP>
428 
429  template<class ScalarType, class MV, class OP>
431 
432  //
433  // Method definitions
434  //
435 
436  template<class ScalarType, class MV, class OP>
438  init();
439  }
440 
441  //Basic Constructor
442  template<class ScalarType, class MV, class OP>
446  // Initialize local pointers to null, and initialize local
447  // variables to default values.
448  init ();
449 
450  TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
451  "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
452  "the linear problem argument 'problem' to be nonnull.");
453 
454  problem_ = problem;
455 
456  // Set the parameters using the list that was passed in. If null, we defer initialization until a non-null list is set (by the
457  // client calling setParameters(), or by calling solve() -- in either case, a null parameter list indicates that default parameters should be used).
458  if (! pl.is_null())
459  setParameters (pl);
460  }
461 
462  template<class ScalarType, class MV, class OP>
464  adaptiveBlockSize_ = adaptiveBlockSize_default_;
465  recycleMethod_ = recycleMethod_default_;
466  isSet_ = false;
467  loaDetected_ = false;
468  builtRecycleSpace_ = false;
469  keff = 0;//Effective Size of Recycle Space
470  //The following are all RCP smart pointers to indicated matrices/vectors.
471  //Some MATLAB notation used in comments.
472  R_ = Teuchos::null;//The Block Residual
473  V_ = Teuchos::null;//Block Arnoldi Vectors
474  U_ = Teuchos::null;//Recycle Space
475  C_ = Teuchos::null;//Image of U Under Action of Operator
476  U1_ = Teuchos::null;//Newly Computed Recycle Space
477  C1_ = Teuchos::null;//Image of Newly Computed Recycle Space
478  PP_ = Teuchos::null;//Coordinates of New Recyc. Vectors in Augmented Space
479  HP_ = Teuchos::null;//H_*PP_ or G_*PP_
480  G_ = Teuchos::null;//G_ such that A*[U V(:,1:numBlocks_*blockSize_)] = [C V_]*G_
481  F_ = Teuchos::null;//Upper Triangular portion of QR factorization for HP_
482  H_ = Teuchos::null;//H_ such that A*V(:,1:numBlocks_*blockSize_) = V_*H_ + C_*C_'*V_
483  B_ = Teuchos::null;//B_ = C_'*V_
484 
485  //THIS BLOCK OF CODE IS COMMENTED OUT AND PLACED ELSEWHERE IN THE CODE
486 /*
487  //WE TEMPORARILY INITIALIZE STATUS TESTS HERE FOR TESTING PURPOSES, BUT THEY SHOULD BE
488  //INITIALIZED SOMEWHERE ELSE, LIKE THE setParameters() FUNCTION
489 
490  //INSTANTIATE AND INITIALIZE TEST OBJECTS AS NEEDED
491  if (maxIterTest_.is_null()){
492  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
493  }
494  //maxIterTest_->setMaxIters (maxIters_);
495 
496  //INSTANTIATE THE PRINTER
497  if (printer_.is_null()) {
498  printer_ = Teuchos::rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
499  }
500 
501  if (ortho_.is_null()) // || changedOrthoType) %%%In other codes, this is also triggered if orthogonalization type changed {
502  // Create orthogonalization manager. This requires that the OutputManager (printer_) already be initialized.
503  Teuchos::RCP<const Teuchos::ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
504  ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, orthoParams);
505  }
506 
507  // Convenience typedefs
508  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
509  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
510 
511  if (impConvTest_.is_null()) {
512  impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
513  impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), Belos::TwoNorm);
514  impConvTest_->setShowMaxResNormOnly( true );
515  }
516 
517  if (expConvTest_.is_null()) {
518  expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
519  expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
520  expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), Belos::TwoNorm);
521  expConvTest_->setShowMaxResNormOnly( true );
522  }
523 
524  if (convTest_.is_null()) {
525  convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ, impConvTest_, expConvTest_));
526  }
527 
528  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
529 
530  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
531  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
532  Passed+Failed+Undefined);
533 
534  */
535 
536  }
537 
538  // This method requires the solver manager to return a string that describes itself.
539  template<class ScalarType, class MV, class OP>
541  std::ostringstream oss;
542  oss << "Belos::BlockGCRODRSolMgr<" << SCT::name() << ", ...>";
543  oss << "{";
544  oss << "Ortho Type='"<<orthoType_ ;
545  oss << ", Num Blocks=" <<numBlocks_;
546  oss << ", Block Size=" <<blockSize_;
547  oss << ", Num Recycle Blocks=" << recycledBlocks_;
548  oss << ", Max Restarts=" << maxRestarts_;
549  oss << "}";
550  return oss.str();
551  }
552 
553  template<class ScalarType, class MV, class OP>
557  using Teuchos::parameterList;
558  using Teuchos::RCP;
559  using Teuchos::rcpFromRef;
560 
561  if (defaultParams_.is_null()) {
562  RCP<ParameterList> pl = parameterList ();
563 
564  const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
565  const int maxRestarts = 1000;
566  const int maxIters = 5000;
567  const int blockSize = 2;
568  const int numBlocks = 100;
569  const int numRecycledBlocks = 25;
571  const Belos::OutputType outputStyle = Belos::General;
572  const int outputFreq = 1;
573  RCP<std::ostream> outputStream = rcpFromRef (std::cout);
574  const std::string impResScale ("Norm of Preconditioned Initial Residual");
575  const std::string expResScale ("Norm of Initial Residual");
576  const std::string timerLabel ("Belos");
577  const std::string orthoType ("ICGS");
578  RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
579  //const MagnitudeType orthoKappa = SCT::magnitude (SCT::eps());
580  //
581  // mfh 31 Dec 2011: Negative values of orthoKappa are ignored.
582  // Setting this to a negative value by default ensures that
583  // this parameter is only _not_ ignored if the user
584  // specifically sets a valid value.
585  const MagnitudeType orthoKappa = -SMT::one();
586 
587  // Set all the valid parameters and their default values.
588  pl->set ("Convergence Tolerance", convTol,
589  "The tolerance that the solver needs to achieve in order for "
590  "the linear system(s) to be declared converged. The meaning "
591  "of this tolerance depends on the convergence test details.");
592  pl->set("Maximum Restarts", maxRestarts,
593  "The maximum number of restart cycles (not counting the first) "
594  "allowed for each set of right-hand sides solved.");
595  pl->set("Maximum Iterations", maxIters,
596  "The maximum number of iterations allowed for each set of "
597  "right-hand sides solved.");
598  pl->set("Block Size", blockSize,
599  "How many linear systems to solve at once.");
600  pl->set("Num Blocks", numBlocks,
601  "The maximum number of blocks allowed in the Krylov subspace "
602  "for each set of right-hand sides solved. This includes the "
603  "initial block. Total Krylov basis storage, not counting the "
604  "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
605  pl->set("Num Recycled Blocks", numRecycledBlocks,
606  "The maximum number of vectors in the recycled subspace." );
607  pl->set("Verbosity", verbosity,
608  "What type(s) of solver information should be written "
609  "to the output stream.");
610  pl->set("Output Style", outputStyle,
611  "What style is used for the solver information to write "
612  "to the output stream.");
613  pl->set("Output Frequency", outputFreq,
614  "How often convergence information should be written "
615  "to the output stream.");
616  pl->set("Output Stream", outputStream,
617  "A reference-counted pointer to the output stream where all "
618  "solver output is sent.");
619  pl->set("Implicit Residual Scaling", impResScale,
620  "The type of scaling used in the implicit residual convergence test.");
621  pl->set("Explicit Residual Scaling", expResScale,
622  "The type of scaling used in the explicit residual convergence test.");
623  pl->set("Timer Label", timerLabel,
624  "The string to use as a prefix for the timer labels.");
625  pl->set("Orthogonalization", orthoType,
626  "The orthogonalization method to use. Valid options: " +
627  orthoFactory_.validNamesString());
628  pl->set ("Orthogonalization Parameters", *orthoParams,
629  "Sublist of parameters specific to the orthogonalization method to use.");
630  pl->set("Orthogonalization Constant", orthoKappa,
631  "When using DGKS orthogonalization: the \"depTol\" constant, used "
632  "to determine whether another step of classical Gram-Schmidt is "
633  "necessary. Otherwise ignored. Nonpositive values are ignored.");
634  defaultParams_ = pl;
635  }
636  return defaultParams_;
637  }
638 
639  template<class ScalarType, class MV, class OP>
640  void
643  using Teuchos::isParameterType;
644  using Teuchos::getParameter;
645  using Teuchos::null;
647  using Teuchos::parameterList;
648  using Teuchos::RCP;
649  using Teuchos::rcp;
650  using Teuchos::rcp_dynamic_cast;
651  using Teuchos::rcpFromRef;
655 
656  // The default parameter list contains all parameters that BlockGCRODRSolMgr understands, and none that it doesn't
657  RCP<const ParameterList> defaultParams = getValidParameters();
658 
659  if (params.is_null()) {
660  // Users really should supply a non-null input ParameterList,
661  // so that we can fill it in. However, if they really did give
662  // us a null list, be nice and use default parameters. In this
663  // case, we don't have to do any validation.
664  params_ = parameterList (*defaultParams);
665  }
666  else {
667  // Validate the user's parameter list and fill in any missing parameters with default values.
668  //
669  // Currently, validation is quite strict. The following line
670  // will throw an exception for misspelled or extra parameters.
671  // However, it will _not_ throw an exception if there are
672  // missing parameters: these will be filled in with their
673  // default values. There is additional discussion of other
674  // validation strategies in the comments of this function for
675  // Belos::GCRODRSolMgr.
676  params->validateParametersAndSetDefaults (*defaultParams);
677  // No side effects on the solver until after validation passes.
678  params_ = params;
679  }
680 
681  //
682  // Validate and set parameters relating to the Krylov subspace
683  // dimensions and the number of blocks.
684  //
685  // We try to read and validate all these parameters' values
686  // before setting them in the solver. This is because the
687  // validity of each may depend on the others' values. Our goal
688  // is to avoid incorrect side effects, so we don't set the values
689  // in the solver until we know they are right.
690  //
691  {
692  const int maxRestarts = params_->get<int> ("Maximum Restarts");
693  TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
694  "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
695  "must be nonnegative, but you specified a negative value of "
696  << maxRestarts << ".");
697 
698  const int maxIters = params_->get<int> ("Maximum Iterations");
699  TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
700  "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
701  "must be positive, but you specified a nonpositive value of "
702  << maxIters << ".");
703 
704  const int numBlocks = params_->get<int> ("Num Blocks");
705  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
706  "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
707  "positive, but you specified a nonpositive value of " << numBlocks
708  << ".");
709 
710  const int blockSize = params_->get<int> ("Block Size");
711  TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
712  "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
713  "positive, but you specified a nonpositive value of " << blockSize
714  << ".");
715 
716  const int recycledBlocks = params_->get<int> ("Num Recycled Blocks");
717  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
718  "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
719  "be positive, but you specified a nonpositive value of "
720  << recycledBlocks << ".");
721  TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
722  std::invalid_argument, "Belos::BlockGCRODRSolMgr: The \"Num Recycled "
723  "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
724  "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
725  << " and \"Num Blocks\" = " << numBlocks << ".");
726 
727  // Now that we've validated the various dimension-related
728  // parameters, we can set them in the solvers.
729  maxRestarts_ = maxRestarts;
730  maxIters_ = maxIters;
731  numBlocks_ = numBlocks;
732  blockSize_ = blockSize;
733  recycledBlocks_ = recycledBlocks;
734  }
735 
736  // Check to see if the timer label changed. If it did, update it in
737  // the parameter list, and create a new timer with that label (if
738  // Belos was compiled with timers enabled).
739  {
740  std::string tempLabel = params_->get<std::string> ("Timer Label");
741  const bool labelChanged = (tempLabel != label_);
742 
743 #ifdef BELOS_TEUCHOS_TIME_MONITOR
744  std::string solveLabel = label_ + ": BlockGCRODRSolMgr total solve time";
745  if (timerSolve_.is_null()) {
746  // We haven't created a timer yet.
747  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
748  } else if (labelChanged) {
749  // We created a timer before with a different label. In that
750  // case, clear the old timer and create a new timer with the
751  // new label. Clearing old timers prevents them from piling
752  // up, since they are stored statically, possibly without
753  // checking for duplicates.
755  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
756  }
757 #endif // BELOS_TEUCHOS_TIME_MONITOR
758  }
759 
760  // Check for a change in verbosity level. Just in case, we allow
761  // the type of this parameter to be either int or MsgType (an
762  // enum).
763  if (params_->isParameter ("Verbosity")) {
764  if (isParameterType<int> (*params_, "Verbosity")) {
765  verbosity_ = params_->get<int> ("Verbosity");
766  }
767  else {
768  verbosity_ = (int) getParameter<MsgType> (*params_, "Verbosity");
769  }
770  }
771 
772  // Check for a change in output style. Just in case, we allow
773  // the type of this parameter to be either int or OutputType (an
774  // enum).
775  if (params_->isParameter ("Output Style")) {
776  if (isParameterType<int> (*params_, "Output Style")) {
777  outputStyle_ = params_->get<int> ("Output Style");
778  }
779  else {
780  outputStyle_ = (int) getParameter<OutputType> (*params_, "Output Style");
781  }
782 
783  // Currently, the only way we can update the output style is to
784  // create a new output test. This is because we must first
785  // instantiate a new StatusTestOutputFactory with the new
786  // style, and then use the factory to make a new output test.
787  // Thus, we set outputTest_ to null for now, so we remember
788  // later to create a new one.
789  outputTest_ = null;
790  }
791 
792  // Get the output stream for the output manager.
793  //
794  // It has been pointed out (mfh 28 Feb 2011 in GCRODRSolMgr code)
795  // that including an RCP<std::ostream> parameter makes it
796  // impossible to serialize the parameter list. If you serialize
797  // the list and read it back in from the serialized
798  // representation, you might not even get a valid
799  // RCP<std::ostream>, let alone the same output stream that you
800  // serialized.
801  //
802  // However, existing Belos users depend on setting the output
803  // stream in the parameter list. We retain this behavior for
804  // backwards compatibility.
805  //
806  // In case the output stream can't be read back in, we default to
807  // stdout (std::cout), just to ensure reasonable behavior.
808  if (params_->isParameter ("Output Stream")) {
809  try {
810  outputStream_ = getParameter<RCP<std::ostream> > (*params_, "Output Stream");
811  }
812  catch (InvalidParameter&) {
813  outputStream_ = rcpFromRef (std::cout);
814  }
815  // We assume that a null output stream indicates that the user
816  // doesn't want to print anything, so we replace it with a
817  // "black hole" stream that prints nothing sent to it. (We
818  // can't use outputStream_ = Teuchos::null, since the output
819  // manager assumes that operator<< works on its output stream.)
820  if (outputStream_.is_null()) {
821  outputStream_ = rcp (new Teuchos::oblackholestream);
822  }
823  }
824 
825  outputFreq_ = params_->get<int> ("Output Frequency");
826  if (verbosity_ & Belos::StatusTestDetails) {
827  // Update parameter in our output status test.
828  if (! outputTest_.is_null()) {
829  outputTest_->setOutputFrequency (outputFreq_);
830  }
831  }
832 
833  // Create output manager if we need to, using the verbosity level
834  // and output stream that we fetched above. We do this here because
835  // instantiating an OrthoManager using OrthoManagerFactory requires
836  // a valid OutputManager.
837  if (printer_.is_null()) {
838  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
839  } else {
840  printer_->setVerbosity (verbosity_);
841  printer_->setOStream (outputStream_);
842  }
843 
844  // Get the orthogonalization manager name ("Orthogonalization").
845  //
846  // Getting default values for the orthogonalization manager
847  // parameters ("Orthogonalization Parameters") requires knowing the
848  // orthogonalization manager name. Save it for later, and also
849  // record whether it's different than before.
850 
851  //bool changedOrthoType = false; // silence warning about not being used
852  if (params_->isParameter ("Orthogonalization")) {
853  const std::string& tempOrthoType =
854  params_->get<std::string> ("Orthogonalization");
855  // Ensure that the specified orthogonalization type is valid.
856  if (! orthoFactory_.isValidName (tempOrthoType)) {
857  std::ostringstream os;
858  os << "Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
859  << tempOrthoType << "\". The following are valid options "
860  << "for the \"Orthogonalization\" name parameter: ";
861  orthoFactory_.printValidNames (os);
862  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
863  }
864  if (tempOrthoType != orthoType_) {
865  //changedOrthoType = true;
866  orthoType_ = tempOrthoType;
867  }
868  }
869 
870  // Get any parameters for the orthogonalization
871  // ("Orthogonalization Parameters"). If not supplied, the
872  // orthogonalization manager factory will supply default values.
873  //
874  // NOTE (mfh 12 Jan 2011, summer 2011, 31 Dec 2011) For backwards
875  // compatibility with other Belos GMRES-type solvers, if params
876  // has an "Orthogonalization Constant" parameter and the DGKS
877  // orthogonalization manager is to be used, the value of this
878  // parameter will override DGKS's "depTol" parameter.
879  //
880  // Users must supply the orthogonalization manager parameters as
881  // a sublist (supplying it as an RCP<ParameterList> would make
882  // the resulting parameter list not serializable).
883  RCP<ParameterList> orthoParams = sublist (params, "Orthogonalization Parameters", true);
884  TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
885  "Failed to get orthogonalization parameters. "
886  "Please report this bug to the Belos developers.");
887 
888  // Check if the desired orthogonalization method changed, or if
889  // the orthogonalization manager has not yet been instantiated.
890  // If either is the case, instantiate a new MatOrthoManager
891  // subclass instance corresponding to the desired
892  // orthogonalization method. We've already fetched the
893  // orthogonalization method name (orthoType_) and its parameters
894  // (orthoParams) above.
895  //
896  // As suggested (by mfh 12 Jan 2011 in Belos::GCRODRSolMgr) In
897  // order to ensure parameter changes get propagated to the
898  // orthomanager we simply reinstantiate the OrthoManager every
899  // time, whether or not the orthogonalization method name or
900  // parameters have changed. This is not efficient for some
901  // orthogonalization managers that keep a lot of internal state.
902  // A more general way to fix this inefficiency would be to
903  //
904  // 1. make each orthogonalization implement
905  // Teuchos::ParameterListAcceptor, and
906  //
907  // 2. make each orthogonalization implement a way to set the
908  // OutputManager instance and timer label.
909  //
910  // Create orthogonalization manager. This requires that the
911  // OutputManager (printer_) already be initialized.
912  ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
913  label_, orthoParams);
914 
915  // The DGKS orthogonalization accepts a "Orthogonalization
916  // Constant" parameter (also called kappa in the code, but not in
917  // the parameter list). If its value is provided in the given
918  // parameter list and its value is positive, use it. Ignore
919  // negative values.
920  //
921  // The "Orthogonalization Constant" parameter is retained only
922  // for backwards compatibility.
923 
924  //bool gotValidOrthoKappa = false; // silence warning about not being used
925  if (params_->isParameter ("Orthogonalization Constant")) {
926  const MagnitudeType orthoKappa =
927  params_->get<MagnitudeType> ("Orthogonalization Constant");
928  if (orthoKappa > 0) {
929  orthoKappa_ = orthoKappa;
930  //gotValidOrthoKappa = true;
931  // Only DGKS currently accepts this parameter.
932  if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
933  typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
934  // This cast should always succeed; it's a bug otherwise.
935  // (If the cast fails, then orthoType_ doesn't correspond
936  // to the OrthoManager subclass instance that we think we
937  // have, so we initialized the wrong subclass somehow.)
938  rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
939  }
940  }
941  }
942 
943  // Convergence
944  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
945  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
946 
947  // Check for convergence tolerance
948  convTol_ = params_->get<MagnitudeType> ("Convergence Tolerance");
949  if (! impConvTest_.is_null()) {
950  impConvTest_->setTolerance (convTol_);
951  }
952  if (! expConvTest_.is_null()) {
953  expConvTest_->setTolerance (convTol_);
954  }
955 
956  // Check for a change in scaling, if so we need to build new residual tests.
957  if (params_->isParameter ("Implicit Residual Scaling")) {
958  std::string tempImpResScale =
959  getParameter<std::string> (*params_, "Implicit Residual Scaling");
960 
961  // Only update the scaling if it's different.
962  if (impResScale_ != tempImpResScale) {
963  ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
964  impResScale_ = tempImpResScale;
965 
966  if (! impConvTest_.is_null()) {
967  try {
968  impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
969  }
970  catch (StatusTestError&) {
971  // Delete the convergence test so it gets constructed again.
972  impConvTest_ = null;
973  convTest_ = null;
974  }
975  }
976  }
977  }
978 
979  if (params_->isParameter("Explicit Residual Scaling")) {
980  std::string tempExpResScale =
981  getParameter<std::string> (*params_, "Explicit Residual Scaling");
982 
983  // Only update the scaling if it's different.
984  if (expResScale_ != tempExpResScale) {
985  ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
986  expResScale_ = tempExpResScale;
987 
988  if (! expConvTest_.is_null()) {
989  try {
990  expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
991  }
992  catch (StatusTestError&) {
993  // Delete the convergence test so it gets constructed again.
994  expConvTest_ = null;
995  convTest_ = null;
996  }
997  }
998  }
999  }
1000 
1001  //
1002  // Create iteration stopping criteria ("status tests") if we need
1003  // to, by combining three different stopping criteria.
1004  //
1005  // First, construct maximum-number-of-iterations stopping criterion.
1006  if (maxIterTest_.is_null()) {
1007  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
1008  } else {
1009  maxIterTest_->setMaxIters (maxIters_);
1010  }
1011 
1012  // Implicit residual test, using the native residual to determine if
1013  // convergence was achieved.
1014  if (impConvTest_.is_null()) {
1015  impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1016  impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
1017  Belos::TwoNorm);
1018  }
1019 
1020  // Explicit residual test once the native residual is below the tolerance
1021  if (expConvTest_.is_null()) {
1022  expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1023  expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
1024  expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
1025  Belos::TwoNorm);
1026  }
1027  // Convergence test first tests the implicit residual, then the
1028  // explicit residual if the implicit residual test passes.
1029  if (convTest_.is_null()) {
1030  convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1031  impConvTest_,
1032  expConvTest_));
1033  }
1034  // Construct the complete iteration stopping criterion:
1035  //
1036  // "Stop iterating if the maximum number of iterations has been
1037  // reached, or if the convergence test passes."
1038  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
1039  maxIterTest_, convTest_));
1040  // Create the status test output class.
1041  // This class manages and formats the output from the status test.
1042  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
1043  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
1045 
1046  // Set the solver string for the output test
1047  std::string solverDesc = "Block GCRODR ";
1048  outputTest_->setSolverDesc (solverDesc);
1049 
1050  // Inform the solver manager that the current parameters were set.
1051  isSet_ = true;
1052  }
1053 
1054  // initializeStateStorage.
1055  template<class ScalarType, class MV, class OP>
1056  void
1058  {
1059 
1060  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1061 
1062  // Check if there is any multivector to clone from.
1063  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1064 
1065  //The Dimension of the Krylov Subspace
1066  int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1067 
1068  //Number of columns in [U_ V_(:,1:KrylSpaDim)]
1069  int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;// + 1 is for possible extra recycle vector
1070 
1071  //Number of columns in [C_ V_]
1072  int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1073 
1074  //TEMPORARY SKELETON DEFINITION OF THIS FUNCTION TO GET THINGS WORKING
1075  //NOT EVERYTHING IS INITIALIZE CORRECTLY YET.
1076 
1077  //INITIALIZE RECYCLE SPACE VARIABLES HERE
1078 
1079  //WE DO NOT ALLOCATE V HERE IN THIS SOLVERMANAGER. If no recycle space exists, then block_gmres_iter
1080  //will allocated V for us. If a recycle space already exists, then we will allocate V after updating the
1081  //recycle space for the current problem.
1082  // If the block Krylov subspace has not been initialized before, generate it using the RHS from lp_.
1083  /*if (V_ == Teuchos::null) {
1084  V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1085  }
1086  else{
1087  // Generate V_ by cloning itself ONLY if more space is needed.
1088  if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1089  Teuchos::RCP<const MV> tmp = V_;
1090  V_ = MVT::Clone( *tmp, numBlocks_+1 );
1091  }
1092  }*/
1093 
1094  //INTITIALIZE SPACE FOR THE NEWLY COMPUTED RECYCLE SPACE VARIABLES HERE
1095 
1096  if (U_ == Teuchos::null) {
1097  U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1098  }
1099  else {
1100  // Generate U_ by cloning itself ONLY if more space is needed.
1101  if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1102  Teuchos::RCP<const MV> tmp = U_;
1103  U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1104  }
1105  }
1106 
1107  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1108  if (C_ == Teuchos::null) {
1109  C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1110  }
1111  else {
1112  // Generate C_ by cloning itself ONLY if more space is needed.
1113  if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1114  Teuchos::RCP<const MV> tmp = C_;
1115  C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1116  }
1117  }
1118 
1119  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1120  if (U1_ == Teuchos::null) {
1121  U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1122  }
1123  else {
1124  // Generate U1_ by cloning itself ONLY if more space is needed.
1125  if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1126  Teuchos::RCP<const MV> tmp = U1_;
1127  U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1128  }
1129  }
1130 
1131  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1132  if (C1_ == Teuchos::null) {
1133  C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1134  }
1135  else {
1136  // Generate C1_ by cloning itself ONLY if more space is needed.
1137  if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1138  Teuchos::RCP<const MV> tmp = C1_;
1139  C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1140  }
1141  }
1142 
1143  // Generate R_ only if it doesn't exist
1144  if (R_ == Teuchos::null){
1145  R_ = MVT::Clone( *rhsMV, blockSize_ );
1146  }
1147 
1148  //INITIALIZE SOME WORK VARIABLES
1149 
1150  // Generate G_ only if it doesn't exist, otherwise resize it.
1151  if (G_ == Teuchos::null){
1152  G_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1153  }
1154  else{
1155  if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1156  {
1157  G_->reshape( augSpaImgDim, augSpaDim );
1158  }
1159  G_->putScalar(zero);
1160  }
1161 
1162  // Generate H_ only if it doesn't exist by pointing it to a view of G_.
1163  if (H_ == Teuchos::null){
1164  H_ = Teuchos::rcp (new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1165  }
1166 
1167  // Generate F_ only if it doesn't exist, otherwise resize it.
1168  if (F_ == Teuchos::null){
1169  F_ = Teuchos::rcp( new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1170  }
1171  else {
1172  if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1173  F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1174  }
1175  }
1176  F_->putScalar(zero);
1177 
1178  // Generate PP_ only if it doesn't exist, otherwise resize it.
1179  if (PP_ == Teuchos::null){
1180  PP_ = Teuchos::rcp( new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1181  }
1182  else{
1183  if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1184  PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1185  }
1186  }
1187 
1188  // Generate HP_ only if it doesn't exist, otherwise resize it.
1189  if (HP_ == Teuchos::null)
1190  HP_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1191  else{
1192  if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1193  HP_->reshape( augSpaImgDim, augSpaDim );
1194  }
1195  }
1196 
1197  // Size of tau_ will change during computation, so just be sure it starts with appropriate size
1198  tau_.resize(recycledBlocks_+1);
1199 
1200  // Size of work_ will change during computation, so just be sure it starts with appropriate size
1201  work_.resize(recycledBlocks_+1);
1202 
1203  // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
1204  ipiv_.resize(recycledBlocks_+1);
1205 
1206  }
1207 
1208 template<class ScalarType, class MV, class OP>
1210 
1211  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1212  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1213 
1214  int p = block_gmres_iter->getState().curDim; //Dimension of the Krylov space generated
1215  std::vector<int> index(keff);//we use this to index certain columns of U, C, and V to
1216  //get views into pieces of these matrices.
1217 
1218  //GET CORRECT PIECE OF MATRIX H APPROPRIATE TO SIZE OF KRYLOV SUBSPACE
1219  SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1220  if(recycledBlocks_ >= p + blockSize_){//keep whole block Krylov subspace
1221  //IF THIS HAS HAPPENED, THIS MEANS WE CONVERGED DURING THIS CYCLE
1222  //THEREFORE, WE DO NOT CONSTRUCT C = A*U;
1223  index.resize(p);
1224  for (int ii=0; ii < p; ++ii) index[ii] = ii;
1225  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1226  MVT::SetBlock(*V_, index, *Utmp);
1227  keff = p;
1228  }
1229  else{ //use a subspace selection method to get recycle space
1230  int info = 0;
1231  Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1232  if(recycleMethod_ == "harmvecs"){
1233  keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1234  printer_->stream(Debug) << "keff = " << keff << std::endl;
1235  }
1236 // Hereafter, only keff columns of PP are needed
1237 PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, keff ) );
1238 // Now get views into C, U, V
1239 index.resize(keff);
1240 for (int ii=0; ii<keff; ++ii) index[ii] = ii;
1241 Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1242 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1243 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1244 index.resize(p);
1245 for (int ii=0; ii < p; ++ii) index[ii] = ii;
1246 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1247 
1248 // Form U (the subspace to recycle)
1249 // U = newstate.V(:,1:p) * PP;
1250 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1251 
1252 // Step #1: Form HP = H*P
1253 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1254 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1255 // Step #1.5: Perform workspace size query for QR factorization of HP
1256 int lwork = -1;
1257 tau_.resize(keff);
1258 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1259 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1260 
1261 // Step #2: Compute QR factorization of HP
1262  lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1263 work_.resize(lwork);
1264 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1265 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1266 
1267 // Step #3: Explicitly construct Q and R factors
1268 // The upper triangular part of HP is copied into R and HP becomes Q.
1269 SDM Rtmp( Teuchos::View, *F_, keff, keff );
1270 for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1271 //lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1272 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1273 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1274  // Now we have [Q,R] = qr(H*P)
1275 
1276  // Now compute C = V(:,1:p+blockSize_) * Q
1277  index.resize( p+blockSize_ );
1278  for (int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1279  Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+blockSize_ vectors now; needed p above)
1280  MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1281 
1282  // Finally, compute U = U*R^{-1}.
1283  // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
1284  // backsolve capabilities don't exist in the Belos::MultiVec class
1285 
1286 // Step #1: First, compute LU factorization of R
1287 ipiv_.resize(Rtmp.numRows());
1288 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1289 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1290 // Step #2: Form inv(R)
1291 lwork = Rtmp.numRows();
1292 work_.resize(lwork);
1293 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1294 //TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1295 // Step #3: Let U = U * R^{-1}
1296 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1297 
1298 } //end else from if(recycledBlocks_ >= p + 1)
1299 return;
1300 } // end buildRecycleSpaceKryl defnition
1301 
1302 template<class ScalarType, class MV, class OP>
1305  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1306 
1307  std::vector<MagnitudeType> d(keff);
1308  std::vector<ScalarType> dscalar(keff);
1309  std::vector<int> index(numBlocks_+1);
1310 
1311  // Get the state
1312  BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1313  int p = oldState.curDim;
1314 
1315  // insufficient new information to update recycle space
1316  if (p<1) return;
1317 
1318  if(recycledBlocks_ >= keff + p){ //we add new Krylov vectors to existing recycle space
1319  // If this has happened, we have converged during this cycle. Therefore, do not construct C = A*C;
1320  index.resize(p);
1321  for (int ii=0; ii < p; ++ii) { index[ii] = keff+ii; } //get a view after current reycle vectors
1322  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1323  for (int ii=0; ii < p; ++ii) { index[ii] =ii; }
1324  MVT::SetBlock(*V_, index, *Utmp);
1325  keff += p;
1326  }
1327 
1328  // Take the norm of the recycled vectors
1329  {
1330  index.resize(keff);
1331  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1332  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1333  d.resize(keff);
1334  dscalar.resize(keff);
1335  MVT::MvNorm( *Utmp, d );
1336  for (int i=0; i<keff; ++i) {
1337  d[i] = one / d[i];
1338  dscalar[i] = (ScalarType)d[i];
1339  }
1340  MVT::MvScale( *Utmp, dscalar );
1341  }
1342 
1343  // Get view into current "full" upper Hessnburg matrix
1344  // note that p describes the dimension of the iter+1 block Krylov space so we have to adjust how we use it to index Gtmp
1345  Teuchos::RCP<SDM> Gtmp = Teuchos::rcp( new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1346 
1347  // Insert D into the leading keff x keff block of G
1348  for (int i=0; i<keff; ++i)
1349  (*Gtmp)(i,i) = d[i];
1350 
1351  // Compute the harmoic Ritz pairs for the generalized eigenproblem
1352  // getHarmonicVecsKryl assumes PP has recycledBlocks_+1 columns available
1353  // See previous block of comments for why we subtract p-blockSize_
1354  int keff_new;
1355  {
1356  SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1357  keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1358  }
1359 
1360  // Code to form new U, C
1361  // U = [U V(:,1:p)] * P; (in two steps)
1362 
1363  // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
1364  Teuchos::RCP<MV> U1tmp;
1365  {
1366  index.resize( keff );
1367  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1368  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1369  index.resize( keff_new );
1370  for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1371  U1tmp = MVT::CloneViewNonConst( *U1_, index );
1372  SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1373  MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1374  }
1375 
1376  // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
1377  {
1378  index.resize(p-blockSize_);
1379  for (int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1380  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1381  SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1382  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1383  }
1384 
1385  // Form GP = G*P
1386  SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1387  {
1388  SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1389  HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1390  }
1391 
1392  // Workspace size query for QR factorization HP= QF (the worksize will be placed in work_[0])
1393  int info = 0, lwork = -1;
1394  tau_.resize(keff_new);
1395  lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1396  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1397 
1398  lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1399  work_.resize(lwork);
1400  lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1401  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1402 
1403  // Explicitly construct Q and F factors
1404  // NOTE: The upper triangular part of HP is copied into F and HP becomes Q.
1405  SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1406  for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1407  //lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1408  lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1409  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1410 
1411  // Form orthonormalized C and adjust U accordingly so that C = A*U
1412  // C = [C V] * Q;
1413 
1414  // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
1415  {
1416  Teuchos::RCP<MV> C1tmp;
1417  {
1418  index.resize(keff);
1419  for (int i=0; i < keff; i++) { index[i] = i; }
1420  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1421  index.resize(keff_new);
1422  for (int i=0; i < keff_new; i++) { index[i] = i; }
1423  C1tmp = MVT::CloneViewNonConst( *C1_, index );
1424  SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1425  MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1426  }
1427  // Now compute C += V(:,1:p+1) * Q
1428  {
1429  index.resize( p );
1430  for (int i=0; i < p; ++i) { index[i] = i; }
1431  Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1432  SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1433  MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1434  }
1435  }
1436 
1437  // C_ = C1_; (via a swap)
1438  std::swap(C_, C1_);
1439 
1440  // Finally, compute U_ = U_*R^{-1}
1441  // First, compute LU factorization of R
1442  ipiv_.resize(Ftmp.numRows());
1443  lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1444  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1445 
1446  // Now, form inv(R)
1447  lwork = Ftmp.numRows();
1448  work_.resize(lwork);
1449  lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1450  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1451 
1452  {
1453  index.resize(keff_new);
1454  for (int i=0; i < keff_new; i++) { index[i] = i; }
1455  Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1456  MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1457  }
1458 
1459  // Set the current number of recycled blocks and subspace
1460  // dimension with the Block GCRO-DR iteration.
1461  if (keff != keff_new) {
1462  keff = keff_new;
1463  block_gcrodr_iter->setSize( keff, numBlocks_ );
1464  // Important to zero this out before next cyle
1465  SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1466  b1.putScalar(zero);
1467  }
1468 
1469 } //end buildRecycleSpaceAugKryl definition
1470 
1471 template<class ScalarType, class MV, class OP>
1473  int i, j;
1474  int m2 = GG.numCols();
1475  bool xtraVec = false;
1476  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1477  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1478  std::vector<int> index;
1479 
1480  // Real and imaginary eigenvalue components
1481  std::vector<MagnitudeType> wr(m2), wi(m2);
1482 
1483  // Magnitude of harmonic Ritz values
1484  std::vector<MagnitudeType> w(m2);
1485 
1486  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1487  SDM vr(m2,m2,false);
1488 
1489  // Sorted order of harmonic Ritz values
1490  std::vector<int> iperm(m2);
1491 
1492  // Set flag indicating recycle space has been generated this solve
1493  builtRecycleSpace_ = true;
1494 
1495  // Form matrices for generalized eigenproblem
1496 
1497  // B = G' * G; Don't zero out matrix when constructing
1498  SDM B(m2,m2,false);
1499  B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
1500 
1501  // A_tmp = | C'*U 0 |
1502  // | V_{m+1}'*U I |
1503  SDM A_tmp( keff+m+blockSize_, keff+m );
1504 
1505 
1506  // A_tmp(1:keff,1:keff) = C' * U;
1507  index.resize(keff);
1508  for (int i=0; i<keff; ++i) { index[i] = i; }
1509  Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1510  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1511  SDM A11( Teuchos::View, A_tmp, keff, keff );
1512  MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1513 
1514  // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
1515  SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1516  index.resize(m+blockSize_);
1517  for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1518  Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1519  MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1520 
1521  // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
1522  for( i=keff; i<keff+m; i++)
1523  A_tmp(i,i) = one;
1524 
1525  // A = G' * A_tmp;
1526  SDM A( m2, A_tmp.numCols() );
1527  A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1528 
1529  // Compute k smallest harmonic Ritz pairs
1530  // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
1531  // ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
1532  // IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
1533  // RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
1534  // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
1535  char balanc='P', jobvl='N', jobvr='V', sense='N';
1536  int ld = A.numRows();
1537  int lwork = 6*ld;
1538  int ldvl = ld, ldvr = ld;
1539  int info = 0,ilo = 0,ihi = 0;
1540  MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1541  ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
1542  std::vector<ScalarType> beta(ld);
1543  std::vector<ScalarType> work(lwork);
1544  std::vector<MagnitudeType> rwork(lwork);
1545  std::vector<MagnitudeType> lscale(ld), rscale(ld);
1546  std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1547  std::vector<int> iwork(ld+6);
1548  int *bwork = 0; // If sense == 'N', bwork is never referenced
1549  //lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1550  // &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1551  // &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
1552  lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1553  &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1554  &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1555  &iwork[0], bwork, &info);
1556  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1557 
1558  // Construct magnitude of each harmonic Ritz value
1559  // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
1560  for( i=0; i<ld; i++ ) // Construct magnitude of each harmonic Ritz value
1561  w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1562 
1563  this->sort(w,ld,iperm);
1564 
1565  bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1566 
1567  // Select recycledBlocks_ smallest eigenvectors
1568  for( i=0; i<recycledBlocks_; i++ )
1569  for( j=0; j<ld; j++ )
1570  PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1571 
1572  if(scalarTypeIsComplex==false) {
1573 
1574  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1575  if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1576  int countImag = 0;
1577  for ( i=ld-recycledBlocks_; i<ld; i++ )
1578  if (wi[iperm[i]] != 0.0) countImag++;
1579  // Check to see if this count is even or odd:
1580  if (countImag % 2) xtraVec = true;
1581  }
1582 
1583  if (xtraVec) { // we need to store one more vector
1584  if (wi[iperm[ld-recycledBlocks_]] > 0.0) { // I picked the "real" component
1585  for( j=0; j<ld; j++ ) // so get the "imag" component
1586  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1587  }
1588  else { // I picked the "imag" component
1589  for( j=0; j<ld; j++ ) // so get the "real" component
1590  PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1591  }
1592  }
1593 
1594  }
1595 
1596  // Return whether we needed to store an additional vector
1597  if (xtraVec)
1598  return recycledBlocks_+1;
1599  else
1600  return recycledBlocks_;
1601 
1602 } //end getHarmonicVecsAugKryl definition
1603 
1604 template<class ScalarType, class MV, class OP>
1606  bool xtraVec = false;
1607  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1608  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1609 
1610  // Real and imaginary eigenvalue components
1611  std::vector<MagnitudeType> wr(m), wi(m);
1612 
1613  // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1614  SDM vr(m,m,false);
1615 
1616  // Magnitude of harmonic Ritz values
1617  std::vector<MagnitudeType> w(m);
1618 
1619  // Sorted order of harmonic Ritz values, also used for DGEEV
1620  std::vector<int> iperm(m);
1621 
1622  // Output info
1623  int info = 0;
1624 
1625  // Set flag indicating recycle space has been generated this solve
1626  builtRecycleSpace_ = true;
1627 
1628  // Solve linear system: H_m^{-H}*E_m where E_m is the last blockSize_ columns of the identity matrix
1629  SDM HHt( HH, Teuchos::TRANS );
1630  Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp( new SDM( m, blockSize_));
1631 
1632  //Initialize harmRitzMatrix as E_m
1633  for(int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1634 
1635  //compute harmRitzMatrix <- H_m^{-H}*E_m
1636  lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1637 
1638  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1639  // Compute H_m + H_m^{-H}*E_m*H_lbl^{H}*H_lbl
1640  // H_lbl is bottom-right block of H_, which is a blockSize_ x blockSize_ matrix
1641 
1642  Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1644 
1645  { //So that HTemp will fall out of scope
1646 
1647  // HH_lbl_t <- H_lbl_t*H_lbl
1649  Htemp = Teuchos::rcp(new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1650  Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1651  H_lbl_t.assign(*Htemp);
1652  // harmRitzMatrix <- harmRitzMatrix*HH_lbl_t
1653  Htemp = Teuchos::rcp(new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1654  Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1655  harmRitzMatrix -> assign(*Htemp);
1656 
1657  // We need to add harmRitzMatrix to the last blockSize_ columns of HH and store the result harmRitzMatrix
1658  int harmColIndex, HHColIndex;
1659  Htemp = Teuchos::rcp(new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1660  for(int i = 0; i<blockSize_; i++){
1661  harmColIndex = harmRitzMatrix -> numCols() - i -1;
1662  HHColIndex = m-i-1;
1663  for(int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1664  }
1665  harmRitzMatrix = Htemp;
1666  }
1667 
1668  // Revise to do query for optimal workspace first
1669  // Create simple storage for the left eigenvectors, which we don't care about.
1670 
1671  // Size of workspace and workspace for DGEEV
1672  int lwork = -1;
1673  std::vector<ScalarType> work(1);
1674  std::vector<MagnitudeType> rwork(2*m);
1675 
1676  const int ldvl = 1;
1677  ScalarType* vl = 0;
1678 
1679  // First query GEEV for the optimal workspace size
1680  lapack.GEEV('N', 'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1681  vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1682 
1683  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
1684  work.resize( lwork );
1685 
1686  lapack.GEEV('N', 'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1687  vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1688 
1689  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1690 
1691  // Construct magnitude of each harmonic Ritz value
1692  for( int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1693 
1694  this->sort(w, m, iperm);
1695 
1696  bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1697 
1698  // Select recycledBlocks_ smallest eigenvectors
1699  for( int i=0; i<recycledBlocks_; ++i )
1700  for(int j=0; j<m; j++ )
1701  PP(j,i) = vr(j,iperm[i]);
1702 
1703  if(scalarTypeIsComplex==false) {
1704 
1705  // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1706  if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1707  int countImag = 0;
1708  for (int i=0; i<recycledBlocks_; ++i )
1709  if (wi[iperm[i]] != 0.0) countImag++;
1710  // Check to see if this count is even or odd:
1711  if (countImag % 2) xtraVec = true;
1712  }
1713 
1714  if (xtraVec) { // we need to store one more vector
1715  if (wi[iperm[recycledBlocks_-1]] > 0.0) { // I picked the "real" component
1716  for(int j=0; j<m; ++j ) // so get the "imag" component
1717  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1718  }
1719  else{ // I picked the "imag" component
1720  for(int j=0; j<m; ++j ) // so get the "real" component
1721  PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1722  }
1723  }
1724 
1725  }
1726 
1727  // Return whether we needed to store an additional vector
1728  if (xtraVec) {
1729  printer_->stream(Debug) << "Recycled " << recycledBlocks_+1 << " vectors" << std::endl;
1730  return recycledBlocks_+1;
1731  }
1732  else {
1733  printer_->stream(Debug) << "Recycled " << recycledBlocks_ << " vectors" << std::endl;
1734  return recycledBlocks_;
1735  }
1736 } //end getHarmonicVecsKryl
1737 
1738 // This method sorts list of n floating-point numbers and return permutation vector
1739 template<class ScalarType, class MV, class OP>
1740 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm) {
1741  int l, r, j, i, flag;
1742  int RR2;
1743  MagnitudeType dRR, dK;
1744 
1745  // Initialize the permutation vector.
1746  for(j=0;j<n;j++)
1747  iperm[j] = j;
1748 
1749  if (n <= 1) return;
1750 
1751  l = n / 2 + 1;
1752  r = n - 1;
1753  l = l - 1;
1754  dRR = dlist[l - 1];
1755  dK = dlist[l - 1];
1756 
1757  RR2 = iperm[l - 1];
1758  while (r != 0) {
1759  j = l;
1760  flag = 1;
1761  while (flag == 1) {
1762  i = j;
1763  j = j + j;
1764  if (j > r + 1)
1765  flag = 0;
1766  else {
1767  if (j < r + 1)
1768  if (dlist[j] > dlist[j - 1]) j = j + 1;
1769  if (dlist[j - 1] > dK) {
1770  dlist[i - 1] = dlist[j - 1];
1771  iperm[i - 1] = iperm[j - 1];
1772  }
1773  else {
1774  flag = 0;
1775  }
1776  }
1777  }
1778  dlist[i - 1] = dRR;
1779  iperm[i - 1] = RR2;
1780  if (l == 1) {
1781  dRR = dlist [r];
1782  RR2 = iperm[r];
1783  dK = dlist[r];
1784  dlist[r] = dlist[0];
1785  iperm[r] = iperm[0];
1786  r = r - 1;
1787  }
1788  else {
1789  l = l - 1;
1790  dRR = dlist[l - 1];
1791  RR2 = iperm[l - 1];
1792  dK = dlist[l - 1];
1793  }
1794  }
1795  dlist[0] = dRR;
1796  iperm[0] = RR2;
1797 } //end sort() definition
1798 
1799 template<class ScalarType, class MV, class OP>
1801  using Teuchos::RCP;
1802  using Teuchos::rcp;
1803  using Teuchos::rcp_const_cast;
1804 
1805  // MLP: NEED TO ADD CHECK IF PARAMETERS ARE SET LATER
1806 
1807  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1808  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1809  std::vector<int> index(numBlocks_+1);
1810 
1811  // MLP: EXCEPTION TESTING SHOULD GO HERE
1812 
1813  //THE FOLLOWING BLOCK OF CODE IS TO INFORM THE PROBLEM HOW MANY RIGHT HAND
1814  //SIDES WILL BE SOLVED. IF THE CHOSEN BLOCK SIZE IS THE SAME AS THE NUMBER
1815  //OF RIGHT HAND SIDES IN THE PROBLEM THEN WE PASS THE PROBLEM
1816  //INDICES FOR ALL RIGHT HAND SIDES. IF BLOCK SIZE IS GREATER THAN
1817  //THE NUMBER OF RIGHT HAND SIDES AND THE USER HAS ADAPTIVE BLOCK SIZING
1818  //TURNED OFF, THEN THE PROBLEM WILL GENERATE RANDOM RIGHT HAND SIDES SO WE HAVE AS MANY
1819  //RIGHT HAND SIDES AS BLOCK SIZE INDICATES. THIS MAY NEED TO BE CHANGED
1820  //LATER TO ACCOMADATE SMARTER CHOOSING OF FICTITIOUS RIGHT HAND SIDES.
1821  //THIS CODE IS DIRECTLY LIFTED FROM BelosBlockGmresSolMgr.hpp.
1822 
1823  // Create indices for the linear systems to be solved.
1824  int startPtr = 0;
1825  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1826  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1827 
1828  //currIdx holds indices to all right-hand sides to be solved
1829  //or -1's to indicate that random right-hand sides should be generated
1830  std::vector<int> currIdx;
1831 
1832  if ( adaptiveBlockSize_ ) {
1833  blockSize_ = numCurrRHS;
1834  currIdx.resize( numCurrRHS );
1835  for (int i=0; i<numCurrRHS; ++i)
1836  currIdx[i] = startPtr+i;
1837  }
1838  else {
1839  currIdx.resize( blockSize_ );
1840  for (int i=0; i<numCurrRHS; ++i)
1841  currIdx[i] = startPtr+i;
1842  for (int i=numCurrRHS; i<blockSize_; ++i)
1843  currIdx[i] = -1;
1844  }
1845 
1846  // Inform the linear problem of the linear systems to solve/generate.
1847  problem_->setLSIndex( currIdx );
1848 
1849  //ADD ERROR CHECKING TO MAKE SURE SIZE OF BLOCK KRYLOV SUBSPACE NOT LARGER THAN dim
1850  //ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1851 
1852  // reset loss of orthogonality flag
1853  loaDetected_ = false;
1854 
1855  // Assume convergence is achieved, then let any failed convergence set this to false.
1856  bool isConverged = true;
1857 
1858  // Initialize storage for all state variables
1859  initializeStateStorage();
1860 
1861  // Parameter list MLP: WE WILL NEED TO ASSIGN SOME PARAMETER VALUES LATER
1862  Teuchos::ParameterList plist;
1863 
1864  while (numRHS2Solve > 0){ //This loops through each block of RHS's being solved
1865  // **************************************************************************************************************************
1866  // Begin initial solver preparations. Either update U,C for new operator or generate an initial space using a cycle of GMRES
1867  // **************************************************************************************************************************
1868  //int prime_iterations; // silence warning about not being used
1869  if(keff > 0){//If there is already a subspace to recycle, then use it
1870  // Update U, C for the new operator
1871 
1872 // TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,BlockGCRODRSolMgrRecyclingFailure,
1873 // "Belos::BlockGCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1874  printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
1875 
1876  // Compute image of U_ under the new operator
1877  index.resize(keff);
1878  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1879  RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1880  RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1881  problem_->apply( *Utmp, *Ctmp );
1882 
1883  RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1884 
1885  // Orthogonalize this block
1886  // Get a matrix to hold the orthonormalization coefficients.
1887  SDM Ftmp( Teuchos::View, *F_, keff, keff );
1888  int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,false));
1889  // Throw an error if we could not orthogonalize this block
1890  TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,BlockGCRODRSolMgrOrthoFailure,"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1891 
1892  // U_ = U_*F^{-1}
1893  // First, compute LU factorization of R
1894  int info = 0;
1895  ipiv_.resize(Ftmp.numRows());
1896  lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1897  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1898 
1899  // Now, form inv(F)
1900  int lwork = Ftmp.numRows();
1901  work_.resize(lwork);
1902  lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1903  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1904 
1905  // U_ = U1_; (via a swap)
1906  MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1907  std::swap(U_, U1_);
1908 
1909  // Must reinitialize after swap
1910  index.resize(keff);
1911  for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1912  Ctmp = MVT::CloneViewNonConst( *C_, index );
1913  Utmp = MVT::CloneView( *U_, index );
1914 
1915  // Compute C_'*R_
1916  SDM Ctr(keff,blockSize_);
1917  problem_->computeCurrPrecResVec( &*R_ );
1918  MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1919 
1920  // Update solution ( x += U_*C_'*R_ )
1921  RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1922  MVT::MvInit( *update, 0.0 );
1923  MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1924  problem_->updateSolution( update, true );
1925 
1926  // Update residual norm ( r -= C_*C_'*R_ )
1927  MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1928 
1929  // We recycled space from previous call
1930  //prime_iterations = 0;
1931 
1932  // Since we have a previous recycle space, we do not need block_gmres_iter
1933  // and we must allocate V ourselves. See comments in initialize() routine for
1934  // further explanation.
1935  if (V_ == Teuchos::null) {
1936  // Check if there is a multivector to clone from.
1937  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1938  V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1939  }
1940  else{
1941  // Generate V_ by cloning itself ONLY if more space is needed.
1942  if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1943  Teuchos::RCP<const MV> tmp = V_;
1944  V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1945  }
1946  }
1947  } //end if(keff > 0)
1948  else{ //if there was no subspace to recycle, then build one
1949  printer_->stream(Debug) << " No recycled subspace available for RHS index " << std::endl << std::endl;
1950 
1951  Teuchos::ParameterList primeList;
1952  //we set this as numBlocks-1 because that is the Krylov dimension we want
1953  //so that the size of the Hessenberg created matches that of what we were expecting.
1954  primeList.set("Num Blocks",numBlocks_-1);
1955  primeList.set("Block Size",blockSize_);
1956  primeList.set("Recycled Blocks",0);
1957  primeList.set("Keep Hessenberg",true);
1958  primeList.set("Initialize Hessenberg",true);
1959 
1960  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1961  if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {//if user has selected a total subspace dimension larger than system dimension
1962  ptrdiff_t tmpNumBlocks = 0;
1963  if (blockSize_ == 1)
1964  tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
1965  else{
1966  tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
1967  printer_->stream(Warnings) << "Belos::BlockGCRODRSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
1968  << std::endl << "The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1969  primeList.set("Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1970  }
1971  }
1972  else{
1973  primeList.set("Num Blocks",numBlocks_-1);
1974  }
1975  //Create Block GMRES iteration object to perform one cycle of GMRES
1977  block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
1978 
1979  // MLP: ADD LOGIC TO DEAL WITH USER ASKING TO GENERATE A LARGER SPACE THAN dim AS IN HEIDI'S BlockGmresSolMgr CODE (DIDN'T WE ALREADY DO THIS SOMEWHERE?)
1980  block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1981 
1982  //Create the first block in the current BLOCK Krylov basis (using residual)
1983  Teuchos::RCP<MV> V_0;
1984  if (currIdx[blockSize_-1] == -1) {//if augmented with random RHS's
1985  V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1986  problem_->computeCurrPrecResVec( &*V_0 );
1987  }
1988  else {
1989  V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1990  }
1991 
1992  // Get a matrix to hold the orthonormalization coefficients.
1993  Teuchos::RCP<SDM > z_0 = Teuchos::rcp( new SDM( blockSize_, blockSize_ ) );
1994 
1995  // Orthonormalize the new V_0
1996  int rank = ortho_->normalize( *V_0, z_0 );
1997  (void) rank; // silence compiler warning for unused variable
1998  // MLP: ADD EXCEPTION IF INITIAL BLOCK IS RANK DEFFICIENT
1999 
2000  // Set the new state and initialize the iteration.
2002  newstate.V = V_0;
2003  newstate.z = z_0;
2004  newstate.curDim = 0;
2005  block_gmres_iter->initializeGmres(newstate);
2006 
2007  bool primeConverged = false;
2008 
2009  try {
2010  printer_->stream(Debug) << " Preparing to Iterate!!!!" << std::endl << std::endl;
2011  block_gmres_iter->iterate();
2012 
2013  // *********************
2014  // check for convergence
2015  // *********************
2016  if ( convTest_->getStatus() == Passed ) {
2017  printer_->stream(Debug) << "We converged during the prime the pump stage" << std::endl << std::endl;
2018  primeConverged = !(expConvTest_->getLOADetected());
2019  if ( expConvTest_->getLOADetected() ) {
2020  // we don't have convergence
2021  loaDetected_ = true;
2022  printer_->stream(Warnings) << "Belos::BlockGCRODRSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2023  }
2024  }
2025  // *******************************************
2026  // check for maximum iterations of block GMRES
2027  // *******************************************
2028  else if( maxIterTest_->getStatus() == Passed ) {
2029  // we don't have convergence
2030  primeConverged = false;
2031  }
2032  // ****************************************************************
2033  // We need to recycle and continue, print a message indicating this
2034  // ****************************************************************
2035  else{ //debug statements
2036  printer_->stream(Debug) << " We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2037  }
2038  } //end try
2039  catch (const GmresIterationOrthoFailure &e) {
2040  // If the block size is not one, it's not considered a lucky breakdown.
2041  if (blockSize_ != 1) {
2042  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2043  << block_gmres_iter->getNumIters() << std::endl
2044  << e.what() << std::endl;
2045  if (convTest_->getStatus() != Passed)
2046  primeConverged = false;
2047  }
2048  else {
2049  // If the block size is one, try to recover the most recent least-squares solution
2050  block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2051  // Check to see if the most recent least-squares solution yielded convergence.
2052  sTest_->checkStatus( &*block_gmres_iter );
2053  if (convTest_->getStatus() != Passed)
2054  isConverged = false;
2055  }
2056  } // end catch (const GmresIterationOrthoFailure &e)
2057  catch (const StatusTestNaNError& e) {
2058  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
2059  achievedTol_ = MT::one();
2060  Teuchos::RCP<MV> X = problem_->getLHS();
2061  MVT::MvInit( *X, SCT::zero() );
2062  printer_->stream(Warnings) << "Belos::BlockGCRODRSolMgr::solve(): Warning! NaN has been detected!"
2063  << std::endl;
2064  return Unconverged;
2065  }
2066  catch (const std::exception &e) {
2067  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2068  << block_gmres_iter->getNumIters() << std::endl
2069  << e.what() << std::endl;
2070  throw;
2071  }
2072 
2073  // Record number of iterations in generating initial recycle spacec
2074  //prime_iterations = block_gmres_iter->getNumIters();//instantiated here because it is not needed outside of else{} scope; we'll see if this is true or not
2075 
2076  // Update the linear problem.
2077  RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2078  problem_->updateSolution( update, true );
2079 
2080  // Update the block residual
2081 
2082  problem_->computeCurrPrecResVec( &*R_ );
2083 
2084  // Get the state.
2085  newstate = block_gmres_iter->getState();
2086  int p = newstate.curDim;
2087  (void) p; // silence compiler warning for unused variable
2088  H_->assign(*(newstate.H));//IS THIS A GOOD IDEA? I DID IT SO MEMBERS WOULD HAVE ACCESS TO H.
2089  // Get new Krylov vectors, store in V_
2090 
2091  // Since the block_gmres_iter returns the state
2092  // with const RCP's we need to cast newstate.V as
2093  // a non const RCP. This is okay because block_gmres_iter
2094  // is about to fall out of scope, so we are simply
2095  // rescuing *newstate.V from being destroyed so we can
2096  // use it for future block recycled GMRES cycles
2097  V_ = rcp_const_cast<MV>(newstate.V);
2098  newstate.V.release();
2099  //COMPUTE NEW RECYCLE SPACE SOMEHOW
2100  buildRecycleSpaceKryl(keff, block_gmres_iter);
2101  printer_->stream(Debug) << "Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2102 
2103  // Return to outer loop if the priming solve
2104  // converged, set the next linear system.
2105  if (primeConverged) {
2106  /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2107  // Inform the linear problem that we are
2108  // finished with this block linear system.
2109  problem_->setCurrLS();
2110 
2111  // Update indices for the linear systems to be solved.
2112  // KMS: Fix the numRHS2Solve; Copy from Heidi's BlockGmres code
2113  numRHS2Solve -= 1;
2114  printer_->stream(Debug) << numRHS2Solve << " Right hand sides left to solve" << std::endl;
2115  if ( numRHS2Solve > 0 ) {
2116  currIdx[0]++;
2117  // Set the next indices
2118  problem_->setLSIndex( currIdx );
2119  }
2120  else {
2121  currIdx.resize( numRHS2Solve );
2122  }
2123  *****************************************************************************/
2124 
2125  // Inform the linear problem that we are finished with this block linear system.
2126  problem_->setCurrLS();
2127 
2128  // Update indices for the linear systems to be solved.
2129  startPtr += numCurrRHS;
2130  numRHS2Solve -= numCurrRHS;
2131  if ( numRHS2Solve > 0 ) {
2132  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2133  if ( adaptiveBlockSize_ ) {
2134  blockSize_ = numCurrRHS;
2135  currIdx.resize( numCurrRHS );
2136  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2137  }
2138  else {
2139  currIdx.resize( blockSize_ );
2140  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2141  for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2142  }
2143  // Set the next indices.
2144  problem_->setLSIndex( currIdx );
2145  }
2146  else {
2147  currIdx.resize( numRHS2Solve );
2148  }
2149  continue;//Begin solving the next block of RHS's
2150  } //end if (primeConverged)
2151  } //end else [if(keff > 0)]
2152 
2153  // *****************************************
2154  // Initial subspace update/construction done
2155  // Start cycles of recycled GMRES
2156  // *****************************************
2157  Teuchos::ParameterList blockgcrodrList;
2158  blockgcrodrList.set("Num Blocks",numBlocks_);
2159  blockgcrodrList.set("Block Size",blockSize_);
2160  blockgcrodrList.set("Recycled Blocks",keff);
2161 
2163  block_gcrodr_iter = Teuchos::rcp( new BlockGCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,blockgcrodrList) );
2165 
2166  index.resize( blockSize_ );
2167  for(int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2168  Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2169 
2170  // MLP: MVT::SetBlock(*R_,index,*V0);
2171  MVT::Assign(*R_,*V0);
2172 
2173  index.resize(keff);//resize to get appropriate recycle space vectors
2174  for(int i=0; i < keff; i++){ index[i] = i;};
2175  B_ = rcp(new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2176  H_ = rcp(new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2177 
2178  newstate.V = V_;
2179  newstate.B= B_;
2180  newstate.U = MVT::CloneViewNonConst(*U_, index);
2181  newstate.C = MVT::CloneViewNonConst(*C_, index);
2182  newstate.H = H_;
2183  newstate.curDim = blockSize_;
2184  block_gcrodr_iter -> initialize(newstate);
2185 
2186  int numRestarts = 0;
2187 
2188  while(1){//Each execution of this loop is a cycle of block GCRODR
2189  try{
2190  block_gcrodr_iter -> iterate();
2191 
2192  // ***********************
2193  // Check Convergence First
2194  // ***********************
2195  if( convTest_->getStatus() == Passed ) {
2196  // we have convergence
2197  break;//from while(1)
2198  } //end if converged
2199 
2200  // ***********************************
2201  // Check if maximum iterations reached
2202  // ***********************************
2203  else if(maxIterTest_->getStatus() == Passed ){
2204  // no convergence; hit maxit
2205  isConverged = false;
2206  break; // from while(1)
2207  } //end elseif reached maxit
2208 
2209  // **********************************************
2210  // Check if subspace full; do we need to restart?
2211  // **********************************************
2212  else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2213 
2214  // Update recycled space even if we have reached max number of restarts
2215 
2216  // Update linear problem
2217  Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2218  problem_->updateSolution(update, true);
2219  buildRecycleSpaceAugKryl(block_gcrodr_iter);
2220 
2221  printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2222  // NOTE: If we have hit the maximum number of restarts, then we will quit
2223  if(numRestarts >= maxRestarts_) {
2224  isConverged = false;
2225  break; //from while(1)
2226  } //end if max restarts
2227 
2228  numRestarts++;
2229 
2230  printer_ -> stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
2231 
2232  // Create the restart vector (first block in the current Krylov basis)
2233  problem_->computeCurrPrecResVec( &*R_ );
2234  index.resize( blockSize_ );
2235  for (int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2236  Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2237  MVT::SetBlock(*R_,index,*V0);
2238 
2239  // Set the new state and initialize the solver.
2241  index.resize( numBlocks_*blockSize_ );
2242  for (int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2243  restartState.V = MVT::CloneViewNonConst( *V_, index );
2244  index.resize( keff );
2245  for (int ii=0; ii<keff; ++ii) index[ii] = ii;
2246  restartState.U = MVT::CloneViewNonConst( *U_, index );
2247  restartState.C = MVT::CloneViewNonConst( *C_, index );
2248  B_ = rcp(new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2249  H_ = rcp(new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2250  restartState.B = B_;
2251  restartState.H = H_;
2252  restartState.curDim = blockSize_;
2253  block_gcrodr_iter->initialize(restartState);
2254 
2255  } //end else if need to restart
2256 
2257  // ****************************************************************
2258  // We returned from iterate(), but none of our status tests passed.
2259  // Something is wrong, and it is probably our fault.
2260  // ****************************************************************
2261  else {
2262  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2263  } //end else (no status test passed)
2264 
2265  }//end try
2266  catch(const BlockGCRODRIterOrthoFailure &e){ //there was an exception
2267  // Try to recover the most recent least-squares solution
2268  block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2269  // Check to see if the most recent least-squares solution yielded convergence.
2270  sTest_->checkStatus( &*block_gcrodr_iter );
2271  if (convTest_->getStatus() != Passed) isConverged = false;
2272  break;
2273  } // end catch orthogonalization failure
2274  catch(const std::exception &e){
2275  printer_->stream(Errors) << "Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
2276  << block_gcrodr_iter->getNumIters() << std::endl
2277  << e.what() << std::endl;
2278  throw;
2279  } //end catch standard exception
2280  } //end while(1)
2281 
2282  // Compute the current solution.
2283  // Update the linear problem.
2284  Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2285  problem_->updateSolution( update, true );
2286 
2287  /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2288  // Inform the linear problem that we are finished with this block linear system.
2289  problem_->setCurrLS();
2290 
2291  // Update indices for the linear systems to be solved.
2292  numRHS2Solve -= 1;
2293  if ( numRHS2Solve > 0 ) {
2294  currIdx[0]++;
2295  // Set the next indices.
2296  problem_->setLSIndex( currIdx );
2297  }
2298  else {
2299  currIdx.resize( numRHS2Solve );
2300  }
2301  ******************************************************************************/
2302 
2303  // Inform the linear problem that we are finished with this block linear system.
2304  problem_->setCurrLS();
2305 
2306  // Update indices for the linear systems to be solved.
2307  startPtr += numCurrRHS;
2308  numRHS2Solve -= numCurrRHS;
2309  if ( numRHS2Solve > 0 ) {
2310  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2311  if ( adaptiveBlockSize_ ) {
2312  blockSize_ = numCurrRHS;
2313  currIdx.resize( numCurrRHS );
2314  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2315  }
2316  else {
2317  currIdx.resize( blockSize_ );
2318  for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2319  for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2320  }
2321  // Set the next indices.
2322  problem_->setLSIndex( currIdx );
2323  }
2324  else {
2325  currIdx.resize( numRHS2Solve );
2326  }
2327 
2328  // If we didn't build a recycle space this solve but ran at least k iterations, force build of new recycle space
2329  if (!builtRecycleSpace_) {
2330  buildRecycleSpaceAugKryl(block_gcrodr_iter);
2331  printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2332  }
2333  }//end while (numRHS2Solve > 0)
2334 
2335  // print final summary
2336  sTest_->print( printer_->stream(FinalSummary) );
2337 
2338  // print timing information
2339  #ifdef BELOS_TEUCHOS_TIME_MONITOR
2340  // Calling summarize() can be expensive, so don't call unless the user wants to print out timing details.
2341  // summarize() will do all the work even if it's passed a "black hole" output stream.
2342  if (verbosity_ & TimingDetails) Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
2343  #endif
2344  // get iteration information for this solve
2345  numIters_ = maxIterTest_->getNumIters();
2346 
2347  // get residual information for this solve
2348  const std::vector<MagnitudeType>* pTestValues = NULL;
2349  pTestValues = impConvTest_->getTestValue();
2350  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2351  "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2352  "getTestValue() method returned NULL. Please report this bug to the "
2353  "Belos developers.");
2354  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2355  "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2356  "getTestValue() method returned a vector of length zero. Please report "
2357  "this bug to the Belos developers.");
2358  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
2359  // achieved tolerances for all vectors in the current solve(), or
2360  // just for the vectors from the last deflation?
2361  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
2362 
2363  if (!isConverged) return Unconverged; // return from BlockGCRODRSolMgr::solve()
2364  return Converged; // return from BlockGCRODRSolMgr::solve()
2365 } //end solve()
2366 
2367 } //End Belos Namespace
2368 
2369 #endif /* BELOS_BLOCK_GCRODR_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
ScalarType * values() const
Collection of types and exceptions used within the Belos solvers.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
OperatorTraits< ScalarType, MV, OP > OPT
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
int getHarmonicVecsAugKryl(int keff, int m, const SDM &HH, const Teuchos::RCP< const MV > &VV, SDM &PP)
Teuchos::LAPACK< int, ScalarType > lapack
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
static const bool adaptiveBlockSize_default_
virtual ~BlockGCRODRSolMgr()
Destructor.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
static void clearCounter(const std::string &name)
bool loaDetected_
Whether a loss of accuracy was detected during the solve.
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
void buildRecycleSpaceAugKryl(Teuchos::RCP< BlockGCRODRIter< ScalarType, MV, OP > > gcrodr_iter)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
T & get(ParameterList &l, const std::string &name)
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
Teuchos::ScalarTraits< MagnitudeType > SMT
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Thrown when the linear problem was not set up correctly.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
static magnitudeType real(T a)
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set the parameters the solver should use to solve the linear problem.
bool isSet_
Whether setParameters() successfully finished setting parameters.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::SerialDenseVector< int, ScalarType > SDV
Structure to contain pointers to GmresIteration state variables.
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::RCP< Teuchos::ParameterList > params_
This solver&#39;s current parameter list.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
Thrown if any problem occurs in using or creating the recycle subspace.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< MV > V
The current Krylov basis.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
int curDim
The current dimension of the reduction.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< Teuchos::Time > timerSolve_
Timer for solve().
Teuchos::ScalarTraits< MagnitudeType > MT
void sort(std::vector< MagnitudeType > &dlist, int n, std::vector< int > &iperm)
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The current linear problem to solve.
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
int curDim
The current dimension of the reduction.
static const std::string recycleMethod_default_
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
bool is_null(const RCP< T > &p)
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
std::vector< ScalarType > work_
BlockGCRODRSolMgr()
Default constructor.
A linear system to solve, and its associated information.
std::string description() const
A description of the Block GCRODR solver manager.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< OutputManager< ScalarType > > printer_
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Thrown when the solution manager was unable to orthogonalize the basis vectors.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > impConvTest_
Belos concrete class for performing the block GMRES iteration.
Ptr< T > release()
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
Teuchos::RCP< std::ostream > outputStream_
void buildRecycleSpaceKryl(int &keff, Teuchos::RCP< BlockGmresIter< ScalarType, MV, OP > > block_gmres_iter)
static magnitudeType magnitude(T a)
OrdinalType numCols() const
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
int getHarmonicVecsKryl(int m, const SDM &HH, SDM &PP)
Class which defines basic traits for the operator type.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
std::vector< ScalarType > tau_
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos::RCP< const Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
ReturnType solve()
Solve the current linear problem.
bool builtRecycleSpace_
Whether we have generated or regenerated a recycle space yet this solve.
Thrown when an LAPACK call fails.
OrdinalType stride() const
ortho_factory_type orthoFactory_
Factory for creating MatOrthoManager subclass instances.
OrdinalType numRows() const
Belos concrete class for performing the block, flexible GMRES iteration.
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
OutputType
Style of output used to display status test information.
Definition: BelosTypes.hpp:141
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.