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