BelosPCPGSolMgr.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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_PCPG_SOLMGR_HPP
00030 #define BELOS_PCPG_SOLMGR_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041 
00042 #include "BelosPCPGIter.hpp"
00043 
00044 #include "BelosDGKSOrthoManager.hpp"
00045 #include "BelosICGSOrthoManager.hpp"
00046 #include "BelosIMGSOrthoManager.hpp"
00047 #include "BelosStatusTestMaxIters.hpp"
00048 #include "BelosStatusTestGenResNorm.hpp"
00049 #include "BelosStatusTestCombo.hpp"
00050 #include "BelosStatusTestOutputFactory.hpp"
00051 #include "BelosOutputManager.hpp"
00052 #include "Teuchos_BLAS.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055 
00056 //#include <vector>              //getPermThatSorts
00057 //#include <algorithm>           //getPermThatSorts
00058 
00059 
00076 namespace Belos {
00077   
00079 
00080   
00087   class PCPGSolMgrLinearProblemFailure : public BelosError {public:
00088     PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00089     {}};
00090   
00096   class PCPGSolMgrOrthoFailure : public BelosError {public:
00097     PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00098     {}};
00099   
00106   class PCPGSolMgrLAPACKFailure : public BelosError {public:
00107     PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00108     {}};
00109 
00116   class PCPGSolMgrRecyclingFailure : public BelosError {public:
00117     PCPGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
00118     {}};
00119 
00121 
00122 
00123   template<class ScalarType, class MV, class OP>
00124   class PCPGSolMgr : public SolverManager<ScalarType,MV,OP> {
00125     
00126   private:
00127     typedef MultiVecTraits<ScalarType,MV> MVT;
00128     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00129     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00130     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00131     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00132     
00133   public:
00134     
00136 
00137    
00144      PCPGSolMgr();
00145  
00181     PCPGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00182           const Teuchos::RCP<Teuchos::ParameterList> &pl );
00183     
00185     virtual ~PCPGSolMgr() {};
00187     
00189 
00190     
00193     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00194       return *problem_;
00195     }
00196 
00199     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00200 
00203     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00204  
00210     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00211       return tuple(timerSolve_);
00212     }
00213 
00215     int getNumIters() const {
00216       return numIters_;
00217     }
00218  
00221     bool isLOADetected() const { return false; }
00222  
00224     
00226 
00227    
00229     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00230    
00232     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00233     
00235    
00237 
00238 
00242     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00244  
00246 
00247     
00265     ReturnType solve();
00266     
00268     
00271     
00273     std::string description() const;
00274     
00276     
00277   private:
00278 
00279     // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary.  Given 
00280     // the seed space U and C = A U, find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU.
00281     int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
00282 
00283     // Linear problem.
00284     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00285     
00286     // Output manager.
00287     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00288     Teuchos::RCP<std::ostream> outputStream_;
00289 
00290     // Status test.
00291     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00292     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00293     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00294     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00295 
00296     // Orthogonalization manager.
00297     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00298     
00299     // Current parameter list.
00300     Teuchos::RCP<ParameterList> params_;
00301 
00302     // Default solver values.
00303     static const MagnitudeType convtol_default_;
00304     static const MagnitudeType orthoKappa_default_;
00305     static const int maxIters_default_;
00306     static const int deflatedBlocks_default_;
00307     static const int savedBlocks_default_;
00308     static const int verbosity_default_;
00309     static const int outputStyle_default_;
00310     static const int outputFreq_default_;
00311     static const std::string label_default_;
00312     static const std::string orthoType_default_;
00313     static const Teuchos::RCP<std::ostream> outputStream_default_;
00314 
00315     // Current solver values.
00316     MagnitudeType convtol_, orthoKappa_;
00317     int numIters_, maxIters_, deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
00318     std::string orthoType_; 
00319 
00320     // Recycled subspace, its image and the residual
00321     Teuchos::RCP<MV> U_, C_, R_;
00322 
00323     // Actual dimension of current recycling subspace (<= savedBlocks_ )
00324     int dimU_; 
00325 
00326     // Timers.
00327     std::string label_;
00328     Teuchos::RCP<Teuchos::Time> timerSolve_;
00329 
00330     // Internal state variables.
00331     bool isSet_;
00332   };
00333 
00334 
00335 // Default solver values.
00336 template<class ScalarType, class MV, class OP>
00337 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00338 
00339 template<class ScalarType, class MV, class OP>
00340 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00341 
00342 template<class ScalarType, class MV, class OP>
00343 const int PCPGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00344 
00345 template<class ScalarType, class MV, class OP>
00346 const int PCPGSolMgr<ScalarType,MV,OP>::deflatedBlocks_default_ = 2;
00347 
00348 template<class ScalarType, class MV, class OP>
00349 const int PCPGSolMgr<ScalarType,MV,OP>::savedBlocks_default_ = 16;
00350 
00351 template<class ScalarType, class MV, class OP>
00352 const int PCPGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00353 
00354 template<class ScalarType, class MV, class OP>
00355 const int PCPGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00356 
00357 template<class ScalarType, class MV, class OP>
00358 const int PCPGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00359 
00360 template<class ScalarType, class MV, class OP>
00361 const std::string PCPGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00362 
00363 template<class ScalarType, class MV, class OP>
00364 const std::string PCPGSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00365 
00366 template<class ScalarType, class MV, class OP>
00367 const Teuchos::RCP<std::ostream> PCPGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00368 
00369 
00370 // Empty Constructor
00371 template<class ScalarType, class MV, class OP>
00372 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr() :
00373   outputStream_(outputStream_default_),
00374   convtol_(convtol_default_),
00375   orthoKappa_(orthoKappa_default_),
00376   maxIters_(maxIters_default_),
00377   deflatedBlocks_(deflatedBlocks_default_),
00378   savedBlocks_(savedBlocks_default_),
00379   verbosity_(verbosity_default_),
00380   outputStyle_(outputStyle_default_),
00381   outputFreq_(outputFreq_default_),
00382   orthoType_(orthoType_default_),
00383   label_(label_default_),
00384   isSet_(false)
00385 {}
00386 
00387 
00388 // Basic Constructor
00389 template<class ScalarType, class MV, class OP>
00390 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr( 
00391                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00392                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00393   problem_(problem),
00394   outputStream_(outputStream_default_),
00395   convtol_(convtol_default_),
00396   orthoKappa_(orthoKappa_default_),
00397   maxIters_(maxIters_default_),
00398   deflatedBlocks_(deflatedBlocks_default_),
00399   savedBlocks_(savedBlocks_default_),
00400   verbosity_(verbosity_default_),
00401   outputStyle_(outputStyle_default_),
00402   outputFreq_(outputFreq_default_),
00403   orthoType_(orthoType_default_),
00404   label_(label_default_),
00405   isSet_(false)
00406 {
00407   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00408 
00409   if (!is_null(pl)) {
00410     // Set the parameters using the list that was passed in.
00411     setParameters( pl );  
00412   }
00413 }
00414 
00415 
00416 template<class ScalarType, class MV, class OP>
00417 void PCPGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00418 {
00419   // Create the internal parameter list if ones doesn't already exist.
00420   if (params_ == Teuchos::null) {
00421     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00422   }
00423   else {
00424     params->validateParameters(*getValidParameters());
00425   }
00426 
00427   // Check for maximum number of iterations
00428   if (params->isParameter("Maximum Iterations")) {
00429     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00430 
00431     // Update parameter in our list and in status test.
00432     params_->set("Maximum Iterations", maxIters_);
00433     if (maxIterTest_!=Teuchos::null)
00434       maxIterTest_->setMaxIters( maxIters_ );
00435   }
00436 
00437   // Check for the maximum numbers of saved and deflated blocks.
00438   if (params->isParameter("Num Saved Blocks")) {
00439     savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
00440     TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
00441            "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
00442 
00443     // savedBlocks > number of matrix rows and columns, not known in parameters.
00444     //TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
00445     //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
00446 
00447     // Update parameter in our list.
00448     params_->set("Num Saved Blocks", savedBlocks_);
00449   }
00450   if (params->isParameter("Num Deflated Blocks")) {
00451     deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
00452     TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
00453            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
00454 
00455     TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
00456            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
00457 
00458     // Update parameter in our list.
00459     params_->set("Num Deflated Blocks", deflatedBlocks_);
00460   }
00461 
00462   // Check to see if the timer label changed.
00463   if (params->isParameter("Timer Label")) {
00464     std::string tempLabel = params->get("Timer Label", label_default_);
00465 
00466     // Update parameter in our list and solver timer
00467     if (tempLabel != label_) {
00468       label_ = tempLabel;
00469       params_->set("Timer Label", label_);
00470       std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00471       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00472     }
00473   }
00474 
00475   // Check if the orthogonalization changed.
00476   if (params->isParameter("Orthogonalization")) {
00477     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00478     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00479       std::invalid_argument,
00480       "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00481     if (tempOrthoType != orthoType_) {
00482       orthoType_ = tempOrthoType;
00483       // Create orthogonalization manager
00484       if (orthoType_=="DGKS") {
00485   if (orthoKappa_ <= 0) {
00486     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00487   }
00488   else {
00489     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00490     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00491   }
00492       }
00493       else if (orthoType_=="ICGS") {
00494   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00495       } 
00496       else if (orthoType_=="IMGS") {
00497   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00498       } 
00499     }  
00500   }
00501 
00502   // Check which orthogonalization constant to use.
00503   if (params->isParameter("Orthogonalization Constant")) {
00504     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00505 
00506     // Update parameter in our list.
00507     params_->set("Orthogonalization Constant",orthoKappa_);
00508     if (orthoType_=="DGKS") {
00509       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00510   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00511       }
00512     } 
00513   }
00514 
00515   // Check for a change in verbosity level
00516   if (params->isParameter("Verbosity")) {
00517     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00518       verbosity_ = params->get("Verbosity", verbosity_default_);
00519     } else {
00520       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00521     }
00522 
00523     // Update parameter in our list.
00524     params_->set("Verbosity", verbosity_);
00525     if (printer_ != Teuchos::null)
00526       printer_->setVerbosity(verbosity_);
00527   }
00528 
00529   // Check for a change in output style
00530   if (params->isParameter("Output Style")) {
00531     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00532       outputStyle_ = params->get("Output Style", outputStyle_default_);
00533     } else {
00534       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00535     }
00536 
00537     // Reconstruct the convergence test if the explicit residual test is not being used.
00538     params_->set("Output Style", outputStyle_);
00539     outputTest_ = Teuchos::null;
00540   }
00541 
00542   // output stream
00543   if (params->isParameter("Output Stream")) {
00544     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00545 
00546     // Update parameter in our list.
00547     params_->set("Output Stream", outputStream_);
00548     if (printer_ != Teuchos::null)
00549       printer_->setOStream( outputStream_ );
00550   }
00551 
00552   // frequency level
00553   if (verbosity_ & Belos::StatusTestDetails) {
00554     if (params->isParameter("Output Frequency")) {
00555       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00556     }
00557 
00558     // Update parameter in out list and output status test.
00559     params_->set("Output Frequency", outputFreq_);
00560     if (outputTest_ != Teuchos::null)
00561       outputTest_->setOutputFrequency( outputFreq_ );
00562   }
00563 
00564   // Create output manager if we need to.
00565   if (printer_ == Teuchos::null) {
00566     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00567   }  
00568   
00569   // Convergence
00570   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00571   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00572 
00573   // Check for convergence tolerance
00574   if (params->isParameter("Convergence Tolerance")) {
00575     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00576 
00577     // Update parameter in our list and residual tests.
00578     params_->set("Convergence Tolerance", convtol_);
00579     if (convTest_ != Teuchos::null)
00580       convTest_->setTolerance( convtol_ );
00581   }
00582 
00583   // Create status tests if we need to.
00584 
00585   // Basic test checks maximum iterations and native residual.
00586   if (maxIterTest_ == Teuchos::null)
00587     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00588 
00589   if (convTest_ == Teuchos::null)
00590     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00591 
00592   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00593   
00594   // Create the status test output class.
00595   // This class manages and formats the output from the status test.
00596   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00597   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00598 
00599   // Set the solver string for the output test
00600   std::string solverDesc = " PCPG ";
00601   outputTest_->setSolverDesc( solverDesc );
00602 
00603 
00604   // Create orthogonalization manager if we need to.
00605   if (ortho_ == Teuchos::null) {
00606     if (orthoType_=="DGKS") {
00607       if (orthoKappa_ <= 0) {
00608   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00609       }
00610       else {
00611   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00612   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00613       }
00614     }
00615     else if (orthoType_=="ICGS") {
00616       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00617     } 
00618     else if (orthoType_=="IMGS") {
00619       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00620     } 
00621     else {
00622       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00623        "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
00624     }  
00625   }
00626 
00627   // Create the timer if we need to.
00628   if (timerSolve_ == Teuchos::null) {
00629     std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00630     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00631   }
00632 
00633   // Inform the solver manager that the current parameters were set.
00634   isSet_ = true;
00635 }
00636 
00637     
00638 template<class ScalarType, class MV, class OP>
00639 Teuchos::RCP<const Teuchos::ParameterList>
00640 PCPGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00641 {
00642   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00643   if (is_null(validPL)) {
00644     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00645     // Set all the valid parameters and their default values.
00646     pl->set("Convergence Tolerance", convtol_default_,
00647       "The relative residual tolerance that needs to be achieved by the\n"
00648       "iterative solver in order for the linear system to be declared converged.");
00649     pl->set("Maximum Iterations", maxIters_default_,
00650       "The maximum number of iterations allowed for each\n"
00651       "set of RHS solved.");
00652     pl->set("Num Deflated Blocks", deflatedBlocks_default_,
00653       "The maximum number of vectors in the seed subspace." );
00654     pl->set("Num Saved Blocks", savedBlocks_default_,
00655       "The maximum number of vectors saved from old Krylov subspaces." );
00656     pl->set("Verbosity", verbosity_default_,
00657       "What type(s) of solver information should be outputted\n"
00658       "to the output stream.");
00659     pl->set("Output Style", outputStyle_default_,
00660       "What style is used for the solver information outputted\n"
00661       "to the output stream.");
00662     pl->set("Output Frequency", outputFreq_default_,
00663       "How often convergence information should be outputted\n"
00664       "to the output stream.");  
00665     pl->set("Output Stream", outputStream_default_,
00666       "A reference-counted pointer to the output stream where all\n"
00667       "solver output is sent.");
00668     pl->set("Timer Label", label_default_,
00669       "The string to use as a prefix for the timer labels.");
00670     //  pl->set("Restart Timers", restartTimers_);
00671     pl->set("Orthogonalization", orthoType_default_,
00672       "The type of orthogonalization to use: DGKS, ICGS, IMGS");
00673     pl->set("Orthogonalization Constant",orthoKappa_default_,
00674       "The constant used by DGKS orthogonalization to determine\n"
00675       "whether another step of classical Gram-Schmidt is necessary.");
00676     validPL = pl;
00677   }
00678   return validPL;
00679 }
00680 
00681   
00682 // solve()
00683 template<class ScalarType, class MV, class OP>
00684 ReturnType PCPGSolMgr<ScalarType,MV,OP>::solve() {
00685 
00686   // Set the current parameters if are not set already.
00687   if (!isSet_) { setParameters( params_ ); }
00688 
00689   Teuchos::BLAS<int,ScalarType> blas;
00690   Teuchos::LAPACK<int,ScalarType> lapack;
00691   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00692   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00693   
00694   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PCPGSolMgrLinearProblemFailure,
00695                      "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
00696 
00697   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PCPGSolMgrLinearProblemFailure,
00698                      "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00699 
00700   // Create indices for the linear systems to be solved.
00701   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00702   std::vector<int> currIdx(1);
00703   currIdx[0] = 0;
00704 
00705    bool debug = false;
00706 
00707   // Inform the linear problem of the current linear system to solve.
00708   problem_->setLSIndex( currIdx ); // block size == 1
00709 
00710   // Assume convergence is achieved, then let any failed convergence set this to false.
00711   bool isConverged = true;  
00712 
00714   // PCPG iteration parameter list
00715   Teuchos::ParameterList plist;
00716   plist.set("Saved Blocks", savedBlocks_);
00717   plist.set("Block Size", 1);
00718   plist.set("Keep Diagonal", true);
00719   plist.set("Initialize Diagonal", true);
00720   
00722   // PCPG solver
00723   
00724   Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
00725   pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00726   // Number of iterations required to generate initial recycle space (if needed)
00727 
00728   // Enter solve() iterations
00729   {
00730     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00731     while ( numRHS2Solve > 0 ) {  // test for quick return
00732       
00733       // Reset the status test.
00734       outputTest_->reset();
00735 
00736       // Create the first block in the current Krylov basis (residual).
00737       if (R_ == Teuchos::null)
00738         R_ = MVT::Clone( *(problem_->getRHS()), 1 ); 
00739 
00740       problem_->computeCurrResVec( &*R_ );
00741 
00742       
00743       // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00744       // TODO: ensure hypothesis right here ... I have to think about use cases.
00745 
00746       if( U_ != Teuchos::null ){
00747         // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00748 
00749         // possibly over solved equation ...  I want residual norms
00750         // relative to the initial residual, not what I am about to compute.
00751         Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
00752         std::vector<MagnitudeType> rnorm0(1);
00753         MVT::MvNorm( *R_, rnorm0 ); // rnorm0  = norm(R_);
00754 
00755         // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
00756         std::cout  << "Solver Manager:  dimU_ = " << dimU_ << std::endl;
00757         Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
00758 
00759         Teuchos::RCP<const MV> Uactive, Cactive;
00760         std::vector<int> active_columns( dimU_ );
00761         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
00762         Uactive = MVT::CloneView(*U_, active_columns);
00763         Cactive = MVT::CloneView(*C_, active_columns);
00764 
00765         if( debug ){
00766           std::cout << " Solver Manager : check duality of seed basis " << std::endl;
00767           Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
00768           MVT::MvTransMv( one, *Uactive, *Cactive, H );
00769           H.print( std::cout );
00770         }
00771         
00772         MVT::MvTransMv( one, *Uactive, *R_, Z );
00773         Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
00774         MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );  // UZ
00775         MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += tmp;
00776         MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );  // CZ
00777         MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= tmp;
00778         std::vector<MagnitudeType> rnorm(1);
00779         MVT::MvNorm( *R_, rnorm );
00780         if( rnorm[0] < rnorm0[0] * .001 ){  //reorthogonalize
00781           MVT::MvTransMv( one, *Uactive, *R_, Z );
00782           MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
00783           MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += UZ;
00784           MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
00785           MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= CZ;
00786         }
00787         Uactive = Teuchos::null;
00788         Cactive = Teuchos::null;
00789         tempU = Teuchos::null;
00790       }
00791       else { 
00792         dimU_ = 0;
00793       }
00794 
00795 
00796       // Set the new state and initialize the solver.
00797       PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
00798 
00799       pcpgState.R = R_;
00800       if( U_ != Teuchos::null ) pcpgState.U = U_;
00801       if( C_ != Teuchos::null ) pcpgState.C = C_;
00802       if( dimU_ > 0 ) pcpgState.curDim = dimU_;
00803       pcpg_iter->initialize(pcpgState);
00804 
00805       // treat initialize() exceptions here?  how to use try-catch-throw? DMD
00806 
00807       // Get the current number of deflated blocks with the PCPG iteration
00808       dimU_ = pcpgState.curDim;
00809       if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
00810       pcpg_iter->resetNumIters();
00811 
00812       if( dimU_ > savedBlocks_ )
00813         std::cout << "Error: dimU_  = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl; 
00814 
00815       while(1) { // dummy loop for break
00816 
00817         // tell pcpg_iter to iterate
00818         try {
00819           if( debug ) printf("********** Calling iterate...\n");
00820           pcpg_iter->iterate();
00821 
00823           //
00824           // check convergence first
00825           //
00827           if ( convTest_->getStatus() == Passed ) {
00828             // we have convergence
00829             break;  // break from while(1){pcpg_iter->iterate()}
00830           }
00832           //
00833           // check for maximum iterations
00834           //
00836           else if ( maxIterTest_->getStatus() == Passed ) {
00837             // we don't have convergence
00838             isConverged = false;
00839             break;  // break from while(1){pcpg_iter->iterate()}
00840           }
00841           else {
00842 
00844           //
00845           // we returned from iterate(), but none of our status tests Passed.
00846           // Something is wrong, and it is probably the developers fault.
00847           //
00849 
00850             TEST_FOR_EXCEPTION(true,std::logic_error,
00851                                "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
00852           } // end if
00853         } // end try
00854         catch (const PCPGIterOrthoFailure &e) {
00855 
00856           // Check to see if the most recent solution yielded convergence.
00857           sTest_->checkStatus( &*pcpg_iter );
00858           if (convTest_->getStatus() != Passed)
00859             isConverged = false;
00860           break;
00861         }
00862         catch (const std::exception &e) {
00863           printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
00864                                    << pcpg_iter->getNumIters() << std::endl
00865                                    << e.what() << std::endl;
00866           throw;
00867         }
00868       } // end of while(1)
00869 
00870       // Update the linear problem.
00871       Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
00872       problem_->updateSolution( update, true );
00873 
00874       // Inform the linear problem that we are finished with this block linear system.
00875       problem_->setCurrLS();
00876 
00877       // Get the state.   How did pcpgState die?
00878       PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
00879 
00880       dimU_ = oldState.curDim;
00881       int q = oldState.prevUdim;
00882 
00883       std::cout << "SolverManager: dimU_ " << dimU_ << "   prevUdim= " << q << std::endl;
00884 
00885       if( q > deflatedBlocks_ ) 
00886         std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
00887 
00888       int rank; 
00889       if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
00890         //Given the seed space U and C = A U for some symmetric positive definite A,
00891         //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
00892 
00893         //oldState.D->print( std::cout ); D = diag( C'*U )
00894 
00895         U_ = oldState.U; //MVT::MvPrint( *U_, std::cout ); 
00896         C_ = oldState.C; //MVT::MvPrint( *C_, std::cout ); 
00897         rank = ARRQR(dimU_,q, *oldState.D );
00898         if( rank < dimU_ ) {
00899                 std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
00900         }  
00901         dimU_ = rank;
00902 
00903       } // Now U_ and C_ = AU are dual bases.
00904 
00905       if( dimU_ > deflatedBlocks_ ){
00906 
00907         if( !deflatedBlocks_ ){
00908            U_ = Teuchos::null;
00909            C_ = Teuchos::null;
00910            dimU_ = deflatedBlocks_;
00911            break;
00912         }
00913 
00914         bool Harmonic = false; // (Harmonic) Ritz vectors
00915      
00916         Teuchos::RCP<MV> Uorth;
00917 
00918         std::vector<int> active_cols( dimU_ ); 
00919         for (int i=0; i < dimU_; ++i) active_cols[i] = i;
00920 
00921         if( Harmonic ){
00922           Uorth = MVT::CloneCopy(*C_, active_cols); 
00923         }
00924         else{
00925           Uorth = MVT::CloneCopy(*U_, active_cols); 
00926         }
00927 
00928         // Explicitly construct Q and R factors 
00929         Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
00930         rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
00931         Uorth = Teuchos::null;
00932         // TODO:  During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.  
00933         // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
00934 
00935         // throw an error if U is both A-orthonormal and rank deficient
00936         TEST_FOR_EXCEPTION(rank < dimU_,PCPGSolMgrOrthoFailure,
00937                            "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
00938 
00939 
00940         // R VT' = Ur S,   
00941         Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
00942         Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
00943         int lwork = 5*dimU_;                           // minimal, extra computation < 67*dimU_
00944         int info = 0;  // Hermite
00945         int lrwork = 1;
00946         if( problem_->isHermitian() ) lrwork = dimU_;
00947         std::vector<ScalarType> work(lwork); // 
00948         std::vector<ScalarType> Svec(dimU_); // 
00949         std::vector<ScalarType> rwork(lrwork); 
00950         lapack.GESVD('N', 'O',
00951                    R.numRows(),R.numCols(),R.values(), R.numRows(),
00952                    &Svec[0],
00953                    Ur.values(),1,
00954                    VT.values(),1, // Output: VT stored in R
00955                    &work[0], lwork,
00956                    &rwork[0], &info);
00957 
00958         TEST_FOR_EXCEPTION(info != 0, PCPGSolMgrLAPACKFailure,
00959            "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
00960 
00961         if( work[0] !=  67. * dimU_ )
00962            std::cout << " SVD " << dimU_ <<  " lwork " << work[0]  << std::endl;
00963         for( int i=0; i< dimU_; i++)
00964            std::cout << i << " " << Svec[i] << std::endl;
00965            
00966         Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
00967 
00968         int startRow = 0, startCol = 0;
00969         if( Harmonic )
00970           startCol = dimU_ - deflatedBlocks_;
00971 
00972         Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
00973                                                      wholeV,
00974                                                      wholeV.numRows(),
00975                                                      deflatedBlocks_,
00976                                          startRow,
00977                                          startCol);
00978         std::vector<int> active_columns( dimU_ );
00979         std::vector<int> def_cols( deflatedBlocks_ );
00980         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
00981         for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
00982 
00983         Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
00984         Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
00985         MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); //  U:= U*V
00986         Ucopy   = Teuchos::null;
00987         Uactive = Teuchos::null;
00988         Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
00989         Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
00990         MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); //  C:= C*V
00991         Ccopy  = Teuchos::null;
00992         Cactive = Teuchos::null;
00993         dimU_ = deflatedBlocks_;
00994       }
00995       printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
00996 
00997       // Inform the linear problem that we are finished with this block linear system.
00998       problem_->setCurrLS();
00999 
01000       // Update indices for the linear systems to be solved.
01001       numRHS2Solve -= 1;
01002       if ( numRHS2Solve > 0 ) {
01003   currIdx[0]++;
01004 
01005         // Set the next indices.
01006         problem_->setLSIndex( currIdx );
01007       }
01008       else {
01009         currIdx.resize( numRHS2Solve );
01010       }
01011     }// while ( numRHS2Solve > 0 )
01012   }
01013     
01014   // print final summary
01015   sTest_->print( printer_->stream(FinalSummary) );
01016   
01017   // print timing information
01018   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01019 
01020   // get iteration information for this solve
01021   numIters_ = maxIterTest_->getNumIters();
01022  
01023   if (!isConverged) {
01024     return Unconverged; // return from PCPGSolMgr::solve() 
01025   }
01026   return Converged; // return from PCPGSolMgr::solve() 
01027 }
01028 
01029 // A-orthogonalize the Seed Space
01030 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
01031 // that are not rank revealing, and are not designed for PCPG in other ways too.
01032 template<class ScalarType, class MV, class OP>
01033 int PCPGSolMgr<ScalarType,MV,OP>::ARRQR(int p, int q, const Teuchos::SerialDenseMatrix<int,ScalarType>& D)
01034 {
01035   using Teuchos::RCP;
01036   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01037   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01038 
01039   // Allocate memory for scalars.
01040   Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
01041   Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
01042   Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
01043   std::vector<int> curind(1);
01044   std::vector<int> ipiv(p - q); // RRQR Pivot indices
01045   std::vector<ScalarType> Pivots(p); // RRQR Pivots
01046   int i, imax, j, k, l;
01047   ScalarType rteps = 1.5e-8;
01048 
01049   // Scale such that diag( U'C) = I
01050   for( int i = q ; i < p ; i++ ){
01051     ipiv[i-q] = i;
01052     curind[0] = i;
01053     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01054     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01055     anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
01056     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
01057     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
01058     Pivots[i]  = one;
01059   }
01060 
01061   for( i = q ; i < p ; i++ ){  
01062     if( q < i && i < p-1 ){ // Find the largest pivot
01063       imax = i;
01064       l = ipiv[imax-q];
01065       for( j = i+1 ; j < p ; j++ ){  
01066          k = ipiv[j-q];  
01067          if( Pivots[k] > Pivots[l] ){
01068            imax = j;  
01069            l = k;
01070          }
01071       } // end for
01072       if( imax > i ){ 
01073           l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
01074           ipiv[imax-q] = ipiv[i-q];
01075           ipiv[i-q] = l;
01076       }
01077     } // largest pivot found
01078     int k = ipiv[i-q];
01079  
01080     if( Pivots[k]  > 1.5625e-2 ){
01081       anorm(0,0) =  Pivots[k]; // A-norm of u
01082     }
01083     else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
01084       curind[0] = k;
01085       RCP<const MV> P = MVT::CloneView(*U_,curind);
01086       RCP<const MV> AP = MVT::CloneView(*C_,curind);
01087       MVT::MvTransMv( one, *P, *AP, anorm );
01088       anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
01089     }
01090     if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
01091        /*
01092        C(:,k) = A*U(:,k);  % Change C
01093        fixC = U(:, ipiv(1:i-1) )'*C(:,k);
01094        U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
01095        C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
01096        anorm = sqrt( U(:,k)'*C(:,k) );
01097        */
01098        std::cout << "ARRQR: Bad case not implemented" << std::endl;
01099     }
01100     if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
01101        std::cout << "ARRQR : deficient case not implemented " << std::endl;
01102        //U = U(:, ipiv(1:i-1) );
01103        //C = C(:, ipiv(1:i-1) );
01104        p = q + i; 
01105        // update curDim_ in State
01106        break; 
01107     }
01108     curind[0] = k;
01109     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01110     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01111     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
01112     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
01113     P = Teuchos::null;
01114     AP = Teuchos::null;
01115     Pivots[k] = one;                 // delete,  for diagonostic purposes
01116     P = MVT::CloneViewNonConst(*U_,curind);  // U(:,k)
01117     AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
01118     for( j = i+1 ; j < p ; j++ ){
01119       l = ipiv[j-q];   // ahhh
01120       curind[0] = l;
01121       RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault,  j=i+1=5
01122       MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
01123       MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
01124       Q = Teuchos::null;
01125       RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
01126       MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
01127       AQ = Teuchos::null;
01128       gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
01129       if( gamma(0,0) > 0){
01130         Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
01131       }
01132       else {
01133         Pivots[l] = zero; //rank deficiency revealed
01134       }
01135     }
01136   }
01137   return p;
01138 }
01139 
01140 //  The method returns a string describing the solver manager.
01141 template<class ScalarType, class MV, class OP>
01142 std::string PCPGSolMgr<ScalarType,MV,OP>::description() const
01143 {
01144   std::ostringstream oss;
01145   oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01146   oss << "{";
01147   oss << "Ortho Type='"<<orthoType_;
01148   oss << "}";
01149   return oss.str();
01150 }
01151   
01152 } // end Belos namespace
01153 
01154 #endif /* BELOS_PCPG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:06 2011 for Belos by  doxygen 1.6.3