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