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::getNewTimer(solveLabel);
00495 #endif
00496     }
00497   }
00498 
00499   // Check if the orthogonalization changed.
00500   if (params->isParameter("Orthogonalization")) {
00501     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00502     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00503                         std::invalid_argument,
00504       "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00505     if (tempOrthoType != orthoType_) {
00506       orthoType_ = tempOrthoType;
00507       // Create orthogonalization manager
00508       if (orthoType_=="DGKS") {
00509   if (orthoKappa_ <= 0) {
00510     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00511   }
00512   else {
00513     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00514     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00515   }
00516       }
00517       else if (orthoType_=="ICGS") {
00518   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00519       } 
00520       else if (orthoType_=="IMGS") {
00521   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00522       } 
00523     }  
00524   }
00525 
00526   // Check which orthogonalization constant to use.
00527   if (params->isParameter("Orthogonalization Constant")) {
00528     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00529 
00530     // Update parameter in our list.
00531     params_->set("Orthogonalization Constant",orthoKappa_);
00532     if (orthoType_=="DGKS") {
00533       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00534   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00535       }
00536     } 
00537   }
00538 
00539   // Check for a change in verbosity level
00540   if (params->isParameter("Verbosity")) {
00541     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00542       verbosity_ = params->get("Verbosity", verbosity_default_);
00543     } else {
00544       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00545     }
00546 
00547     // Update parameter in our list.
00548     params_->set("Verbosity", verbosity_);
00549     if (printer_ != Teuchos::null)
00550       printer_->setVerbosity(verbosity_);
00551   }
00552 
00553   // Check for a change in output style
00554   if (params->isParameter("Output Style")) {
00555     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00556       outputStyle_ = params->get("Output Style", outputStyle_default_);
00557     } else {
00558       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00559     }
00560 
00561     // Update parameter in our list.
00562     params_->set("Output Style", outputStyle_);
00563     outputTest_ = Teuchos::null;
00564   }
00565 
00566   // output stream
00567   if (params->isParameter("Output Stream")) {
00568     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00569 
00570     // Update parameter in our list.
00571     params_->set("Output Stream", outputStream_);
00572     if (printer_ != Teuchos::null)
00573       printer_->setOStream( outputStream_ );
00574   }
00575 
00576   // frequency level
00577   if (verbosity_ & Belos::StatusTestDetails) {
00578     if (params->isParameter("Output Frequency")) {
00579       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00580     }
00581 
00582     // Update parameter in out list and output status test.
00583     params_->set("Output Frequency", outputFreq_);
00584     if (outputTest_ != Teuchos::null)
00585       outputTest_->setOutputFrequency( outputFreq_ );
00586   }
00587 
00588   // Create output manager if we need to.
00589   if (printer_ == Teuchos::null) {
00590     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00591   }  
00592   
00593   // Convergence
00594   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00595   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00596 
00597   // Check for convergence tolerance
00598   if (params->isParameter("Convergence Tolerance")) {
00599     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00600 
00601     // Update parameter in our list and residual tests.
00602     params_->set("Convergence Tolerance", convtol_);
00603     if (convTest_ != Teuchos::null)
00604       convTest_->setTolerance( convtol_ );
00605   }
00606   
00607   if (params->isParameter("Show Maximum Residual Norm Only")) {
00608     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00609 
00610     // Update parameter in our list and residual tests
00611     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00612     if (convTest_ != Teuchos::null)
00613       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00614   }
00615 
00616   // Create status tests if we need to.
00617 
00618   // Basic test checks maximum iterations and native residual.
00619   if (maxIterTest_ == Teuchos::null)
00620     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00621   
00622   // Implicit residual test, using the native residual to determine if convergence was achieved.
00623   if (convTest_ == Teuchos::null)
00624     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00625   
00626   if (sTest_ == Teuchos::null)
00627     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00628   
00629   if (outputTest_ == Teuchos::null) {
00630 
00631     // Create the status test output class.
00632     // This class manages and formats the output from the status test.
00633     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00634     outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00635 
00636     // Set the solver string for the output test
00637     std::string solverDesc = " Block CG ";
00638     outputTest_->setSolverDesc( solverDesc );
00639 
00640   }
00641 
00642   // Create orthogonalization manager if we need to.
00643   if (ortho_ == Teuchos::null) {
00644     if (orthoType_=="DGKS") {
00645       if (orthoKappa_ <= 0) {
00646   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00647       }
00648       else {
00649   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00650   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00651       }
00652     }
00653     else if (orthoType_=="ICGS") {
00654       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00655     } 
00656     else if (orthoType_=="IMGS") {
00657       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00658     } 
00659     else {
00660       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00661        "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
00662     }  
00663   }
00664 
00665   // Create the timer if we need to.
00666   if (timerSolve_ == Teuchos::null) {
00667     std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
00668 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00669     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00670 #endif
00671   }
00672 
00673   // Inform the solver manager that the current parameters were set.
00674   isSet_ = true;
00675 }
00676 
00677     
00678 template<class ScalarType, class MV, class OP>
00679 Teuchos::RCP<const Teuchos::ParameterList> 
00680 BlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00681 {
00682   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00683   
00684   // Set all the valid parameters and their default values.
00685   if(is_null(validPL)) {
00686     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00687     pl->set("Convergence Tolerance", convtol_default_,
00688       "The relative residual tolerance that needs to be achieved by the\n"
00689       "iterative solver in order for the linear system to be declared converged.");
00690     pl->set("Maximum Iterations", maxIters_default_,
00691       "The maximum number of block iterations allowed for each\n"
00692       "set of RHS solved.");
00693     pl->set("Block Size", blockSize_default_,
00694       "The number of vectors in each block.");
00695     pl->set("Adaptive Block Size", adaptiveBlockSize_default_,
00696       "Whether the solver manager should adapt to the block size\n"
00697       "based on the number of RHS to solve.");
00698     pl->set("Verbosity", verbosity_default_,
00699       "What type(s) of solver information should be outputted\n"
00700       "to the output stream.");
00701     pl->set("Output Style", outputStyle_default_,
00702       "What style is used for the solver information outputted\n"
00703       "to the output stream.");
00704     pl->set("Output Frequency", outputFreq_default_,
00705       "How often convergence information should be outputted\n"
00706       "to the output stream.");  
00707     pl->set("Output Stream", outputStream_default_,
00708       "A reference-counted pointer to the output stream where all\n"
00709       "solver output is sent.");
00710     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00711       "When convergence information is printed, only show the maximum\n"
00712       "relative residual norm when the block size is greater than one.");
00713     pl->set("Timer Label", label_default_,
00714       "The string to use as a prefix for the timer labels.");
00715     //  pl->set("Restart Timers", restartTimers_);
00716     pl->set("Orthogonalization", orthoType_default_,
00717       "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
00718     pl->set("Orthogonalization Constant",orthoKappa_default_,
00719       "The constant used by DGKS orthogonalization to determine\n"
00720       "whether another step of classical Gram-Schmidt is necessary.");
00721     validPL = pl;
00722   }
00723   return validPL;
00724 }
00725 
00726   
00727 // solve()
00728 template<class ScalarType, class MV, class OP>
00729 ReturnType BlockCGSolMgr<ScalarType,MV,OP>::solve() {
00730   using Teuchos::RCP;
00731   using Teuchos::rcp;
00732   using Teuchos::rcp_const_cast;
00733   using Teuchos::rcp_dynamic_cast;
00734 
00735   // Set the current parameters if they were not set before.  NOTE:
00736   // This may occur if the user generated the solver manager with the
00737   // default constructor and then didn't set any parameters using
00738   // setParameters().
00739   if (!isSet_) {
00740     setParameters(Teuchos::parameterList(*getValidParameters()));
00741   }
00742 
00743   Teuchos::BLAS<int,ScalarType> blas;
00744   Teuchos::LAPACK<int,ScalarType> lapack;
00745   
00746   TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(), 
00747     BlockCGSolMgrLinearProblemFailure,
00748     "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
00749     "has not been called.");
00750 
00751   // Create indices for the linear systems to be solved.
00752   int startPtr = 0;
00753   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00754   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00755 
00756   std::vector<int> currIdx, currIdx2;
00757   //  If an adaptive block size is allowed then only the linear
00758   //  systems that need to be solved are solved.  Otherwise, the index
00759   //  set is generated that informs the linear problem that some
00760   //  linear systems are augmented.
00761   if ( adaptiveBlockSize_ ) {
00762     blockSize_ = numCurrRHS;
00763     currIdx.resize( numCurrRHS  );
00764     currIdx2.resize( numCurrRHS  );
00765     for (int i=0; i<numCurrRHS; ++i) 
00766       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00767     
00768   }
00769   else {
00770     currIdx.resize( blockSize_ );
00771     currIdx2.resize( blockSize_ );
00772     for (int i=0; i<numCurrRHS; ++i) 
00773       { currIdx[i] = startPtr+i; currIdx2[i]=i; }
00774     for (int i=numCurrRHS; i<blockSize_; ++i) 
00775       { currIdx[i] = -1; currIdx2[i] = i; }
00776   }
00777 
00778   // Inform the linear problem of the current linear system to solve.
00779   problem_->setLSIndex( currIdx );
00780 
00782   // Set up the parameter list for the Iteration subclass.
00783   Teuchos::ParameterList plist;
00784   plist.set("Block Size",blockSize_);
00785   
00786   // Reset the output status test (controls all the other status tests).
00787   outputTest_->reset();
00788 
00789   // Assume convergence is achieved, then let any failed convergence
00790   // set this to false.  "Innocent until proven guilty."
00791   bool isConverged = true;  
00792 
00794   // Set up the BlockCG Iteration subclass.
00795 
00796   RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
00797   if (blockSize_ == 1) {
00798     // Standard (nonblock) CG is faster for the special case of a
00799     // block size of 1.
00800     block_cg_iter = 
00801       rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, 
00802            outputTest_, plist));
00803   } else {
00804     block_cg_iter = 
00805       rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, 
00806                 ortho_, plist));
00807   }
00808   
00809 
00810   // Enter solve() iterations
00811   {
00812 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00813     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00814 #endif
00815 
00816     while ( numRHS2Solve > 0 ) {
00817       //
00818       // Reset the active / converged vectors from this block
00819       std::vector<int> convRHSIdx;
00820       std::vector<int> currRHSIdx( currIdx );
00821       currRHSIdx.resize(numCurrRHS);
00822 
00823       // Reset the number of iterations.
00824       block_cg_iter->resetNumIters();
00825 
00826       // Reset the number of calls that the status test output knows about.
00827       outputTest_->resetNumCalls();
00828 
00829       // Get the current residual for this block of linear systems.
00830       RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00831 
00832       // Set the new state and initialize the solver.
00833       CGIterationState<ScalarType,MV> newstate;
00834       newstate.R = R_0;
00835       block_cg_iter->initializeCG(newstate);
00836 
00837       while(1) {
00838   
00839   // tell block_cg_iter to iterate
00840   try {
00841     block_cg_iter->iterate();
00842     //
00843     // Check whether any of the linear systems converged.
00844     //
00845     if (convTest_->getStatus() == Passed) {
00846       // At least one of the linear system(s) converged.  
00847       //
00848       // Get the column indices of the linear systems that converged.
00849       typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
00850       std::vector<int> convIdx = 
00851         rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
00852 
00853       // If the number of converged linear systems equals the
00854             // number of linear systems currently being solved, then
00855             // we are done with this block.
00856       if (convIdx.size() == currRHSIdx.size())
00857         break;  // break from while(1){block_cg_iter->iterate()}
00858 
00859       // Inform the linear problem that we are finished with
00860       // this current linear system.
00861       problem_->setCurrLS();
00862 
00863       // Reset currRHSIdx to contain the right-hand sides that
00864       // are left to converge for this block.
00865       int have = 0;
00866             std::vector<int> unconvIdx(currRHSIdx.size());
00867       for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00868         bool found = false;
00869         for (unsigned int j=0; j<convIdx.size(); ++j) {
00870     if (currRHSIdx[i] == convIdx[j]) {
00871       found = true;
00872       break;
00873     }
00874         }
00875         if (!found) {
00876                 currIdx2[have] = currIdx2[i];
00877     currRHSIdx[have++] = currRHSIdx[i];
00878         }
00879               else {
00880               } 
00881       }
00882       currRHSIdx.resize(have);
00883       currIdx2.resize(have);
00884 
00885       // Set the remaining indices after deflation.
00886       problem_->setLSIndex( currRHSIdx );
00887 
00888       // Get the current residual vector.
00889       std::vector<MagnitudeType> norms;
00890             R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00891       for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00892 
00893       // Set the new blocksize for the solver.
00894       block_cg_iter->setBlockSize( have );
00895 
00896       // Set the new state and initialize the solver.
00897       CGIterationState<ScalarType,MV> defstate;
00898       defstate.R = R_0;
00899       block_cg_iter->initializeCG(defstate);
00900     }
00901     //
00902     // None of the linear systems converged.  Check whether the
00903     // maximum iteration count was reached.
00904     //
00905     else if (maxIterTest_->getStatus() == Passed) {
00906       isConverged = false; // None of the linear systems converged.
00907       break;  // break from while(1){block_cg_iter->iterate()}
00908     }
00909     //
00910     // iterate() returned, but none of our status tests Passed.
00911     // This indicates a bug.
00912     //
00913     else {
00914       TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00915         "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
00916         "the maximum iteration count test passed.  Please report this bug "
00917         "to the Belos developers.");
00918     }
00919   }
00920   catch (const std::exception &e) {
00921     std::ostream& err = printer_->stream (Errors);
00922     err << "Error! Caught std::exception in CGIteration::iterate() at "
00923         << "iteration " << block_cg_iter->getNumIters() << std::endl 
00924         << e.what() << std::endl;
00925     throw;
00926   }
00927       }
00928       
00929       // Inform the linear problem that we are finished with this
00930       // block linear system.
00931       problem_->setCurrLS();
00932       
00933       // Update indices for the linear systems to be solved.
00934       startPtr += numCurrRHS;
00935       numRHS2Solve -= numCurrRHS;
00936       if ( numRHS2Solve > 0 ) {
00937   numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
00938 
00939   if ( adaptiveBlockSize_ ) {
00940           blockSize_ = numCurrRHS;
00941     currIdx.resize( numCurrRHS  );
00942     currIdx2.resize( numCurrRHS  );
00943     for (int i=0; i<numCurrRHS; ++i) 
00944       { currIdx[i] = startPtr+i; currIdx2[i] = i; }   
00945   }
00946   else {
00947     currIdx.resize( blockSize_ );
00948     currIdx2.resize( blockSize_ );
00949     for (int i=0; i<numCurrRHS; ++i) 
00950       { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00951     for (int i=numCurrRHS; i<blockSize_; ++i) 
00952       { currIdx[i] = -1; currIdx2[i] = i; }
00953   }
00954   // Set the next indices.
00955   problem_->setLSIndex( currIdx );
00956 
00957   // Set the new blocksize for the solver.
00958   block_cg_iter->setBlockSize( blockSize_ );  
00959       }
00960       else {
00961         currIdx.resize( numRHS2Solve );
00962       }
00963       
00964     }// while ( numRHS2Solve > 0 )
00965     
00966   }
00967 
00968   // print final summary
00969   sTest_->print( printer_->stream(FinalSummary) );
00970  
00971   // print timing information
00972 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00973   // Calling summarize() requires communication in general, so don't
00974   // call it unless the user wants to print out timing details.
00975   // summarize() will do all the work even if it's passed a "black
00976   // hole" output stream.
00977   if (verbosity_ & TimingDetails) {
00978     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00979   }
00980 #endif
00981  
00982   // Save the iteration count for this solve.
00983   numIters_ = maxIterTest_->getNumIters();
00984 
00985   // Save the convergence test value ("achieved tolerance") for this solve.
00986   {
00987     typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
00988     // testValues is nonnull and not persistent.
00989     const std::vector<MagnitudeType>* pTestValues = 
00990       rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
00991     
00992     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
00993       "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
00994       "method returned NULL.  Please report this bug to the Belos developers.");
00995     
00996     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
00997       "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
00998       "method returned a vector of length zero.  Please report this bug to the "
00999       "Belos developers.");
01000 
01001     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01002     // achieved tolerances for all vectors in the current solve(), or
01003     // just for the vectors from the last deflation?
01004     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01005   }
01006  
01007   if (!isConverged) {
01008     return Unconverged; // return from BlockCGSolMgr::solve() 
01009   }
01010   return Converged; // return from BlockCGSolMgr::solve() 
01011 }
01012 
01013 //  This method requires the solver manager to return a std::string that describes itself.
01014 template<class ScalarType, class MV, class OP>
01015 std::string BlockCGSolMgr<ScalarType,MV,OP>::description() const
01016 {
01017   std::ostringstream oss;
01018   oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01019   oss << "{";
01020   oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
01021   oss << "}";
01022   return oss.str();
01023 }
01024   
01025 } // end Belos namespace
01026 
01027 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines