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