Belos Package Browser (Single Doxygen Collection) Development
BelosBlockCGSolMgr.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_BLOCK_CG_SOLMGR_HPP
00043 #define BELOS_BLOCK_CG_SOLMGR_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosCGIter.hpp"
00056 #include "BelosBlockCGIter.hpp"
00057 #include "BelosDGKSOrthoManager.hpp"
00058 #include "BelosICGSOrthoManager.hpp"
00059 #include "BelosIMGSOrthoManager.hpp"
00060 #include "BelosStatusTestMaxIters.hpp"
00061 #include "BelosStatusTestGenResNorm.hpp"
00062 #include "BelosStatusTestCombo.hpp"
00063 #include "BelosStatusTestOutputFactory.hpp"
00064 #include "BelosOutputManager.hpp"
00065 #include "Teuchos_BLAS.hpp"
00066 #include "Teuchos_LAPACK.hpp"
00067 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00068 #  include "Teuchos_TimeMonitor.hpp"
00069 #endif
00070 #include <algorithm>
00071 
00088 namespace Belos {
00089   
00091 
00092   
00099   class BlockCGSolMgrLinearProblemFailure : public BelosError {public:
00100     BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00101     {}};
00102   
00109   class BlockCGSolMgrOrthoFailure : public BelosError {public:
00110     BlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00111     {}};
00112   
00113   template<class ScalarType, class MV, class OP>
00114   class BlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00115     
00116   private:
00117     typedef MultiVecTraits<ScalarType,MV> MVT;
00118     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00119     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00120     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00121     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00122     
00123   public:
00124     
00126 
00127 
00133      BlockCGSolMgr();
00134 
00163     BlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00164        const Teuchos::RCP<Teuchos::ParameterList> &pl );
00165     
00167     virtual ~BlockCGSolMgr() {};
00169     
00171 
00172     
00173     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00174       return *problem_;
00175     }
00176 
00179     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00180     
00183     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00184     
00190     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00191       return Teuchos::tuple(timerSolve_);
00192     }
00193 
00199     MagnitudeType achievedTol() const {
00200       return achievedTol_;
00201     }
00202 
00204     int getNumIters() const {
00205       return numIters_;
00206     }
00207 
00210     bool isLOADetected() const { return false; }
00211  
00213     
00215 
00216    
00218     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00219    
00221     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00222     
00224    
00226 
00227 
00231     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00233  
00235 
00236     
00254     ReturnType solve();
00255     
00257     
00260     
00262     std::string description() const;
00263     
00265     
00266   private:
00267 
00269     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00270     
00272     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00274     Teuchos::RCP<std::ostream> outputStream_;
00275 
00280     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00281 
00283     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00284 
00286     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00287 
00289     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00290 
00292     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00293     
00295     Teuchos::RCP<Teuchos::ParameterList> params_;
00296     
00297     //
00298     // Default solver parameters.
00299     //
00300     static const MagnitudeType convtol_default_;
00301     static const MagnitudeType orthoKappa_default_;
00302     static const int maxIters_default_;
00303     static const bool adaptiveBlockSize_default_;
00304     static const bool showMaxResNormOnly_default_;
00305     static const int blockSize_default_;
00306     static const int verbosity_default_;
00307     static const int outputStyle_default_;
00308     static const int outputFreq_default_;
00309     static const std::string label_default_;
00310     static const std::string orthoType_default_;
00311     static const Teuchos::RCP<std::ostream> outputStream_default_;
00312 
00313     //
00314     // Current solver parameters and other values.
00315     //
00316 
00318     MagnitudeType convtol_;
00319 
00321     MagnitudeType orthoKappa_;
00322 
00328     MagnitudeType achievedTol_;
00329 
00331     int maxIters_;
00332 
00334     int numIters_;
00335 
00336     int blockSize_, verbosity_, outputStyle_, outputFreq_;
00337     bool adaptiveBlockSize_, showMaxResNormOnly_;
00338     std::string orthoType_; 
00339     
00341     std::string label_;
00342 
00344     Teuchos::RCP<Teuchos::Time> timerSolve_;
00345 
00347     bool isSet_;
00348   };
00349 
00350 
00351 // Default solver values.
00352 template<class ScalarType, class MV, class OP>
00353 const typename BlockCGSolMgr<ScalarType,MV,OP>::MagnitudeType BlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00354 
00355 template<class ScalarType, class MV, class OP>
00356 const typename BlockCGSolMgr<ScalarType,MV,OP>::MagnitudeType BlockCGSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00357 
00358 template<class ScalarType, class MV, class OP>
00359 const int BlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00360 
00361 template<class ScalarType, class MV, class OP>
00362 const bool BlockCGSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ = true;
00363 
00364 template<class ScalarType, class MV, class OP>
00365 const bool BlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00366 
00367 template<class ScalarType, class MV, class OP>
00368 const int BlockCGSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00369 
00370 template<class ScalarType, class MV, class OP>
00371 const int BlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00372 
00373 template<class ScalarType, class MV, class OP>
00374 const int BlockCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00375 
00376 template<class ScalarType, class MV, class OP>
00377 const int BlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00378 
00379 template<class ScalarType, class MV, class OP>
00380 const std::string BlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00381 
00382 template<class ScalarType, class MV, class OP>
00383 const std::string BlockCGSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00384 
00385 template<class ScalarType, class MV, class OP>
00386 const Teuchos::RCP<std::ostream> BlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00387 
00388 
00389 // Empty Constructor
00390 template<class ScalarType, class MV, class OP>
00391 BlockCGSolMgr<ScalarType,MV,OP>::BlockCGSolMgr() :
00392   outputStream_(outputStream_default_),
00393   convtol_(convtol_default_),
00394   orthoKappa_(orthoKappa_default_),
00395   achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00396   maxIters_(maxIters_default_),
00397   numIters_(0),
00398   blockSize_(blockSize_default_),
00399   verbosity_(verbosity_default_),
00400   outputStyle_(outputStyle_default_),
00401   outputFreq_(outputFreq_default_),
00402   adaptiveBlockSize_(adaptiveBlockSize_default_),
00403   showMaxResNormOnly_(showMaxResNormOnly_default_),
00404   orthoType_(orthoType_default_),
00405   label_(label_default_),
00406   isSet_(false)
00407 {}
00408 
00409 
00410 // Basic Constructor
00411 template<class ScalarType, class MV, class OP>
00412 BlockCGSolMgr<ScalarType,MV,OP>::
00413 BlockCGSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00414         const Teuchos::RCP<Teuchos::ParameterList> &pl) : 
00415   problem_(problem),
00416   outputStream_(outputStream_default_),
00417   convtol_(convtol_default_),
00418   orthoKappa_(orthoKappa_default_),
00419   achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00420   maxIters_(maxIters_default_),
00421   numIters_(0),
00422   blockSize_(blockSize_default_),
00423   verbosity_(verbosity_default_),
00424   outputStyle_(outputStyle_default_),
00425   outputFreq_(outputFreq_default_),
00426   adaptiveBlockSize_(adaptiveBlockSize_default_),
00427   showMaxResNormOnly_(showMaxResNormOnly_default_),
00428   orthoType_(orthoType_default_),
00429   label_(label_default_),
00430   isSet_(false)
00431 {
00432   TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument, 
00433     "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
00434 
00435   // If the user passed in a nonnull parameter list, set parameters.
00436   // Otherwise, the next solve() call will use default parameters,
00437   // unless the user calls setParameters() first.
00438   if (! pl.is_null()) {
00439     setParameters (pl);  
00440   }
00441 }
00442 
00443 template<class ScalarType, class MV, class OP>
00444 void 
00445 BlockCGSolMgr<ScalarType,MV,OP>::
00446 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00447 {
00448   // Create the internal parameter list if one doesn't already exist.
00449   if (params_ == Teuchos::null) {
00450     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00451   }
00452   else {
00453     params->validateParameters(*getValidParameters());
00454   }
00455 
00456   // Check for maximum number of iterations
00457   if (params->isParameter("Maximum Iterations")) {
00458     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00459 
00460     // Update parameter in our list and in status test.
00461     params_->set("Maximum Iterations", maxIters_);
00462     if (maxIterTest_!=Teuchos::null)
00463       maxIterTest_->setMaxIters( maxIters_ );
00464   }
00465 
00466   // Check for blocksize
00467   if (params->isParameter("Block Size")) {
00468     blockSize_ = params->get("Block Size",blockSize_default_);    
00469     TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00470            "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
00471 
00472     // Update parameter in our list.
00473     params_->set("Block Size", blockSize_);
00474   }
00475 
00476   // Check if the blocksize should be adaptive
00477   if (params->isParameter("Adaptive Block Size")) {
00478     adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
00479     
00480     // Update parameter in our list.
00481     params_->set("Adaptive Block Size", adaptiveBlockSize_);
00482   }
00483 
00484   // Check to see if the timer label changed.
00485   if (params->isParameter("Timer Label")) {
00486     std::string tempLabel = params->get("Timer Label", label_default_);
00487 
00488     // Update parameter in our list and solver timer
00489     if (tempLabel != label_) {
00490       label_ = tempLabel;
00491       params_->set("Timer Label", label_);
00492       std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
00493 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00494       timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00495 #endif
00496       if (ortho_ != Teuchos::null) {
00497         ortho_->setLabel( label_ );
00498       }
00499     }
00500   }
00501 
00502   // Check if the orthogonalization changed.
00503   if (params->isParameter("Orthogonalization")) {
00504     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00505     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00506                         std::invalid_argument,
00507       "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00508     if (tempOrthoType != orthoType_) {
00509       orthoType_ = tempOrthoType;
00510       // Create orthogonalization manager
00511       if (orthoType_=="DGKS") {
00512   if (orthoKappa_ <= 0) {
00513     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00514   }
00515   else {
00516     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00517     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00518   }
00519       }
00520       else if (orthoType_=="ICGS") {
00521   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00522       } 
00523       else if (orthoType_=="IMGS") {
00524   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00525       } 
00526     }  
00527   }
00528 
00529   // Check which orthogonalization constant to use.
00530   if (params->isParameter("Orthogonalization Constant")) {
00531     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00532 
00533     // Update parameter in our list.
00534     params_->set("Orthogonalization Constant",orthoKappa_);
00535     if (orthoType_=="DGKS") {
00536       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00537   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00538       }
00539     } 
00540   }
00541 
00542   // Check for a change in verbosity level
00543   if (params->isParameter("Verbosity")) {
00544     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00545       verbosity_ = params->get("Verbosity", verbosity_default_);
00546     } else {
00547       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00548     }
00549 
00550     // Update parameter in our list.
00551     params_->set("Verbosity", verbosity_);
00552     if (printer_ != Teuchos::null)
00553       printer_->setVerbosity(verbosity_);
00554   }
00555 
00556   // Check for a change in output style
00557   if (params->isParameter("Output Style")) {
00558     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00559       outputStyle_ = params->get("Output Style", outputStyle_default_);
00560     } else {
00561       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00562     }
00563 
00564     // Update parameter in our list.
00565     params_->set("Output Style", outputStyle_);
00566     outputTest_ = Teuchos::null;
00567   }
00568 
00569   // output stream
00570   if (params->isParameter("Output Stream")) {
00571     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00572 
00573     // Update parameter in our list.
00574     params_->set("Output Stream", outputStream_);
00575     if (printer_ != Teuchos::null)
00576       printer_->setOStream( outputStream_ );
00577   }
00578 
00579   // frequency level
00580   if (verbosity_ & Belos::StatusTestDetails) {
00581     if (params->isParameter("Output Frequency")) {
00582       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00583     }
00584 
00585     // Update parameter in out list and output status test.
00586     params_->set("Output Frequency", outputFreq_);
00587     if (outputTest_ != Teuchos::null)
00588       outputTest_->setOutputFrequency( outputFreq_ );
00589   }
00590 
00591   // Create output manager if we need to.
00592   if (printer_ == Teuchos::null) {
00593     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00594   }  
00595   
00596   // Convergence
00597   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00598   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00599 
00600   // Check for convergence tolerance
00601   if (params->isParameter("Convergence Tolerance")) {
00602     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00603 
00604     // Update parameter in our list and residual tests.
00605     params_->set("Convergence Tolerance", convtol_);
00606     if (convTest_ != Teuchos::null)
00607       convTest_->setTolerance( convtol_ );
00608   }
00609   
00610   if (params->isParameter("Show Maximum Residual Norm Only")) {
00611     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00612 
00613     // Update parameter in our list and residual tests
00614     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00615     if (convTest_ != Teuchos::null)
00616       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00617   }
00618 
00619   // Create status tests if we need to.
00620 
00621   // Basic test checks maximum iterations and native residual.
00622   if (maxIterTest_ == Teuchos::null)
00623     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00624   
00625   // Implicit residual test, using the native residual to determine if convergence was achieved.
00626   if (convTest_ == Teuchos::null)
00627     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00628   
00629   if (sTest_ == Teuchos::null)
00630     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00631   
00632   if (outputTest_ == Teuchos::null) {
00633 
00634     // Create the status test output class.
00635     // This class manages and formats the output from the status test.
00636     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00637     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00638 
00639     // Set the solver string for the output test
00640     std::string solverDesc = " Block CG ";
00641     outputTest_->setSolverDesc( solverDesc );
00642 
00643   }
00644 
00645   // Create orthogonalization manager if we need to.
00646   if (ortho_ == Teuchos::null) {
00647     if (orthoType_=="DGKS") {
00648       if (orthoKappa_ <= 0) {
00649   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00650       }
00651       else {
00652   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00653   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00654       }
00655     }
00656     else if (orthoType_=="ICGS") {
00657       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00658     } 
00659     else if (orthoType_=="IMGS") {
00660       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00661     } 
00662     else {
00663       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00664        "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
00665     }  
00666   }
00667 
00668   // Create the timer if we need to.
00669   if (timerSolve_ == Teuchos::null) {
00670     std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
00671 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00672     timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
00673 #endif
00674   }
00675 
00676   // Inform the solver manager that the current parameters were set.
00677   isSet_ = true;
00678 }
00679 
00680     
00681 template<class ScalarType, class MV, class OP>
00682 Teuchos::RCP<const Teuchos::ParameterList> 
00683 BlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00684 {
00685   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00686   
00687   // Set all the valid parameters and their default values.
00688   if(is_null(validPL)) {
00689     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00690     pl->set("Convergence Tolerance", convtol_default_,
00691       "The relative residual tolerance that needs to be achieved by the\n"
00692       "iterative solver in order for the linear system to be declared converged.");
00693     pl->set("Maximum Iterations", maxIters_default_,
00694       "The maximum number of block iterations allowed for each\n"
00695       "set of RHS solved.");
00696     pl->set("Block Size", blockSize_default_,
00697       "The number of vectors in each block.");
00698     pl->set("Adaptive Block Size", adaptiveBlockSize_default_,
00699       "Whether the solver manager should adapt to the block size\n"
00700       "based on the number of RHS to solve.");
00701     pl->set("Verbosity", verbosity_default_,
00702       "What type(s) of solver information should be outputted\n"
00703       "to the output stream.");
00704     pl->set("Output Style", outputStyle_default_,
00705       "What style is used for the solver information outputted\n"
00706       "to the output stream.");
00707     pl->set("Output Frequency", outputFreq_default_,
00708       "How often convergence information should be outputted\n"
00709       "to the output stream.");  
00710     pl->set("Output Stream", outputStream_default_,
00711       "A reference-counted pointer to the output stream where all\n"
00712       "solver output is sent.");
00713     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00714       "When convergence information is printed, only show the maximum\n"
00715       "relative residual norm when the block size is greater than one.");
00716     pl->set("Timer Label", label_default_,
00717       "The string to use as a prefix for the timer labels.");
00718     //  pl->set("Restart Timers", restartTimers_);
00719     pl->set("Orthogonalization", orthoType_default_,
00720       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00721     pl->set("Orthogonalization Constant",orthoKappa_default_,
00722       "The constant used by DGKS orthogonalization to determine\n"
00723       "whether another step of classical Gram-Schmidt is necessary.");
00724     validPL = pl;
00725   }
00726   return validPL;
00727 }
00728 
00729   
00730 // solve()
00731 template<class ScalarType, class MV, class OP>
00732 ReturnType BlockCGSolMgr<ScalarType,MV,OP>::solve() {
00733   using Teuchos::RCP;
00734   using Teuchos::rcp;
00735   using Teuchos::rcp_const_cast;
00736   using Teuchos::rcp_dynamic_cast;
00737 
00738   // Set the current parameters if they were not set before.  NOTE:
00739   // This may occur if the user generated the solver manager with the
00740   // default constructor and then didn't set any parameters using
00741   // setParameters().
00742   if (!isSet_) {
00743     setParameters(Teuchos::parameterList(*getValidParameters()));
00744   }
00745 
00746   Teuchos::BLAS<int,ScalarType> blas;
00747   Teuchos::LAPACK<int,ScalarType> lapack;
00748   
00749   TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(), 
00750     BlockCGSolMgrLinearProblemFailure,
00751     "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
00752     "has not been called.");
00753 
00754   // Create indices for the linear systems to be solved.
00755   int startPtr = 0;
00756   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00757   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00758 
00759   std::vector<int> currIdx, currIdx2;
00760   //  If an adaptive block size is allowed then only the linear
00761   //  systems that need to be solved are solved.  Otherwise, the index
00762   //  set is generated that informs the linear problem that some
00763   //  linear systems are augmented.
00764   if ( adaptiveBlockSize_ ) {
00765     blockSize_ = numCurrRHS;
00766     currIdx.resize( numCurrRHS  );
00767     currIdx2.resize( numCurrRHS  );
00768     for (int i=0; i<numCurrRHS; ++i) 
00769       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00770     
00771   }
00772   else {
00773     currIdx.resize( blockSize_ );
00774     currIdx2.resize( blockSize_ );
00775     for (int i=0; i<numCurrRHS; ++i) 
00776       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00777     for (int i=numCurrRHS; i<blockSize_; ++i) 
00778       { currIdx[i] = -1; currIdx2[i] = i; }
00779   }
00780 
00781   // Inform the linear problem of the current linear system to solve.
00782   problem_->setLSIndex( currIdx );
00783 
00785   // Set up the parameter list for the Iteration subclass.
00786   Teuchos::ParameterList plist;
00787   plist.set("Block Size",blockSize_);
00788   
00789   // Reset the output status test (controls all the other status tests).
00790   outputTest_->reset();
00791 
00792   // Assume convergence is achieved, then let any failed convergence
00793   // set this to false.  "Innocent until proven guilty."
00794   bool isConverged = true;  
00795 
00797   // Set up the BlockCG Iteration subclass.
00798 
00799   RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
00800   if (blockSize_ == 1) {
00801     // Standard (nonblock) CG is faster for the special case of a
00802     // block size of 1.
00803     block_cg_iter = 
00804       rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, 
00805            outputTest_, plist));
00806   } else {
00807     block_cg_iter = 
00808       rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, 
00809                 ortho_, plist));
00810   }
00811   
00812 
00813   // Enter solve() iterations
00814   {
00815 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00816     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00817 #endif
00818 
00819     while ( numRHS2Solve > 0 ) {
00820       //
00821       // Reset the active / converged vectors from this block
00822       std::vector<int> convRHSIdx;
00823       std::vector<int> currRHSIdx( currIdx );
00824       currRHSIdx.resize(numCurrRHS);
00825 
00826       // Reset the number of iterations.
00827       block_cg_iter->resetNumIters();
00828 
00829       // Reset the number of calls that the status test output knows about.
00830       outputTest_->resetNumCalls();
00831 
00832       // Get the current residual for this block of linear systems.
00833       RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00834 
00835       // Set the new state and initialize the solver.
00836       CGIterationState<ScalarType,MV> newstate;
00837       newstate.R = R_0;
00838       block_cg_iter->initializeCG(newstate);
00839 
00840       while(1) {
00841   
00842   // tell block_cg_iter to iterate
00843   try {
00844     block_cg_iter->iterate();
00845     //
00846     // Check whether any of the linear systems converged.
00847     //
00848     if (convTest_->getStatus() == Passed) {
00849       // At least one of the linear system(s) converged.  
00850       //
00851       // Get the column indices of the linear systems that converged.
00852       typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
00853       std::vector<int> convIdx = 
00854         rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
00855 
00856       // If the number of converged linear systems equals the
00857             // number of linear systems currently being solved, then
00858             // we are done with this block.
00859       if (convIdx.size() == currRHSIdx.size())
00860         break;  // break from while(1){block_cg_iter->iterate()}
00861 
00862       // Inform the linear problem that we are finished with
00863       // this current linear system.
00864       problem_->setCurrLS();
00865 
00866       // Reset currRHSIdx to contain the right-hand sides that
00867       // are left to converge for this block.
00868       int have = 0;
00869             std::vector<int> unconvIdx(currRHSIdx.size());
00870       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00871         bool found = false;
00872         for (unsigned int j=0; j<convIdx.size(); ++j) {
00873     if (currRHSIdx[i] == convIdx[j]) {
00874       found = true;
00875       break;
00876     }
00877         }
00878         if (!found) {
00879                 currIdx2[have] = currIdx2[i];
00880     currRHSIdx[have++] = currRHSIdx[i];
00881         }
00882               else {
00883               } 
00884       }
00885       currRHSIdx.resize(have);
00886       currIdx2.resize(have);
00887 
00888       // Set the remaining indices after deflation.
00889       problem_->setLSIndex( currRHSIdx );
00890 
00891       // Get the current residual vector.
00892       std::vector<MagnitudeType> norms;
00893             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00894       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00895 
00896       // Set the new blocksize for the solver.
00897       block_cg_iter->setBlockSize( have );
00898 
00899       // Set the new state and initialize the solver.
00900       CGIterationState<ScalarType,MV> defstate;
00901       defstate.R = R_0;
00902       block_cg_iter->initializeCG(defstate);
00903     }
00904     //
00905     // None of the linear systems converged.  Check whether the
00906     // maximum iteration count was reached.
00907     //
00908     else if (maxIterTest_->getStatus() == Passed) {
00909       isConverged = false; // None of the linear systems converged.
00910       break;  // break from while(1){block_cg_iter->iterate()}
00911     }
00912     //
00913     // iterate() returned, but none of our status tests Passed.
00914     // This indicates a bug.
00915     //
00916     else {
00917       TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00918         "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
00919         "the maximum iteration count test passed.  Please report this bug "
00920         "to the Belos developers.");
00921     }
00922   }
00923   catch (const std::exception &e) {
00924     std::ostream& err = printer_->stream (Errors);
00925     err << "Error! Caught std::exception in CGIteration::iterate() at "
00926         << "iteration " << block_cg_iter->getNumIters() << std::endl 
00927         << e.what() << std::endl;
00928     throw;
00929   }
00930       }
00931       
00932       // Inform the linear problem that we are finished with this
00933       // block linear system.
00934       problem_->setCurrLS();
00935       
00936       // Update indices for the linear systems to be solved.
00937       startPtr += numCurrRHS;
00938       numRHS2Solve -= numCurrRHS;
00939       if ( numRHS2Solve > 0 ) {
00940   numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00941 
00942   if ( adaptiveBlockSize_ ) {
00943           blockSize_ = numCurrRHS;
00944     currIdx.resize( numCurrRHS  );
00945     currIdx2.resize( numCurrRHS  );
00946     for (int i=0; i<numCurrRHS; ++i) 
00947       { currIdx[i] = startPtr+i; currIdx2[i] = i; }   
00948   }
00949   else {
00950     currIdx.resize( blockSize_ );
00951     currIdx2.resize( blockSize_ );
00952     for (int i=0; i<numCurrRHS; ++i) 
00953       { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00954     for (int i=numCurrRHS; i<blockSize_; ++i) 
00955       { currIdx[i] = -1; currIdx2[i] = i; }
00956   }
00957   // Set the next indices.
00958   problem_->setLSIndex( currIdx );
00959 
00960   // Set the new blocksize for the solver.
00961   block_cg_iter->setBlockSize( blockSize_ );  
00962       }
00963       else {
00964         currIdx.resize( numRHS2Solve );
00965       }
00966       
00967     }// while ( numRHS2Solve > 0 )
00968     
00969   }
00970 
00971   // print final summary
00972   sTest_->print( printer_->stream(FinalSummary) );
00973  
00974   // print timing information
00975 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00976   // Calling summarize() requires communication in general, so don't
00977   // call it unless the user wants to print out timing details.
00978   // summarize() will do all the work even if it's passed a "black
00979   // hole" output stream.
00980   if (verbosity_ & TimingDetails) {
00981     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00982   }
00983 #endif
00984  
00985   // Save the iteration count for this solve.
00986   numIters_ = maxIterTest_->getNumIters();
00987 
00988   // Save the convergence test value ("achieved tolerance") for this solve.
00989   {
00990     typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
00991     // testValues is nonnull and not persistent.
00992     const std::vector<MagnitudeType>* pTestValues = 
00993       rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
00994     
00995     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
00996       "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
00997       "method returned NULL.  Please report this bug to the Belos developers.");
00998     
00999     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
01000       "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
01001       "method returned a vector of length zero.  Please report this bug to the "
01002       "Belos developers.");
01003 
01004     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01005     // achieved tolerances for all vectors in the current solve(), or
01006     // just for the vectors from the last deflation?
01007     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01008   }
01009  
01010   if (!isConverged) {
01011     return Unconverged; // return from BlockCGSolMgr::solve() 
01012   }
01013   return Converged; // return from BlockCGSolMgr::solve() 
01014 }
01015 
01016 //  This method requires the solver manager to return a std::string that describes itself.
01017 template<class ScalarType, class MV, class OP>
01018 std::string BlockCGSolMgr<ScalarType,MV,OP>::description() const
01019 {
01020   std::ostringstream oss;
01021   oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01022   oss << "{";
01023   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
01024   oss << "}";
01025   return oss.str();
01026 }
01027   
01028 } // end Belos namespace
01029 
01030 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines