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