Belos Package Browser (Single Doxygen Collection)  Development
BelosBlockCGSolMgr.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_BLOCK_CG_SOLMGR_HPP
43 #define BELOS_BLOCK_CG_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosCGIter.hpp"
56 #include "BelosBlockCGIter.hpp"
62 #include "BelosStatusTestCombo.hpp"
64 #include "BelosOutputManager.hpp"
65 #include "Teuchos_BLAS.hpp"
66 #include "Teuchos_LAPACK.hpp"
67 #ifdef BELOS_TEUCHOS_TIME_MONITOR
68 # include "Teuchos_TimeMonitor.hpp"
69 #endif
70 #if defined(HAVE_TEUCHOSCORE_CXX11)
71 # include <type_traits>
72 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
73 #include <algorithm>
74 
91 namespace Belos {
92 
94 
95 
103  BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
104  {}};
105 
112  class BlockCGSolMgrOrthoFailure : public BelosError {public:
113  BlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
114  {}};
115 
116 
117  template<class ScalarType, class MV, class OP,
118  const bool lapackSupportsScalarType =
121  public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
122  {
123  static const bool requiresLapack =
125  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
127 
128  public:
130  base_type ()
131  {}
134  base_type ()
135  {}
136  virtual ~BlockCGSolMgr () {}
137  };
138 
139 
140  // Partial specialization for ScalarType types for which
141  // Teuchos::LAPACK has a valid implementation. This contains the
142  // actual working implementation of BlockCGSolMgr.
143  template<class ScalarType, class MV, class OP>
144  class BlockCGSolMgr<ScalarType, MV, OP, true> :
145  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
146  {
147  // This partial specialization is already chosen for those scalar types
148  // that support lapack, so we don't need to have an additional compile-time
149  // check that the scalar type is float/double/complex.
150 // #if defined(HAVE_TEUCHOSCORE_CXX11)
151 // # if defined(HAVE_TEUCHOS_COMPLEX)
152 // static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
153 // std::is_same<ScalarType, std::complex<double> >::value ||
154 // std::is_same<ScalarType, float>::value ||
155 // std::is_same<ScalarType, double>::value,
156 // "Belos::GCRODRSolMgr: ScalarType must be one of the four "
157 // "types (S,D,C,Z) supported by LAPACK.");
158 // # else
159 // static_assert (std::is_same<ScalarType, float>::value ||
160 // std::is_same<ScalarType, double>::value,
161 // "Belos::GCRODRSolMgr: ScalarType must be float or double. "
162 // "Complex arithmetic support is currently disabled. To "
163 // "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
164 // # endif // defined(HAVE_TEUCHOS_COMPLEX)
165 // #endif // defined(HAVE_TEUCHOSCORE_CXX11)
166 
167  private:
173 
174  public:
175 
177 
178 
184  BlockCGSolMgr();
185 
222 
224  virtual ~BlockCGSolMgr() {};
226 
228 
229 
231  return *problem_;
232  }
233 
236  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
237 
241 
248  return Teuchos::tuple(timerSolve_);
249  }
250 
257  return achievedTol_;
258  }
259 
261  int getNumIters() const {
262  return numIters_;
263  }
264 
267  bool isLOADetected() const { return false; }
268 
270 
272 
273 
275  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
276 
278  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
279 
281 
283 
284 
288  void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
290 
292 
293 
311  ReturnType solve();
312 
314 
317 
319  std::string description() const;
320 
322 
323  private:
324 
327 
332 
338 
341 
344 
347 
350 
353 
354  //
355  // Default solver parameters.
356  //
359  static const int maxIters_default_;
360  static const bool adaptiveBlockSize_default_;
361  static const bool showMaxResNormOnly_default_;
362  static const int blockSize_default_;
363  static const int verbosity_default_;
364  static const int outputStyle_default_;
365  static const int outputFreq_default_;
366  static const std::string label_default_;
367  static const std::string orthoType_default_;
369 
370  //
371  // Current solver parameters and other values.
372  //
373 
376 
379 
386 
389 
392 
393  int blockSize_, verbosity_, outputStyle_, outputFreq_;
394  bool adaptiveBlockSize_, showMaxResNormOnly_;
395  std::string orthoType_;
396 
398  std::string label_;
399 
402 
404  bool isSet_;
405  };
406 
407 
408 // Default solver values.
409 template<class ScalarType, class MV, class OP>
411 
412 template<class ScalarType, class MV, class OP>
414 
415 template<class ScalarType, class MV, class OP>
417 
418 template<class ScalarType, class MV, class OP>
420 
421 template<class ScalarType, class MV, class OP>
423 
424 template<class ScalarType, class MV, class OP>
426 
427 template<class ScalarType, class MV, class OP>
429 
430 template<class ScalarType, class MV, class OP>
432 
433 template<class ScalarType, class MV, class OP>
435 
436 template<class ScalarType, class MV, class OP>
437 const std::string BlockCGSolMgr<ScalarType,MV,OP,true>::label_default_ = "Belos";
438 
439 template<class ScalarType, class MV, class OP>
441 
442 template<class ScalarType, class MV, class OP>
444 
445 
446 // Empty Constructor
447 template<class ScalarType, class MV, class OP>
449  outputStream_(outputStream_default_),
450  convtol_(convtol_default_),
451  orthoKappa_(orthoKappa_default_),
452  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
453  maxIters_(maxIters_default_),
454  numIters_(0),
455  blockSize_(blockSize_default_),
456  verbosity_(verbosity_default_),
457  outputStyle_(outputStyle_default_),
458  outputFreq_(outputFreq_default_),
459  adaptiveBlockSize_(adaptiveBlockSize_default_),
460  showMaxResNormOnly_(showMaxResNormOnly_default_),
461  orthoType_(orthoType_default_),
462  label_(label_default_),
463  isSet_(false)
464 {}
465 
466 
467 // Basic Constructor
468 template<class ScalarType, class MV, class OP>
472  problem_(problem),
473  outputStream_(outputStream_default_),
474  convtol_(convtol_default_),
475  orthoKappa_(orthoKappa_default_),
476  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
477  maxIters_(maxIters_default_),
478  numIters_(0),
479  blockSize_(blockSize_default_),
480  verbosity_(verbosity_default_),
481  outputStyle_(outputStyle_default_),
482  outputFreq_(outputFreq_default_),
483  adaptiveBlockSize_(adaptiveBlockSize_default_),
484  showMaxResNormOnly_(showMaxResNormOnly_default_),
485  orthoType_(orthoType_default_),
486  label_(label_default_),
487  isSet_(false)
488 {
489  TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
490  "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
491 
492  // If the user passed in a nonnull parameter list, set parameters.
493  // Otherwise, the next solve() call will use default parameters,
494  // unless the user calls setParameters() first.
495  if (! pl.is_null()) {
496  setParameters (pl);
497  }
498 }
499 
500 template<class ScalarType, class MV, class OP>
501 void
504 {
505  // Create the internal parameter list if one doesn't already exist.
506  if (params_ == Teuchos::null) {
507  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
508  }
509  else {
510  params->validateParameters(*getValidParameters());
511  }
512 
513  // Check for maximum number of iterations
514  if (params->isParameter("Maximum Iterations")) {
515  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
516 
517  // Update parameter in our list and in status test.
518  params_->set("Maximum Iterations", maxIters_);
519  if (maxIterTest_!=Teuchos::null)
520  maxIterTest_->setMaxIters( maxIters_ );
521  }
522 
523  // Check for blocksize
524  if (params->isParameter("Block Size")) {
525  blockSize_ = params->get("Block Size",blockSize_default_);
526  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
527  "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
528 
529  // Update parameter in our list.
530  params_->set("Block Size", blockSize_);
531  }
532 
533  // Check if the blocksize should be adaptive
534  if (params->isParameter("Adaptive Block Size")) {
535  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
536 
537  // Update parameter in our list.
538  params_->set("Adaptive Block Size", adaptiveBlockSize_);
539  }
540 
541  // Check to see if the timer label changed.
542  if (params->isParameter("Timer Label")) {
543  std::string tempLabel = params->get("Timer Label", label_default_);
544 
545  // Update parameter in our list and solver timer
546  if (tempLabel != label_) {
547  label_ = tempLabel;
548  params_->set("Timer Label", label_);
549  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
550 #ifdef BELOS_TEUCHOS_TIME_MONITOR
551  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
552 #endif
553  if (ortho_ != Teuchos::null) {
554  ortho_->setLabel( label_ );
555  }
556  }
557  }
558 
559  // Check if the orthogonalization changed.
560  if (params->isParameter("Orthogonalization")) {
561  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
562  TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
563  std::invalid_argument,
564  "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
565  if (tempOrthoType != orthoType_) {
566  orthoType_ = tempOrthoType;
567  // Create orthogonalization manager
568  if (orthoType_=="DGKS") {
569  if (orthoKappa_ <= 0) {
570  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
571  }
572  else {
573  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
574  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
575  }
576  }
577  else if (orthoType_=="ICGS") {
578  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
579  }
580  else if (orthoType_=="IMGS") {
581  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
582  }
583  }
584  }
585 
586  // Check which orthogonalization constant to use.
587  if (params->isParameter("Orthogonalization Constant")) {
588  orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
589 
590  // Update parameter in our list.
591  params_->set("Orthogonalization Constant",orthoKappa_);
592  if (orthoType_=="DGKS") {
593  if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
594  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
595  }
596  }
597  }
598 
599  // Check for a change in verbosity level
600  if (params->isParameter("Verbosity")) {
601  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
602  verbosity_ = params->get("Verbosity", verbosity_default_);
603  } else {
604  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
605  }
606 
607  // Update parameter in our list.
608  params_->set("Verbosity", verbosity_);
609  if (printer_ != Teuchos::null)
610  printer_->setVerbosity(verbosity_);
611  }
612 
613  // Check for a change in output style
614  if (params->isParameter("Output Style")) {
615  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
616  outputStyle_ = params->get("Output Style", outputStyle_default_);
617  } else {
618  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
619  }
620 
621  // Update parameter in our list.
622  params_->set("Output Style", outputStyle_);
623  outputTest_ = Teuchos::null;
624  }
625 
626  // output stream
627  if (params->isParameter("Output Stream")) {
628  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
629 
630  // Update parameter in our list.
631  params_->set("Output Stream", outputStream_);
632  if (printer_ != Teuchos::null)
633  printer_->setOStream( outputStream_ );
634  }
635 
636  // frequency level
637  if (verbosity_ & Belos::StatusTestDetails) {
638  if (params->isParameter("Output Frequency")) {
639  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
640  }
641 
642  // Update parameter in out list and output status test.
643  params_->set("Output Frequency", outputFreq_);
644  if (outputTest_ != Teuchos::null)
645  outputTest_->setOutputFrequency( outputFreq_ );
646  }
647 
648  // Create output manager if we need to.
649  if (printer_ == Teuchos::null) {
650  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
651  }
652 
653  // Convergence
654  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
655  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
656 
657  // Check for convergence tolerance
658  if (params->isParameter("Convergence Tolerance")) {
659  convtol_ = params->get("Convergence Tolerance",convtol_default_);
660 
661  // Update parameter in our list and residual tests.
662  params_->set("Convergence Tolerance", convtol_);
663  if (convTest_ != Teuchos::null)
664  convTest_->setTolerance( convtol_ );
665  }
666 
667  if (params->isParameter("Show Maximum Residual Norm Only")) {
668  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
669 
670  // Update parameter in our list and residual tests
671  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
672  if (convTest_ != Teuchos::null)
673  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
674  }
675 
676  // Create status tests if we need to.
677 
678  // Basic test checks maximum iterations and native residual.
679  if (maxIterTest_ == Teuchos::null)
680  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
681 
682  // Implicit residual test, using the native residual to determine if convergence was achieved.
683  if (convTest_ == Teuchos::null)
684  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
685 
686  if (sTest_ == Teuchos::null)
687  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
688 
689  if (outputTest_ == Teuchos::null) {
690 
691  // Create the status test output class.
692  // This class manages and formats the output from the status test.
693  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
694  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
695 
696  // Set the solver string for the output test
697  std::string solverDesc = " Block CG ";
698  outputTest_->setSolverDesc( solverDesc );
699 
700  }
701 
702  // Create orthogonalization manager if we need to.
703  if (ortho_ == Teuchos::null) {
704  if (orthoType_=="DGKS") {
705  if (orthoKappa_ <= 0) {
706  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
707  }
708  else {
709  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
710  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
711  }
712  }
713  else if (orthoType_=="ICGS") {
714  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
715  }
716  else if (orthoType_=="IMGS") {
717  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
718  }
719  else {
720  TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
721  "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
722  }
723  }
724 
725  // Create the timer if we need to.
726  if (timerSolve_ == Teuchos::null) {
727  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
728 #ifdef BELOS_TEUCHOS_TIME_MONITOR
729  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
730 #endif
731  }
732 
733  // Inform the solver manager that the current parameters were set.
734  isSet_ = true;
735 }
736 
737 
738 template<class ScalarType, class MV, class OP>
741 {
743 
744  // Set all the valid parameters and their default values.
745  if(is_null(validPL)) {
746  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
747  pl->set("Convergence Tolerance", convtol_default_,
748  "The relative residual tolerance that needs to be achieved by the\n"
749  "iterative solver in order for the linear system to be declared converged.");
750  pl->set("Maximum Iterations", maxIters_default_,
751  "The maximum number of block iterations allowed for each\n"
752  "set of RHS solved.");
753  pl->set("Block Size", blockSize_default_,
754  "The number of vectors in each block.");
755  pl->set("Adaptive Block Size", adaptiveBlockSize_default_,
756  "Whether the solver manager should adapt to the block size\n"
757  "based on the number of RHS to solve.");
758  pl->set("Verbosity", verbosity_default_,
759  "What type(s) of solver information should be outputted\n"
760  "to the output stream.");
761  pl->set("Output Style", outputStyle_default_,
762  "What style is used for the solver information outputted\n"
763  "to the output stream.");
764  pl->set("Output Frequency", outputFreq_default_,
765  "How often convergence information should be outputted\n"
766  "to the output stream.");
767  pl->set("Output Stream", outputStream_default_,
768  "A reference-counted pointer to the output stream where all\n"
769  "solver output is sent.");
770  pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
771  "When convergence information is printed, only show the maximum\n"
772  "relative residual norm when the block size is greater than one.");
773  pl->set("Timer Label", label_default_,
774  "The string to use as a prefix for the timer labels.");
775  // pl->set("Restart Timers", restartTimers_);
776  pl->set("Orthogonalization", orthoType_default_,
777  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
778  pl->set("Orthogonalization Constant",orthoKappa_default_,
779  "The constant used by DGKS orthogonalization to determine\n"
780  "whether another step of classical Gram-Schmidt is necessary.");
781  validPL = pl;
782  }
783  return validPL;
784 }
785 
786 
787 // solve()
788 template<class ScalarType, class MV, class OP>
790  using Teuchos::RCP;
791  using Teuchos::rcp;
792  using Teuchos::rcp_const_cast;
793  using Teuchos::rcp_dynamic_cast;
794 
795  // Set the current parameters if they were not set before. NOTE:
796  // This may occur if the user generated the solver manager with the
797  // default constructor and then didn't set any parameters using
798  // setParameters().
799  if (!isSet_) {
800  setParameters(Teuchos::parameterList(*getValidParameters()));
801  }
802 
805 
806  TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
808  "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
809  "has not been called.");
810 
811  // Create indices for the linear systems to be solved.
812  int startPtr = 0;
813  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
814  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
815 
816  std::vector<int> currIdx, currIdx2;
817  // If an adaptive block size is allowed then only the linear
818  // systems that need to be solved are solved. Otherwise, the index
819  // set is generated that informs the linear problem that some
820  // linear systems are augmented.
821  if ( adaptiveBlockSize_ ) {
822  blockSize_ = numCurrRHS;
823  currIdx.resize( numCurrRHS );
824  currIdx2.resize( numCurrRHS );
825  for (int i=0; i<numCurrRHS; ++i)
826  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
827 
828  }
829  else {
830  currIdx.resize( blockSize_ );
831  currIdx2.resize( blockSize_ );
832  for (int i=0; i<numCurrRHS; ++i)
833  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
834  for (int i=numCurrRHS; i<blockSize_; ++i)
835  { currIdx[i] = -1; currIdx2[i] = i; }
836  }
837 
838  // Inform the linear problem of the current linear system to solve.
839  problem_->setLSIndex( currIdx );
840 
842  // Set up the parameter list for the Iteration subclass.
844  plist.set("Block Size",blockSize_);
845 
846  // Reset the output status test (controls all the other status tests).
847  outputTest_->reset();
848 
849  // Assume convergence is achieved, then let any failed convergence
850  // set this to false. "Innocent until proven guilty."
851  bool isConverged = true;
852 
854  // Set up the BlockCG Iteration subclass.
855 
856  RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
857  if (blockSize_ == 1) {
858  // Standard (nonblock) CG is faster for the special case of a
859  // block size of 1.
860  block_cg_iter =
861  rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
862  outputTest_, plist));
863  } else {
864  block_cg_iter =
865  rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
866  ortho_, plist));
867  }
868 
869 
870  // Enter solve() iterations
871  {
872 #ifdef BELOS_TEUCHOS_TIME_MONITOR
873  Teuchos::TimeMonitor slvtimer(*timerSolve_);
874 #endif
875 
876  while ( numRHS2Solve > 0 ) {
877  //
878  // Reset the active / converged vectors from this block
879  std::vector<int> convRHSIdx;
880  std::vector<int> currRHSIdx( currIdx );
881  currRHSIdx.resize(numCurrRHS);
882 
883  // Reset the number of iterations.
884  block_cg_iter->resetNumIters();
885 
886  // Reset the number of calls that the status test output knows about.
887  outputTest_->resetNumCalls();
888 
889  // Get the current residual for this block of linear systems.
890  RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
891 
892  // Set the new state and initialize the solver.
894  newstate.R = R_0;
895  block_cg_iter->initializeCG(newstate);
896 
897  while(1) {
898 
899  // tell block_cg_iter to iterate
900  try {
901  block_cg_iter->iterate();
902  //
903  // Check whether any of the linear systems converged.
904  //
905  if (convTest_->getStatus() == Passed) {
906  // At least one of the linear system(s) converged.
907  //
908  // Get the column indices of the linear systems that converged.
909  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
910  std::vector<int> convIdx =
911  rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
912 
913  // If the number of converged linear systems equals the
914  // number of linear systems currently being solved, then
915  // we are done with this block.
916  if (convIdx.size() == currRHSIdx.size())
917  break; // break from while(1){block_cg_iter->iterate()}
918 
919  // Inform the linear problem that we are finished with
920  // this current linear system.
921  problem_->setCurrLS();
922 
923  // Reset currRHSIdx to contain the right-hand sides that
924  // are left to converge for this block.
925  int have = 0;
926  std::vector<int> unconvIdx(currRHSIdx.size());
927  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
928  bool found = false;
929  for (unsigned int j=0; j<convIdx.size(); ++j) {
930  if (currRHSIdx[i] == convIdx[j]) {
931  found = true;
932  break;
933  }
934  }
935  if (!found) {
936  currIdx2[have] = currIdx2[i];
937  currRHSIdx[have++] = currRHSIdx[i];
938  }
939  else {
940  }
941  }
942  currRHSIdx.resize(have);
943  currIdx2.resize(have);
944 
945  // Set the remaining indices after deflation.
946  problem_->setLSIndex( currRHSIdx );
947 
948  // Get the current residual vector.
949  std::vector<MagnitudeType> norms;
950  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
951  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
952 
953  // Set the new blocksize for the solver.
954  block_cg_iter->setBlockSize( have );
955 
956  // Set the new state and initialize the solver.
958  defstate.R = R_0;
959  block_cg_iter->initializeCG(defstate);
960  }
961  //
962  // None of the linear systems converged. Check whether the
963  // maximum iteration count was reached.
964  //
965  else if (maxIterTest_->getStatus() == Passed) {
966  isConverged = false; // None of the linear systems converged.
967  break; // break from while(1){block_cg_iter->iterate()}
968  }
969  //
970  // iterate() returned, but none of our status tests Passed.
971  // This indicates a bug.
972  //
973  else {
974  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
975  "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
976  "the maximum iteration count test passed. Please report this bug "
977  "to the Belos developers.");
978  }
979  }
980  catch (const std::exception &e) {
981  std::ostream& err = printer_->stream (Errors);
982  err << "Error! Caught std::exception in CGIteration::iterate() at "
983  << "iteration " << block_cg_iter->getNumIters() << std::endl
984  << e.what() << std::endl;
985  throw;
986  }
987  }
988 
989  // Inform the linear problem that we are finished with this
990  // block linear system.
991  problem_->setCurrLS();
992 
993  // Update indices for the linear systems to be solved.
994  startPtr += numCurrRHS;
995  numRHS2Solve -= numCurrRHS;
996  if ( numRHS2Solve > 0 ) {
997  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
998 
999  if ( adaptiveBlockSize_ ) {
1000  blockSize_ = numCurrRHS;
1001  currIdx.resize( numCurrRHS );
1002  currIdx2.resize( numCurrRHS );
1003  for (int i=0; i<numCurrRHS; ++i)
1004  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1005  }
1006  else {
1007  currIdx.resize( blockSize_ );
1008  currIdx2.resize( blockSize_ );
1009  for (int i=0; i<numCurrRHS; ++i)
1010  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1011  for (int i=numCurrRHS; i<blockSize_; ++i)
1012  { currIdx[i] = -1; currIdx2[i] = i; }
1013  }
1014  // Set the next indices.
1015  problem_->setLSIndex( currIdx );
1016 
1017  // Set the new blocksize for the solver.
1018  block_cg_iter->setBlockSize( blockSize_ );
1019  }
1020  else {
1021  currIdx.resize( numRHS2Solve );
1022  }
1023 
1024  }// while ( numRHS2Solve > 0 )
1025 
1026  }
1027 
1028  // print final summary
1029  sTest_->print( printer_->stream(FinalSummary) );
1030 
1031  // print timing information
1032 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1033  // Calling summarize() requires communication in general, so don't
1034  // call it unless the user wants to print out timing details.
1035  // summarize() will do all the work even if it's passed a "black
1036  // hole" output stream.
1037  if (verbosity_ & TimingDetails) {
1038  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1039  }
1040 #endif
1041 
1042  // Save the iteration count for this solve.
1043  numIters_ = maxIterTest_->getNumIters();
1044 
1045  // Save the convergence test value ("achieved tolerance") for this solve.
1046  {
1047  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1048  // testValues is nonnull and not persistent.
1049  const std::vector<MagnitudeType>* pTestValues =
1050  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1051 
1052  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1053  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1054  "method returned NULL. Please report this bug to the Belos developers.");
1055 
1056  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1057  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1058  "method returned a vector of length zero. Please report this bug to the "
1059  "Belos developers.");
1060 
1061  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1062  // achieved tolerances for all vectors in the current solve(), or
1063  // just for the vectors from the last deflation?
1064  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1065  }
1066 
1067  if (!isConverged) {
1068  return Unconverged; // return from BlockCGSolMgr::solve()
1069  }
1070  return Converged; // return from BlockCGSolMgr::solve()
1071 }
1072 
1073 // This method requires the solver manager to return a std::string that describes itself.
1074 template<class ScalarType, class MV, class OP>
1076 {
1077  std::ostringstream oss;
1078  oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1079  oss << "{";
1080  oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1081  oss << "}";
1082  return oss.str();
1083 }
1084 
1085 } // end Belos namespace
1086 
1087 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
Teuchos::RCP< const MV > R
The current residual.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Collection of types and exceptions used within the Belos solvers.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
int maxIters_
Maximum iteration count (read from parameter list).
Class which manages the output and verbosity of the Belos solvers.
int getNumIters() const
Get the iteration count for the most recent call to solve().
bool is_null(const boost::shared_ptr< T > &p)
bool isSet_
Whether or not the parameters have been set (via setParameters()).
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Output "status test" that controls all the other status tests.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
T & get(ParameterList &l, const std::string &name)
Structure to contain pointers to CGIteration state variables.
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Belos concrete class for performing the block conjugate-gradient (CG) iteration.
An implementation of StatusTestResNorm using a family of residual norms.
std::string label_
Prefix label for all the timers.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
BlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
BlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Maximum iteration count stopping criterion.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Aggregate stopping criterion.
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
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)
bool is_null(const RCP< T > &p)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
BlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonor...
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Convergence stopping criterion.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
BlockCGSolMgrOrthoFailure(const std::string &what_arg)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
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.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
Teuchos::RCP< Teuchos::Time > timerSolve_
Solve timer.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Teuchos::RCP< OutputManager< ScalarType > > printer_
Output manager, that handles printing of different kinds of messages.
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
int numIters_
Number of iterations taken by the last solve() invocation.
static const Teuchos::RCP< std::ostream > outputStream_default_
Teuchos::ScalarTraits< MagnitudeType > MT
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
static const bool requiresLapack
Teuchos::RCP< std::ostream > outputStream_
Output stream to which the output manager prints.
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...