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