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

Generated on Mon Jul 15 2024 09:24:24 for Belos by doxygen 1.8.5