Belos  Version of the Day
BelosRCGSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_RCG_SOLMGR_HPP
43 #define BELOS_RCG_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosRCGIter.hpp"
58 #include "BelosStatusTestCombo.hpp"
60 #include "BelosOutputManager.hpp"
61 #include "Teuchos_LAPACK.hpp"
62 #include "Teuchos_as.hpp"
63 #ifdef BELOS_TEUCHOS_TIME_MONITOR
64 #include "Teuchos_TimeMonitor.hpp"
65 #endif
66 
106 namespace Belos {
107 
109 
110 
118  RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
119  {}};
120 
127  class RCGSolMgrLAPACKFailure : public BelosError {public:
128  RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
129  {}};
130 
137  class RCGSolMgrRecyclingFailure : public BelosError {public:
138  RCGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
139  {}};
140 
142 
143 
144  // Partial specialization for unsupported ScalarType types.
145  // This contains a stub implementation.
146  template<class ScalarType, class MV, class OP,
147  const bool supportsScalarType =
149  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
150  class RCGSolMgr :
151  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
152  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
153  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
154  {
155  static const bool scalarTypeIsSupported =
157  ! Teuchos::ScalarTraits<ScalarType>::isComplex;
158  typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
159  scalarTypeIsSupported> base_type;
160 
161  public:
163  base_type ()
164  {}
165  RCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
166  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
167  base_type ()
168  {}
169  virtual ~RCGSolMgr () {}
170 
172  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
173  return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP,supportsScalarType>);
174  }
175  };
176 
177  // Partial specialization for real ScalarType.
178  // This contains the actual working implementation of RCG.
179  // See discussion in the class documentation above.
180  template<class ScalarType, class MV, class OP>
181  class RCGSolMgr<ScalarType, MV, OP, true> :
182  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
183  private:
186  typedef Teuchos::ScalarTraits<ScalarType> SCT;
187  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
188  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
189 
190  public:
191 
193 
194 
200  RCGSolMgr();
201 
223  RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
224  const Teuchos::RCP<Teuchos::ParameterList> &pl );
225 
227  virtual ~RCGSolMgr() {};
228 
230  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
231  return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP>);
232  }
234 
236 
237 
238  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
239  return *problem_;
240  }
241 
243  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
244 
246  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
247 
253  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
254  return Teuchos::tuple(timerSolve_);
255  }
256 
261  MagnitudeType achievedTol() const override {
262  return achievedTol_;
263  }
264 
266  int getNumIters() const override {
267  return numIters_;
268  }
269 
271  bool isLOADetected() const override { return false; }
272 
274 
276 
277 
279  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
280 
282  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
283 
285 
287 
288 
294  void reset( const ResetType type ) override {
295  if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
296  else if (type & Belos::RecycleSubspace) existU_ = false;
297  }
299 
301 
302 
320  ReturnType solve() override;
321 
323 
326 
328  std::string description() const override;
329 
331 
332  private:
333 
334  // Called by all constructors; Contains init instructions common to all constructors
335  void init();
336 
337  // Computes harmonic eigenpairs of projected matrix created during one cycle.
338  // Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude.
339  void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
340  const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
341  Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
342 
343  // Sort list of n floating-point numbers and return permutation vector
344  void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
345 
346  // Initialize solver state storage
347  void initializeStateStorage();
348 
349  // Linear problem.
350  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
351 
352  // Output manager.
353  Teuchos::RCP<OutputManager<ScalarType> > printer_;
354  Teuchos::RCP<std::ostream> outputStream_;
355 
356  // Status test.
357  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
358  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
359  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
360  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
361 
362  // Current parameter list.
363  Teuchos::RCP<Teuchos::ParameterList> params_;
364 
365  // Default solver values.
366  static constexpr int maxIters_default_ = 1000;
367  static constexpr int blockSize_default_ = 1;
368  static constexpr int numBlocks_default_ = 25;
369  static constexpr int recycleBlocks_default_ = 3;
370  static constexpr bool showMaxResNormOnly_default_ = false;
371  static constexpr int verbosity_default_ = Belos::Errors;
372  static constexpr int outputStyle_default_ = Belos::General;
373  static constexpr int outputFreq_default_ = -1;
374  static constexpr const char * label_default_ = "Belos";
375 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
376 #if defined(_WIN32) && defined(__clang__)
377  static constexpr std::ostream * outputStream_default_ =
378  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
379 #else
380  static constexpr std::ostream * outputStream_default_ = &std::cout;
381 #endif
382 
383  //
384  // Current solver values.
385  //
386 
388  MagnitudeType convtol_;
389 
394  MagnitudeType achievedTol_;
395 
397  int maxIters_;
398 
400  int numIters_;
401 
402  int numBlocks_, recycleBlocks_;
403  bool showMaxResNormOnly_;
404  int verbosity_, outputStyle_, outputFreq_;
405 
407  // Solver State Storage
409  // Search vectors
410  Teuchos::RCP<MV> P_;
411  //
412  // A times current search direction
413  Teuchos::RCP<MV> Ap_;
414  //
415  // Residual vector
416  Teuchos::RCP<MV> r_;
417  //
418  // Preconditioned residual
419  Teuchos::RCP<MV> z_;
420  //
421  // Flag indicating that the recycle space should be used
422  bool existU_;
423  //
424  // Flag indicating that the updated recycle space has been created
425  bool existU1_;
426  //
427  // Recycled subspace and its image
428  Teuchos::RCP<MV> U_, AU_;
429  //
430  // Recycled subspace for next system and its image
431  Teuchos::RCP<MV> U1_;
432  //
433  // Coefficients arising in RCG iteration
434  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
435  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
436  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
437  //
438  // Solutions to local least-squares problems
439  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
440  //
441  // The matrix U^T A U
442  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
443  //
444  // LU factorization of U^T A U
445  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
446  //
447  // Data from LU factorization of UTAU
448  Teuchos::RCP<std::vector<int> > ipiv_;
449  //
450  // The matrix (AU)^T AU
451  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
452  //
453  // The scalar r'*z
454  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
455  //
456  // Matrices needed for calculation of harmonic Ritz eigenproblem
457  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
458  //
459  // Matrices needed for updating recycle space
460  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
461  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
462  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
463  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
464  Teuchos::RCP<MV> U1Y1_, PY2_;
465  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
466  ScalarType dold;
468 
469  // Timers.
470  std::string label_;
471  Teuchos::RCP<Teuchos::Time> timerSolve_;
472 
473  // Internal state variables.
474  bool params_Set_;
475  };
476 
477 
478 // Empty Constructor
479 template<class ScalarType, class MV, class OP>
481  achievedTol_(0.0),
482  numIters_(0)
483 {
484  init();
485 }
486 
487 // Basic Constructor
488 template<class ScalarType, class MV, class OP>
490  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
491  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
492  problem_(problem),
493  achievedTol_(0.0),
494  numIters_(0)
495 {
496  init();
497  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
498 
499  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
500  if ( !is_null(pl) ) {
501  setParameters( pl );
502  }
503 }
504 
505 // Common instructions executed in all constructors
506 template<class ScalarType, class MV, class OP>
508 {
509  outputStream_ = Teuchos::rcp(outputStream_default_,false);
511  maxIters_ = maxIters_default_;
512  numBlocks_ = numBlocks_default_;
513  recycleBlocks_ = recycleBlocks_default_;
514  verbosity_ = verbosity_default_;
515  outputStyle_= outputStyle_default_;
516  outputFreq_= outputFreq_default_;
517  showMaxResNormOnly_ = showMaxResNormOnly_default_;
518  label_ = label_default_;
519  params_Set_ = false;
520  P_ = Teuchos::null;
521  Ap_ = Teuchos::null;
522  r_ = Teuchos::null;
523  z_ = Teuchos::null;
524  existU_ = false;
525  existU1_ = false;
526  U_ = Teuchos::null;
527  AU_ = Teuchos::null;
528  U1_ = Teuchos::null;
529  Alpha_ = Teuchos::null;
530  Beta_ = Teuchos::null;
531  D_ = Teuchos::null;
532  Delta_ = Teuchos::null;
533  UTAU_ = Teuchos::null;
534  LUUTAU_ = Teuchos::null;
535  ipiv_ = Teuchos::null;
536  AUTAU_ = Teuchos::null;
537  rTz_old_ = Teuchos::null;
538  F_ = Teuchos::null;
539  G_ = Teuchos::null;
540  Y_ = Teuchos::null;
541  L2_ = Teuchos::null;
542  DeltaL2_ = Teuchos::null;
543  AU1TUDeltaL2_ = Teuchos::null;
544  AU1TAU1_ = Teuchos::null;
545  AU1TU1_ = Teuchos::null;
546  AU1TAP_ = Teuchos::null;
547  FY_ = Teuchos::null;
548  GY_ = Teuchos::null;
549  APTAP_ = Teuchos::null;
550  U1Y1_ = Teuchos::null;
551  PY2_ = Teuchos::null;
552  AUTAP_ = Teuchos::null;
553  AU1TU_ = Teuchos::null;
554  dold = 0.;
555 }
556 
557 template<class ScalarType, class MV, class OP>
558 void RCGSolMgr<ScalarType,MV,OP,true>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
559 {
560  // Create the internal parameter list if ones doesn't already exist.
561  if (params_ == Teuchos::null) {
562  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
563  }
564  else {
565  params->validateParameters(*getValidParameters());
566  }
567 
568  // Check for maximum number of iterations
569  if (params->isParameter("Maximum Iterations")) {
570  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
571 
572  // Update parameter in our list and in status test.
573  params_->set("Maximum Iterations", maxIters_);
574  if (maxIterTest_!=Teuchos::null)
575  maxIterTest_->setMaxIters( maxIters_ );
576  }
577 
578  // Check for the maximum number of blocks.
579  if (params->isParameter("Num Blocks")) {
580  numBlocks_ = params->get("Num Blocks",numBlocks_default_);
581  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
582  "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
583 
584  // Update parameter in our list.
585  params_->set("Num Blocks", numBlocks_);
586  }
587 
588  // Check for the maximum number of blocks.
589  if (params->isParameter("Num Recycled Blocks")) {
590  recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_);
591  TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
592  "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
593 
594  TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
595  "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
596 
597  // Update parameter in our list.
598  params_->set("Num Recycled Blocks", recycleBlocks_);
599  }
600 
601  // Check to see if the timer label changed.
602  if (params->isParameter("Timer Label")) {
603  std::string tempLabel = params->get("Timer Label", label_default_);
604 
605  // Update parameter in our list and solver timer
606  if (tempLabel != label_) {
607  label_ = tempLabel;
608  params_->set("Timer Label", label_);
609  std::string solveLabel = label_ + ": RCGSolMgr total solve time";
610 #ifdef BELOS_TEUCHOS_TIME_MONITOR
611  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
612 #endif
613  }
614  }
615 
616  // Check for a change in verbosity level
617  if (params->isParameter("Verbosity")) {
618  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
619  verbosity_ = params->get("Verbosity", verbosity_default_);
620  } else {
621  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
622  }
623 
624  // Update parameter in our list.
625  params_->set("Verbosity", verbosity_);
626  if (printer_ != Teuchos::null)
627  printer_->setVerbosity(verbosity_);
628  }
629 
630  // Check for a change in output style
631  if (params->isParameter("Output Style")) {
632  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
633  outputStyle_ = params->get("Output Style", outputStyle_default_);
634  } else {
635  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
636  }
637 
638  // Reconstruct the convergence test if the explicit residual test is not being used.
639  params_->set("Output Style", outputStyle_);
640  outputTest_ = Teuchos::null;
641  }
642 
643  // output stream
644  if (params->isParameter("Output Stream")) {
645  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
646 
647  // Update parameter in our list.
648  params_->set("Output Stream", outputStream_);
649  if (printer_ != Teuchos::null)
650  printer_->setOStream( outputStream_ );
651  }
652 
653  // frequency level
654  if (verbosity_ & Belos::StatusTestDetails) {
655  if (params->isParameter("Output Frequency")) {
656  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
657  }
658 
659  // Update parameter in out list and output status test.
660  params_->set("Output Frequency", outputFreq_);
661  if (outputTest_ != Teuchos::null)
662  outputTest_->setOutputFrequency( outputFreq_ );
663  }
664 
665  // Create output manager if we need to.
666  if (printer_ == Teuchos::null) {
667  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
668  }
669 
670  // Convergence
671  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
672  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
673 
674  // Check for convergence tolerance
675  if (params->isParameter("Convergence Tolerance")) {
676  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
677  convtol_ = params->get ("Convergence Tolerance",
678  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
679  }
680  else {
681  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
682  }
683 
684  // Update parameter in our list and residual tests.
685  params_->set("Convergence Tolerance", convtol_);
686  if (convTest_ != Teuchos::null)
687  convTest_->setTolerance( convtol_ );
688  }
689 
690  if (params->isParameter("Show Maximum Residual Norm Only")) {
691  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
692 
693  // Update parameter in our list and residual tests
694  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
695  if (convTest_ != Teuchos::null)
696  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
697  }
698 
699  // Create status tests if we need to.
700 
701  // Basic test checks maximum iterations and native residual.
702  if (maxIterTest_ == Teuchos::null)
703  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
704 
705  // Implicit residual test, using the native residual to determine if convergence was achieved.
706  if (convTest_ == Teuchos::null)
707  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
708 
709  if (sTest_ == Teuchos::null)
710  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
711 
712  if (outputTest_ == Teuchos::null) {
713 
714  // Create the status test output class.
715  // This class manages and formats the output from the status test.
716  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
717  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
718 
719  // Set the solver string for the output test
720  std::string solverDesc = " Recycling CG ";
721  outputTest_->setSolverDesc( solverDesc );
722  }
723 
724  // Create the timer if we need to.
725  if (timerSolve_ == Teuchos::null) {
726  std::string solveLabel = label_ + ": RCGSolMgr total solve time";
727 #ifdef BELOS_TEUCHOS_TIME_MONITOR
728  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
729 #endif
730  }
731 
732  // Inform the solver manager that the current parameters were set.
733  params_Set_ = true;
734 }
735 
736 
737 template<class ScalarType, class MV, class OP>
738 Teuchos::RCP<const Teuchos::ParameterList>
740 {
741  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
742 
743  // Set all the valid parameters and their default values.
744  if(is_null(validPL)) {
745  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
746  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
747  "The relative residual tolerance that needs to be achieved by the\n"
748  "iterative solver in order for the linear system to be declared converged.");
749  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
750  "The maximum number of block iterations allowed for each\n"
751  "set of RHS solved.");
752  pl->set("Block Size", static_cast<int>(blockSize_default_),
753  "Block Size Parameter -- currently must be 1 for RCG");
754  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
755  "The length of a cycle (and this max number of search vectors kept)\n");
756  pl->set("Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
757  "The number of vectors in the recycle subspace.");
758  pl->set("Verbosity", static_cast<int>(verbosity_default_),
759  "What type(s) of solver information should be outputted\n"
760  "to the output stream.");
761  pl->set("Output Style", static_cast<int>(outputStyle_default_),
762  "What style is used for the solver information outputted\n"
763  "to the output stream.");
764  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
765  "How often convergence information should be outputted\n"
766  "to the output stream.");
767  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
768  "A reference-counted pointer to the output stream where all\n"
769  "solver output is sent.");
770  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
771  "When convergence information is printed, only show the maximum\n"
772  "relative residual norm when the block size is greater than one.");
773  pl->set("Timer Label", static_cast<const char *>(label_default_),
774  "The string to use as a prefix for the timer labels.");
775  validPL = pl;
776  }
777  return validPL;
778 }
779 
780 // initializeStateStorage
781 template<class ScalarType, class MV, class OP>
783 
784  // Check if there is any multivector to clone from.
785  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
786  if (rhsMV == Teuchos::null) {
787  // Nothing to do
788  return;
789  }
790  else {
791 
792  // Initialize the state storage
793  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
794  "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
795 
796  // If the subspace has not been initialized before, generate it using the RHS from lp_.
797  if (P_ == Teuchos::null) {
798  P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
799  }
800  else {
801  // Generate P_ by cloning itself ONLY if more space is needed.
802  if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
803  Teuchos::RCP<const MV> tmp = P_;
804  P_ = MVT::Clone( *tmp, numBlocks_+2 );
805  }
806  }
807 
808  // Generate Ap_ only if it doesn't exist
809  if (Ap_ == Teuchos::null)
810  Ap_ = MVT::Clone( *rhsMV, 1 );
811 
812  // Generate r_ only if it doesn't exist
813  if (r_ == Teuchos::null)
814  r_ = MVT::Clone( *rhsMV, 1 );
815 
816  // Generate z_ only if it doesn't exist
817  if (z_ == Teuchos::null)
818  z_ = MVT::Clone( *rhsMV, 1 );
819 
820  // If the recycle space has not been initialized before, generate it using the RHS from lp_.
821  if (U_ == Teuchos::null) {
822  U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
823  }
824  else {
825  // Generate U_ by cloning itself ONLY if more space is needed.
826  if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
827  Teuchos::RCP<const MV> tmp = U_;
828  U_ = MVT::Clone( *tmp, recycleBlocks_ );
829  }
830  }
831 
832  // If the recycle space has not be initialized before, generate it using the RHS from lp_.
833  if (AU_ == Teuchos::null) {
834  AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
835  }
836  else {
837  // Generate AU_ by cloning itself ONLY if more space is needed.
838  if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
839  Teuchos::RCP<const MV> tmp = AU_;
840  AU_ = MVT::Clone( *tmp, recycleBlocks_ );
841  }
842  }
843 
844  // If the recycle space has not been initialized before, generate it using the RHS from lp_.
845  if (U1_ == Teuchos::null) {
846  U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
847  }
848  else {
849  // Generate U1_ by cloning itself ONLY if more space is needed.
850  if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
851  Teuchos::RCP<const MV> tmp = U1_;
852  U1_ = MVT::Clone( *tmp, recycleBlocks_ );
853  }
854  }
855 
856  // Generate Alpha_ only if it doesn't exist, otherwise resize it.
857  if (Alpha_ == Teuchos::null)
858  Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
859  else {
860  if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
861  Alpha_->reshape( numBlocks_, 1 );
862  }
863 
864  // Generate Beta_ only if it doesn't exist, otherwise resize it.
865  if (Beta_ == Teuchos::null)
866  Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
867  else {
868  if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
869  Beta_->reshape( numBlocks_ + 1, 1 );
870  }
871 
872  // Generate D_ only if it doesn't exist, otherwise resize it.
873  if (D_ == Teuchos::null)
874  D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
875  else {
876  if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
877  D_->reshape( numBlocks_, 1 );
878  }
879 
880  // Generate Delta_ only if it doesn't exist, otherwise resize it.
881  if (Delta_ == Teuchos::null)
882  Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
883  else {
884  if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
885  Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
886  }
887 
888  // Generate UTAU_ only if it doesn't exist, otherwise resize it.
889  if (UTAU_ == Teuchos::null)
890  UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
891  else {
892  if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
893  UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
894  }
895 
896  // Generate LUUTAU_ only if it doesn't exist, otherwise resize it.
897  if (LUUTAU_ == Teuchos::null)
898  LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
899  else {
900  if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
901  LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
902  }
903 
904  // Generate ipiv_ only if it doesn't exist, otherwise resize it.
905  if (ipiv_ == Teuchos::null)
906  ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) );
907  else {
908  if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it
909  ipiv_->resize(recycleBlocks_);
910  }
911 
912  // Generate AUTAU_ only if it doesn't exist, otherwise resize it.
913  if (AUTAU_ == Teuchos::null)
914  AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
915  else {
916  if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
917  AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
918  }
919 
920  // Generate rTz_old_ only if it doesn't exist
921  if (rTz_old_ == Teuchos::null)
922  rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
923  else {
924  if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
925  rTz_old_->reshape( 1, 1 );
926  }
927 
928  // Generate F_ only if it doesn't exist
929  if (F_ == Teuchos::null)
930  F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
931  else {
932  if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
933  F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
934  }
935 
936  // Generate G_ only if it doesn't exist
937  if (G_ == Teuchos::null)
938  G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
939  else {
940  if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
941  G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
942  }
943 
944  // Generate Y_ only if it doesn't exist
945  if (Y_ == Teuchos::null)
946  Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
947  else {
948  if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
949  Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
950  }
951 
952  // Generate L2_ only if it doesn't exist
953  if (L2_ == Teuchos::null)
954  L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
955  else {
956  if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
957  L2_->reshape( numBlocks_+1, numBlocks_ );
958  }
959 
960  // Generate DeltaL2_ only if it doesn't exist
961  if (DeltaL2_ == Teuchos::null)
962  DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
963  else {
964  if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
965  DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
966  }
967 
968  // Generate AU1TUDeltaL2_ only if it doesn't exist
969  if (AU1TUDeltaL2_ == Teuchos::null)
970  AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
971  else {
972  if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
973  AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
974  }
975 
976  // Generate AU1TAU1_ only if it doesn't exist
977  if (AU1TAU1_ == Teuchos::null)
978  AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
979  else {
980  if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
981  AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
982  }
983 
984  // Generate GY_ only if it doesn't exist
985  if (GY_ == Teuchos::null)
986  GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
987  else {
988  if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
989  GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
990  }
991 
992  // Generate AU1TU1_ only if it doesn't exist
993  if (AU1TU1_ == Teuchos::null)
994  AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
995  else {
996  if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
997  AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
998  }
999 
1000  // Generate FY_ only if it doesn't exist
1001  if (FY_ == Teuchos::null)
1002  FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
1003  else {
1004  if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1005  FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1006  }
1007 
1008  // Generate AU1TAP_ only if it doesn't exist
1009  if (AU1TAP_ == Teuchos::null)
1010  AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1011  else {
1012  if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1013  AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1014  }
1015 
1016  // Generate APTAP_ only if it doesn't exist
1017  if (APTAP_ == Teuchos::null)
1018  APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1019  else {
1020  if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1021  APTAP_->reshape( numBlocks_, numBlocks_ );
1022  }
1023 
1024  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1025  if (U1Y1_ == Teuchos::null) {
1026  U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1027  }
1028  else {
1029  // Generate U1Y1_ by cloning itself ONLY if more space is needed.
1030  if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1031  Teuchos::RCP<const MV> tmp = U1Y1_;
1032  U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1033  }
1034  }
1035 
1036  // If the subspace has not been initialized before, generate it using the RHS from lp_.
1037  if (PY2_ == Teuchos::null) {
1038  PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1039  }
1040  else {
1041  // Generate PY2_ by cloning itself ONLY if more space is needed.
1042  if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1043  Teuchos::RCP<const MV> tmp = PY2_;
1044  PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1045  }
1046  }
1047 
1048  // Generate AUTAP_ only if it doesn't exist
1049  if (AUTAP_ == Teuchos::null)
1050  AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1051  else {
1052  if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1053  AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1054  }
1055 
1056  // Generate AU1TU_ only if it doesn't exist
1057  if (AU1TU_ == Teuchos::null)
1058  AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1059  else {
1060  if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1061  AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1062  }
1063 
1064 
1065  }
1066 }
1067 
1068 template<class ScalarType, class MV, class OP>
1070 
1071  Teuchos::LAPACK<int,ScalarType> lapack;
1072  std::vector<int> index(recycleBlocks_);
1073  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1074  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1075 
1076  // Count of number of cycles performed on current rhs
1077  int cycle = 0;
1078 
1079  // Set the current parameters if they were not set before.
1080  // NOTE: This may occur if the user generated the solver manager with the default constructor and
1081  // then didn't set any parameters using setParameters().
1082  if (!params_Set_) {
1083  setParameters(Teuchos::parameterList(*getValidParameters()));
1084  }
1085 
1086  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
1087  "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1088  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
1089  "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1090  TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1092  "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1093 
1094  // Grab the preconditioning object
1095  Teuchos::RCP<OP> precObj;
1096  if (problem_->getLeftPrec() != Teuchos::null) {
1097  precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1098  }
1099  else if (problem_->getRightPrec() != Teuchos::null) {
1100  precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1101  }
1102 
1103  // Create indices for the linear systems to be solved.
1104  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1105  std::vector<int> currIdx(1);
1106  currIdx[0] = 0;
1107 
1108  // Inform the linear problem of the current linear system to solve.
1109  problem_->setLSIndex( currIdx );
1110 
1111  // Check the number of blocks and change them if necessary.
1112  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1113  if (numBlocks_ > dim) {
1114  numBlocks_ = Teuchos::asSafe<int>(dim);
1115  params_->set("Num Blocks", numBlocks_);
1116  printer_->stream(Warnings) <<
1117  "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1118  " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1119  }
1120 
1121  // Initialize storage for all state variables
1122  initializeStateStorage();
1123 
1124  // Parameter list
1125  Teuchos::ParameterList plist;
1126  plist.set("Num Blocks",numBlocks_);
1127  plist.set("Recycled Blocks",recycleBlocks_);
1128 
1129  // Reset the status test.
1130  outputTest_->reset();
1131 
1132  // Assume convergence is achieved, then let any failed convergence set this to false.
1133  bool isConverged = true;
1134 
1135  // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
1136  if (existU_) {
1137  index.resize(recycleBlocks_);
1138  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1139  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1140  index.resize(recycleBlocks_);
1141  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1142  Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1143  // Initialize AU
1144  problem_->applyOp( *Utmp, *AUtmp );
1145  // Initialize UTAU
1146  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1147  MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1148  // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1149  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1150  if ( precObj != Teuchos::null ) {
1151  index.resize(recycleBlocks_);
1152  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1153  index.resize(recycleBlocks_);
1154  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1155  Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1156  OPT::Apply( *precObj, *AUtmp, *PCAU );
1157  MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1158  } else {
1159  MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1160  }
1161  }
1162 
1164  // RCG solver
1165 
1166  Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1167  rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
1168 
1169  // Enter solve() iterations
1170  {
1171 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1172  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1173 #endif
1174 
1175  while ( numRHS2Solve > 0 ) {
1176 
1177  // Debugging output to tell use if recycle space exists and will be used
1178  if (printer_->isVerbosity( Debug ) ) {
1179  if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
1180  else printer_->print( Debug, "No recycle space exists." );
1181  }
1182 
1183  // Reset the number of iterations.
1184  rcg_iter->resetNumIters();
1185 
1186  // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
1187  rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1188 
1189  // Reset the number of calls that the status test output knows about.
1190  outputTest_->resetNumCalls();
1191 
1192  // indicate that updated recycle space has not yet been generated for this linear system
1193  existU1_ = false;
1194 
1195  // reset cycle count
1196  cycle = 0;
1197 
1198  // Get the current residual
1199  problem_->computeCurrResVec( &*r_ );
1200 
1201  // If U exists, find best soln over this space first
1202  if (existU_) {
1203  // Solve linear system UTAU * y = (U'*r)
1204  Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1205  index.resize(recycleBlocks_);
1206  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1207  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1208  MVT::MvTransMv( one, *Utmp, *r_, Utr );
1209  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1210  Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1211  LUUTAUtmp.assign(UTAUtmp);
1212  int info = 0;
1213  lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1214  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1215  "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1216 
1217  // Update solution (x = x + U*y)
1218  MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1219 
1220  // Update residual ( r = r - AU*y )
1221  index.resize(recycleBlocks_);
1222  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1223  Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1224  MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1225  }
1226 
1227  if ( precObj != Teuchos::null ) {
1228  OPT::Apply( *precObj, *r_, *z_ );
1229  } else {
1230  z_ = r_;
1231  }
1232 
1233  // rTz_old = r'*z
1234  MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1235 
1236  if ( existU_ ) {
1237  // mu = UTAU\(AU'*z);
1238  Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1239  index.resize(recycleBlocks_);
1240  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1241  Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1242  MVT::MvTransMv( one, *AUtmp, *z_, mu );
1243  Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1244  char TRANS = 'N';
1245  int info;
1246  lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1247  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1248  "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1249  // p = z - U*mu;
1250  index.resize( 1 );
1251  index[0] = 0;
1252  Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1253  MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1254  MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1255  } else {
1256  // p = z;
1257  index.resize( 1 );
1258  index[0] = 0;
1259  Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1260  MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1261  }
1262 
1263  // Set the new state and initialize the solver.
1264  RCGIterState<ScalarType,MV> newstate;
1265 
1266  // Create RCP views here
1267  index.resize( numBlocks_+1 );
1268  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1269  newstate.P = MVT::CloneViewNonConst( *P_, index );
1270  index.resize( recycleBlocks_ );
1271  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1272  newstate.U = MVT::CloneViewNonConst( *U_, index );
1273  index.resize( recycleBlocks_ );
1274  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1275  newstate.AU = MVT::CloneViewNonConst( *AU_, index );
1276  newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1277  newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1278  newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1279  newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1280  newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1281  // assign the rest of the values in the struct
1282  newstate.curDim = 1; // We have initialized the first search vector
1283  newstate.Ap = Ap_;
1284  newstate.r = r_;
1285  newstate.z = z_;
1286  newstate.existU = existU_;
1287  newstate.ipiv = ipiv_;
1288  newstate.rTz_old = rTz_old_;
1289 
1290  rcg_iter->initialize(newstate);
1291 
1292  while(1) {
1293 
1294  // tell rcg_iter to iterate
1295  try {
1296  rcg_iter->iterate();
1297 
1299  //
1300  // check convergence first
1301  //
1303  if ( convTest_->getStatus() == Passed ) {
1304  // We have convergence
1305  break; // break from while(1){rcg_iter->iterate()}
1306  }
1308  //
1309  // check for maximum iterations
1310  //
1312  else if ( maxIterTest_->getStatus() == Passed ) {
1313  // we don't have convergence
1314  isConverged = false;
1315  break; // break from while(1){rcg_iter->iterate()}
1316  }
1318  //
1319  // check if cycle complete; update for next cycle
1320  //
1322  else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1323  // index into P_ of last search vector generated this cycle
1324  int lastp = -1;
1325  // index into Beta_ of last entry generated this cycle
1326  int lastBeta = -1;
1327  if (recycleBlocks_ > 0) {
1328  if (!existU_) {
1329  if (cycle == 0) { // No U, no U1
1330 
1331  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1332  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1333  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1334  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1335  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1336  Ftmp.putScalar(zero);
1337  Gtmp.putScalar(zero);
1338  for (int ii=0;ii<numBlocks_;ii++) {
1339  Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1340  if (ii > 0) {
1341  Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1342  Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1343  }
1344  Ftmp(ii,ii) = Dtmp(ii,0);
1345  }
1346 
1347  // compute harmonic Ritz vectors
1348  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1349  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1350 
1351  // U1 = [P(:,1:end-1)*Y];
1352  index.resize( numBlocks_ );
1353  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1354  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1355  index.resize( recycleBlocks_ );
1356  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1357  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1358  MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1359 
1360  // Precompute some variables for next cycle
1361 
1362  // AU1TAU1 = Y'*G*Y;
1363  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1364  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1365  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1366  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1367 
1368  // AU1TU1 = Y'*F*Y;
1369  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1370  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1371  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1372  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1373 
1374  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1375  // Must reinitialize AU1TAP; can become dense later
1376  AU1TAPtmp.putScalar(zero);
1377  // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1378  ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1379  for (int ii=0; ii<recycleBlocks_; ++ii) {
1380  AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1381  }
1382 
1383  // indicate that updated recycle space now defined
1384  existU1_ = true;
1385 
1386  // Indicate the size of the P, Beta structures generated this cycle
1387  lastp = numBlocks_;
1388  lastBeta = numBlocks_-1;
1389 
1390  } // if (cycle == 0)
1391  else { // No U, but U1 guaranteed to exist now
1392 
1393  // Finish computation of subblocks
1394  // AU1TAP = AU1TAP * D(1);
1395  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1396  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1397  AU1TAPtmp.scale(Dtmp(0,0));
1398 
1399  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1400  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1401  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1402  APTAPtmp.putScalar(zero);
1403  for (int ii=0; ii<numBlocks_; ii++) {
1404  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1405  if (ii > 0) {
1406  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1407  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1408  }
1409  }
1410 
1411  // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1412  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1413  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1414  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1415  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1416  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1417  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1418  F11.assign(AU1TU1tmp);
1419  F12.putScalar(zero);
1420  F21.putScalar(zero);
1421  F22.putScalar(zero);
1422  for(int ii=0;ii<numBlocks_;ii++) {
1423  F22(ii,ii) = Dtmp(ii,0);
1424  }
1425 
1426  // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1427  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1428  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1429  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1430  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1431  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1432  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1433  G11.assign(AU1TAU1tmp);
1434  G12.assign(AU1TAPtmp);
1435  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1436  for (int ii=0;ii<recycleBlocks_;++ii)
1437  for (int jj=0;jj<numBlocks_;++jj)
1438  G21(jj,ii) = G12(ii,jj);
1439  G22.assign(APTAPtmp);
1440 
1441  // compute harmonic Ritz vectors
1442  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1443  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1444 
1445  // U1 = [U1 P(:,2:end-1)]*Y;
1446  index.resize( numBlocks_ );
1447  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1448  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1449  index.resize( recycleBlocks_ );
1450  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1451  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1452  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1453  index.resize( recycleBlocks_ );
1454  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1455  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1456  index.resize( recycleBlocks_ );
1457  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1458  Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1459  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1460  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1461  MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1462  MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1463 
1464  // Precompute some variables for next cycle
1465 
1466  // AU1TAU1 = Y'*G*Y;
1467  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1468  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1469  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1470 
1471  // AU1TAP = zeros(k,m);
1472  // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1473  AU1TAPtmp.putScalar(zero);
1474  ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1475  for (int ii=0; ii<recycleBlocks_; ++ii) {
1476  AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1477  }
1478 
1479  // AU1TU1 = Y'*F*Y;
1480  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1481  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1482  //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1483  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1484 
1485  // Indicate the size of the P, Beta structures generated this cycle
1486  lastp = numBlocks_+1;
1487  lastBeta = numBlocks_;
1488 
1489  } // if (cycle != 1)
1490  } // if (!existU_)
1491  else { // U exists
1492  if (cycle == 0) { // No U1, but U exists
1493  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1494  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1495  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1496  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1497  APTAPtmp.putScalar(zero);
1498  for (int ii=0; ii<numBlocks_; ii++) {
1499  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1500  if (ii > 0) {
1501  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1502  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1503  }
1504  }
1505 
1506  Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1507  L2tmp.putScalar(zero);
1508  for(int ii=0;ii<numBlocks_;ii++) {
1509  L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1510  L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1511  }
1512 
1513  // AUTAP = UTAU*Delta*L2;
1514  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1515  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1516  Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1517  Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1518  DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1519  AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1520 
1521  // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
1522  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1523  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1524  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1525  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1526  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1527  F11.assign(UTAUtmp);
1528  F12.putScalar(zero);
1529  F21.putScalar(zero);
1530  F22.putScalar(zero);
1531  for(int ii=0;ii<numBlocks_;ii++) {
1532  F22(ii,ii) = Dtmp(ii,0);
1533  }
1534 
1535  // G = [AUTAU AUTAP; AUTAP' APTAP];
1536  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1537  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1538  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1539  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1540  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1541  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1542  G11.assign(AUTAUtmp);
1543  G12.assign(AUTAPtmp);
1544  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1545  for (int ii=0;ii<recycleBlocks_;++ii)
1546  for (int jj=0;jj<numBlocks_;++jj)
1547  G21(jj,ii) = G12(ii,jj);
1548  G22.assign(APTAPtmp);
1549 
1550  // compute harmonic Ritz vectors
1551  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1552  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1553 
1554  // U1 = [U P(:,1:end-1)]*Y;
1555  index.resize( recycleBlocks_ );
1556  for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1557  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1558  index.resize( numBlocks_ );
1559  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1560  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1561  index.resize( recycleBlocks_ );
1562  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1563  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1564  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1565  index.resize( recycleBlocks_ );
1566  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1567  Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1568  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1569  index.resize( recycleBlocks_ );
1570  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1571  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1572  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1573  MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1574  MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1575 
1576  // Precompute some variables for next cycle
1577 
1578  // AU1TAU1 = Y'*G*Y;
1579  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1580  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1581  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1582  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1583 
1584  // AU1TU1 = Y'*F*Y;
1585  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1586  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1587  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1588  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1589 
1590  // AU1TU = UTAU;
1591  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1592  AU1TUtmp.assign(UTAUtmp);
1593 
1594  // dold = D(end);
1595  dold = Dtmp(numBlocks_-1,0);
1596 
1597  // indicate that updated recycle space now defined
1598  existU1_ = true;
1599 
1600  // Indicate the size of the P, Beta structures generated this cycle
1601  lastp = numBlocks_;
1602  lastBeta = numBlocks_-1;
1603  }
1604  else { // Have U and U1
1605  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1606  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1607  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1608  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1609  for (int ii=0; ii<numBlocks_; ii++) {
1610  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1611  if (ii > 0) {
1612  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1613  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1614  }
1615  }
1616 
1617  Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1618  for(int ii=0;ii<numBlocks_;ii++) {
1619  L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1620  L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1621  }
1622 
1623  // M(end,1) = dold*(-Beta(1)/Alpha(1));
1624  // AU1TAP = Y'*[AU1TU*Delta*L2; M];
1625  Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1626  Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1627  DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1628  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1629  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1630  AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1631  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1632  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1633  AU1TAPtmp.putScalar(zero);
1634  AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1635  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1636  ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1637  for(int ii=0;ii<recycleBlocks_;ii++) {
1638  AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1639  }
1640 
1641  // AU1TU = Y1'*AU1TU
1642  Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1643  Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1644  AU1TUtmp.assign(Y1TAU1TU);
1645 
1646  // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1647  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1648  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1649  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1650  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1651  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1652  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1653  F11.assign(AU1TU1tmp);
1654  for(int ii=0;ii<numBlocks_;ii++) {
1655  F22(ii,ii) = Dtmp(ii,0);
1656  }
1657 
1658  // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1659  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1660  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1661  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1662  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1663  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1664  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1665  G11.assign(AU1TAU1tmp);
1666  G12.assign(AU1TAPtmp);
1667  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1668  for (int ii=0;ii<recycleBlocks_;++ii)
1669  for (int jj=0;jj<numBlocks_;++jj)
1670  G21(jj,ii) = G12(ii,jj);
1671  G22.assign(APTAPtmp);
1672 
1673  // compute harmonic Ritz vectors
1674  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1675  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1676 
1677  // U1 = [U1 P(:,2:end-1)]*Y;
1678  index.resize( numBlocks_ );
1679  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1680  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1681  index.resize( recycleBlocks_ );
1682  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1683  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1684  index.resize( recycleBlocks_ );
1685  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1686  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1687  index.resize( recycleBlocks_ );
1688  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1689  Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1690  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1691  MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1692  MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1693 
1694  // Precompute some variables for next cycle
1695 
1696  // AU1TAU1 = Y'*G*Y;
1697  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1698  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1699  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1700 
1701  // AU1TU1 = Y'*F*Y;
1702  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1703  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1704  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1705 
1706  // dold = D(end);
1707  dold = Dtmp(numBlocks_-1,0);
1708 
1709  // Indicate the size of the P, Beta structures generated this cycle
1710  lastp = numBlocks_+1;
1711  lastBeta = numBlocks_;
1712 
1713  }
1714  }
1715  } // if (recycleBlocks_ > 0)
1716 
1717  // Cleanup after end of cycle
1718 
1719  // P = P(:,end-1:end);
1720  index.resize( 2 );
1721  index[0] = lastp-1; index[1] = lastp;
1722  Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1723  index[0] = 0; index[1] = 1;
1724  MVT::SetBlock(*Ptmp2,index,*P_);
1725 
1726  // Beta = Beta(end);
1727  (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1728 
1729  // Delta = Delta(:,end);
1730  if (existU_) { // Delta only initialized if U exists
1731  Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1732  Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1733  mu1.assign(mu2);
1734  }
1735 
1736  // Now reinitialize state variables for next cycle
1737  newstate.P = Teuchos::null;
1738  index.resize( numBlocks_+1 );
1739  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1740  newstate.P = MVT::CloneViewNonConst( *P_, index );
1741 
1742  newstate.Beta = Teuchos::null;
1743  newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1744 
1745  newstate.Delta = Teuchos::null;
1746  newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1747 
1748  newstate.curDim = 1; // We have initialized the first search vector
1749 
1750  // Pass to iteration object
1751  rcg_iter->initialize(newstate);
1752 
1753  // increment cycle count
1754  cycle = cycle + 1;
1755 
1756  }
1758  //
1759  // we returned from iterate(), but none of our status tests Passed.
1760  // something is wrong, and it is probably our fault.
1761  //
1763  else {
1764  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1765  "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1766  }
1767  }
1768  catch (const std::exception &e) {
1769  printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration "
1770  << rcg_iter->getNumIters() << std::endl
1771  << e.what() << std::endl;
1772  throw;
1773  }
1774  }
1775 
1776  // Inform the linear problem that we are finished with this block linear system.
1777  problem_->setCurrLS();
1778 
1779  // Update indices for the linear systems to be solved.
1780  numRHS2Solve -= 1;
1781  if ( numRHS2Solve > 0 ) {
1782  currIdx[0]++;
1783  // Set the next indices.
1784  problem_->setLSIndex( currIdx );
1785  }
1786  else {
1787  currIdx.resize( numRHS2Solve );
1788  }
1789 
1790  // Update the recycle space for the next linear system
1791  if (existU1_) { // be sure updated recycle space was created
1792  // U = U1
1793  index.resize(recycleBlocks_);
1794  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1795  MVT::SetBlock(*U1_,index,*U_);
1796  // Set flag indicating recycle space is now defined
1797  existU_ = true;
1798  if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
1799  // Free pointers in newstate
1800  newstate.P = Teuchos::null;
1801  newstate.Ap = Teuchos::null;
1802  newstate.r = Teuchos::null;
1803  newstate.z = Teuchos::null;
1804  newstate.U = Teuchos::null;
1805  newstate.AU = Teuchos::null;
1806  newstate.Alpha = Teuchos::null;
1807  newstate.Beta = Teuchos::null;
1808  newstate.D = Teuchos::null;
1809  newstate.Delta = Teuchos::null;
1810  newstate.LUUTAU = Teuchos::null;
1811  newstate.ipiv = Teuchos::null;
1812  newstate.rTz_old = Teuchos::null;
1813 
1814  // Reinitialize AU, UTAU, AUTAU
1815  index.resize(recycleBlocks_);
1816  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1817  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1818  index.resize(recycleBlocks_);
1819  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1820  Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1821  // Initialize AU
1822  problem_->applyOp( *Utmp, *AUtmp );
1823  // Initialize UTAU
1824  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1825  MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1826  // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1827  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1828  if ( precObj != Teuchos::null ) {
1829  index.resize(recycleBlocks_);
1830  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1831  index.resize(recycleBlocks_);
1832  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1833  Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1834  OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1835  MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1836  } else {
1837  MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1838  }
1839  } // if (numRHS2Solve > 0)
1840 
1841  } // if (existU1)
1842  }// while ( numRHS2Solve > 0 )
1843 
1844  }
1845 
1846  // print final summary
1847  sTest_->print( printer_->stream(FinalSummary) );
1848 
1849  // print timing information
1850 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1851  // Calling summarize() can be expensive, so don't call unless the
1852  // user wants to print out timing details. summarize() will do all
1853  // the work even if it's passed a "black hole" output stream.
1854  if (verbosity_ & TimingDetails)
1855  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1856 #endif
1857 
1858  // get iteration information for this solve
1859  numIters_ = maxIterTest_->getNumIters();
1860 
1861  // Save the convergence test value ("achieved tolerance") for this solve.
1862  {
1863  using Teuchos::rcp_dynamic_cast;
1864  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1865  // testValues is nonnull and not persistent.
1866  const std::vector<MagnitudeType>* pTestValues =
1867  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1868 
1869  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1870  "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1871  "method returned NULL. Please report this bug to the Belos developers.");
1872 
1873  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1874  "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1875  "method returned a vector of length zero. Please report this bug to the "
1876  "Belos developers.");
1877 
1878  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1879  // achieved tolerances for all vectors in the current solve(), or
1880  // just for the vectors from the last deflation?
1881  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1882  }
1883 
1884  if (!isConverged) {
1885  return Unconverged; // return from RCGSolMgr::solve()
1886  }
1887  return Converged; // return from RCGSolMgr::solve()
1888 }
1889 
1890 // Compute the harmonic eigenpairs of the projected, dense system.
1891 template<class ScalarType, class MV, class OP>
1892 void RCGSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F,
1893  const Teuchos::SerialDenseMatrix<int,ScalarType>& /* G */,
1894  Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1895  // order of F,G
1896  int n = F.numCols();
1897 
1898  // The LAPACK interface
1899  Teuchos::LAPACK<int,ScalarType> lapack;
1900 
1901  // Magnitude of harmonic Ritz values
1902  std::vector<MagnitudeType> w(n);
1903 
1904  // Sorted order of harmonic Ritz values
1905  std::vector<int> iperm(n);
1906 
1907  // Compute k smallest harmonic Ritz pairs
1908  // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
1909  int itype = 1; // solve A*x = (lambda)*B*x
1910  char jobz='V'; // compute eigenvalues and eigenvectors
1911  char uplo='U'; // since F,G symmetric, reference only their upper triangular data
1912  std::vector<ScalarType> work(1);
1913  int lwork = -1;
1914  int info = 0;
1915  // since SYGV destroys workspace, create copies of F,G
1916  Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1917  Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1918 
1919  // query for optimal workspace size
1920  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1921  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1922  "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1923  lwork = (int)work[0];
1924  work.resize(lwork);
1925  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1926  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1927  "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1928 
1929 
1930  // Construct magnitude of each harmonic Ritz value
1931  this->sort(w,n,iperm);
1932 
1933  // Select recycledBlocks_ smallest eigenvectors
1934  for( int i=0; i<recycleBlocks_; i++ ) {
1935  for( int j=0; j<n; j++ ) {
1936  Y(j,i) = G2(j,iperm[i]);
1937  }
1938  }
1939 
1940 }
1941 
1942 // This method sorts list of n floating-point numbers and return permutation vector
1943 template<class ScalarType, class MV, class OP>
1944 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
1945 {
1946  int l, r, j, i, flag;
1947  int RR2;
1948  double dRR, dK;
1949 
1950  // Initialize the permutation vector.
1951  for(j=0;j<n;j++)
1952  iperm[j] = j;
1953 
1954  if (n <= 1) return;
1955 
1956  l = n / 2 + 1;
1957  r = n - 1;
1958  l = l - 1;
1959  dRR = dlist[l - 1];
1960  dK = dlist[l - 1];
1961 
1962  RR2 = iperm[l - 1];
1963  while (r != 0) {
1964  j = l;
1965  flag = 1;
1966 
1967  while (flag == 1) {
1968  i = j;
1969  j = j + j;
1970 
1971  if (j > r + 1)
1972  flag = 0;
1973  else {
1974  if (j < r + 1)
1975  if (dlist[j] > dlist[j - 1]) j = j + 1;
1976 
1977  if (dlist[j - 1] > dK) {
1978  dlist[i - 1] = dlist[j - 1];
1979  iperm[i - 1] = iperm[j - 1];
1980  }
1981  else {
1982  flag = 0;
1983  }
1984  }
1985  }
1986  dlist[i - 1] = dRR;
1987  iperm[i - 1] = RR2;
1988  if (l == 1) {
1989  dRR = dlist [r];
1990  RR2 = iperm[r];
1991  dK = dlist[r];
1992  dlist[r] = dlist[0];
1993  iperm[r] = iperm[0];
1994  r = r - 1;
1995  }
1996  else {
1997  l = l - 1;
1998  dRR = dlist[l - 1];
1999  RR2 = iperm[l - 1];
2000  dK = dlist[l - 1];
2001  }
2002  }
2003  dlist[0] = dRR;
2004  iperm[0] = RR2;
2005 }
2006 
2007 // This method requires the solver manager to return a std::string that describes itself.
2008 template<class ScalarType, class MV, class OP>
2010 {
2011  std::ostringstream oss;
2012  oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
2013  return oss.str();
2014 }
2015 
2016 } // end Belos namespace
2017 
2018 #endif /* BELOS_RCG_SOLMGR_HPP */
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
Teuchos::RCP< MV > AU
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
A factory class for generating StatusTestOutput objects.
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Belos::StatusTest class for specifying a maximum number of iterations.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
A factory class for generating StatusTestOutput objects.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
bool existU
Flag to indicate the recycle space should be used.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Pure virtual base class which describes the basic interface for a solver manager. ...
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A class for extending the status testing capabilities of Belos via logical combinations.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
RCGSolMgrRecyclingFailure(const std::string &what_arg)
Belos header file which uses auto-configuration information to include necessary C++ headers...
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...

Generated for Belos by doxygen 1.8.14