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