BelosGCRODRSolMgr.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_GCRODR_SOLMGR_HPP
00030 #define BELOS_GCRODR_SOLMGR_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041 
00042 #include "BelosBlockGmresIter.hpp"
00043 #include "BelosGCRODRIter.hpp"
00044 #include "BelosBlockFGmresIter.hpp"
00045 #include "BelosDGKSOrthoManager.hpp"
00046 #include "BelosICGSOrthoManager.hpp"
00047 #include "BelosIMGSOrthoManager.hpp"
00048 #include "BelosStatusTestMaxIters.hpp"
00049 #include "BelosStatusTestGenResNorm.hpp"
00050 #include "BelosStatusTestCombo.hpp"
00051 #include "BelosStatusTestOutput.hpp"
00052 #include "BelosOutputManager.hpp"
00053 #include "Teuchos_BLAS.hpp"
00054 #include "Teuchos_LAPACK.hpp"
00055 #include "Teuchos_TimeMonitor.hpp"
00056 
00072 namespace Belos {
00073   
00075 
00076   
00083   class GCRODRSolMgrLinearProblemFailure : public BelosError {public:
00084     GCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00085     {}};
00086   
00093   class GCRODRSolMgrOrthoFailure : public BelosError {public:
00094     GCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00095     {}};
00096   
00103   class GCRODRSolMgrLAPACKFailure : public BelosError {public:
00104     GCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00105     {}};
00106 
00113   class GCRODRSolMgrRecyclingFailure : public BelosError {public:
00114     GCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
00115     {}};
00116 
00118 
00119 
00120   template<class ScalarType, class MV, class OP>
00121   class GCRODRSolMgr : public SolverManager<ScalarType,MV,OP> {
00122     
00123   private:
00124     typedef MultiVecTraits<ScalarType,MV> MVT;
00125     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00126     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00127     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00128     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00129     
00130   public:
00131     
00133 
00134    
00140      GCRODRSolMgr();
00141  
00154     GCRODRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00155           const Teuchos::RCP<Teuchos::ParameterList> &pl );
00156     
00158     virtual ~GCRODRSolMgr() {};
00160     
00162 
00163     
00166     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00167       return *problem_;
00168     }
00169 
00172     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00173 
00176     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00177  
00183     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00184       return tuple(timerSolve_);
00185     }
00186   
00188     int getNumIters() const {
00189       return numIters_;
00190     }
00191  
00194     bool isLOADetected() const { return false; }
00195  
00197     
00199 
00200     
00201     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00202     
00203     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00204     
00206     
00208 
00209     
00227     ReturnType solve();
00228     
00230     
00233     
00235     std::string description() const;
00236     
00238     
00239   private:
00240 
00241     //  Computes harmonic eigenpairs of projected matrix created during the priming solve.
00242     //  HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
00243     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00244     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00245     int getHarmonicVecs1(int m, 
00246        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00247        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00248 
00249     //  Computes harmonic eigenpairs of projected matrix created during one cycle.
00250     //  HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
00251     //  VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
00252     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00253     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00254     int getHarmonicVecs2(int keff, int m, 
00255        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00256        const Teuchos::RCP<const MV>& VV,
00257        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00258 
00259     // Sort list of n floating-point numbers and return permutation vector
00260     void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00261 
00262     // Method to convert string to enumerated type for residual.
00263     Belos::ScaleType convertStringToScaleType( string& scaleType ) {
00264       if (scaleType == "Norm of Initial Residual") {
00265   return Belos::NormOfInitRes;
00266       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00267   return Belos::NormOfPrecInitRes;
00268       } else if (scaleType == "Norm of RHS") {
00269   return Belos::NormOfRHS;
00270       } else if (scaleType == "None") {
00271   return Belos::None;
00272       } else 
00273   TEST_FOR_EXCEPTION( true ,std::logic_error,
00274           "Belos::GCRODRSolMgr(): Invalid residual scaling type.");
00275     }
00276 
00277     // Linear problem.
00278     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00279     
00280     // Output manager.
00281     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00282     Teuchos::RCP<ostream> outputStream_;
00283 
00284     // Status test.
00285     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00286     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00287     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00288     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00289     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00290 
00291     // Orthogonalization manager.
00292     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00293     
00294     // Current parameter list.
00295     Teuchos::RCP<ParameterList> params_;
00296 
00297     // Default solver values.
00298     static const MagnitudeType convtol_default_;
00299     static const MagnitudeType orthoKappa_default_;
00300     static const int maxRestarts_default_;
00301     static const int maxIters_default_;
00302     static const int numBlocks_default_;
00303     static const int recycledBlocks_default_;
00304     static const int verbosity_default_;
00305     static const int outputFreq_default_;
00306     static const std::string impResScale_default_; 
00307     static const std::string expResScale_default_; 
00308     static const std::string label_default_;
00309     static const std::string orthoType_default_;
00310     static const Teuchos::RCP<ostream> outputStream_default_;
00311 
00312     // Current solver values.
00313     MagnitudeType convtol_, orthoKappa_;
00314     int maxRestarts_, maxIters_, numIters_;
00315     int numBlocks_, recycledBlocks_, verbosity_, outputFreq_;
00316     std::string orthoType_; 
00317     std::string impResScale_, expResScale_;
00318 
00319     // Recycled subspace and its image.
00320     Teuchos::RCP<MV> U_, C_, r_;
00321 
00322     // Timers.
00323     std::string label_;
00324     Teuchos::RCP<Teuchos::Time> timerSolve_;
00325 
00326     // Internal state variables.
00327     bool isSet_;
00328   };
00329 
00330 
00331 // Default solver values.
00332 template<class ScalarType, class MV, class OP>
00333 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GCRODRSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00334 
00335 template<class ScalarType, class MV, class OP>
00336 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType GCRODRSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00337 
00338 template<class ScalarType, class MV, class OP>
00339 const int GCRODRSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 20;
00340 
00341 template<class ScalarType, class MV, class OP>
00342 const int GCRODRSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00343 
00344 template<class ScalarType, class MV, class OP>
00345 const int GCRODRSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 300;
00346 
00347 template<class ScalarType, class MV, class OP>
00348 const int GCRODRSolMgr<ScalarType,MV,OP>::recycledBlocks_default_ = 5;
00349 
00350 template<class ScalarType, class MV, class OP>
00351 const int GCRODRSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00352 
00353 template<class ScalarType, class MV, class OP>
00354 const int GCRODRSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00355 
00356 template<class ScalarType, class MV, class OP>
00357 const std::string GCRODRSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00358 
00359 template<class ScalarType, class MV, class OP>
00360 const std::string GCRODRSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00361 
00362 template<class ScalarType, class MV, class OP>
00363 const std::string GCRODRSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00364 
00365 template<class ScalarType, class MV, class OP>
00366 const std::string GCRODRSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00367 
00368 template<class ScalarType, class MV, class OP>
00369 const Teuchos::RCP<ostream> GCRODRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00370 
00371 
00372 // Empty Constructor
00373 template<class ScalarType, class MV, class OP>
00374 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr() :
00375   outputStream_(outputStream_default_),
00376   convtol_(convtol_default_),
00377   orthoKappa_(orthoKappa_default_),
00378   maxRestarts_(maxRestarts_default_),
00379   maxIters_(maxIters_default_),
00380   numBlocks_(numBlocks_default_),
00381   recycledBlocks_(recycledBlocks_default_),
00382   verbosity_(verbosity_default_),
00383   outputFreq_(outputFreq_default_),
00384   orthoType_(orthoType_default_),
00385   impResScale_(impResScale_default_),
00386   expResScale_(expResScale_default_),
00387   label_(label_default_),
00388   isSet_(false)
00389 {}
00390 
00391 
00392 // Basic Constructor
00393 template<class ScalarType, class MV, class OP>
00394 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr( 
00395                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00396                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00397   problem_(problem),
00398   outputStream_(outputStream_default_),
00399   convtol_(convtol_default_),
00400   orthoKappa_(orthoKappa_default_),
00401   maxRestarts_(maxRestarts_default_),
00402   maxIters_(maxIters_default_),
00403   numBlocks_(numBlocks_default_),
00404   recycledBlocks_(recycledBlocks_default_),
00405   verbosity_(verbosity_default_),
00406   outputFreq_(outputFreq_default_),
00407   orthoType_(orthoType_default_),
00408   impResScale_(impResScale_default_),
00409   expResScale_(expResScale_default_),
00410   label_(label_default_),
00411   isSet_(false)
00412 {
00413   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00414 
00415   if (!is_null(pl)) {
00416     // Set the parameters using the list that was passed in.
00417     setParameters( pl );  
00418   }
00419 }
00420 
00421 
00422 template<class ScalarType, class MV, class OP>
00423 void GCRODRSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00424 {
00425   // Create the internal parameter list if ones doesn't already exist.
00426   if (params_ == Teuchos::null) {
00427     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00428   }
00429   else {
00430     params->validateParameters(*getValidParameters());
00431   }
00432 
00433   // Check for maximum number of restarts
00434   if (params->isParameter("Maximum Restarts")) {
00435     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00436 
00437     // Update parameter in our list.
00438     params_->set("Maximum Restarts", maxRestarts_);
00439   }
00440 
00441   // Check for maximum number of iterations
00442   if (params->isParameter("Maximum Iterations")) {
00443     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00444 
00445     // Update parameter in our list and in status test.
00446     params_->set("Maximum Iterations", maxIters_);
00447     if (maxIterTest_!=Teuchos::null)
00448       maxIterTest_->setMaxIters( maxIters_ );
00449   }
00450 
00451   // Check for the maximum number of blocks.
00452   if (params->isParameter("Num Blocks")) {
00453     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00454     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00455            "Belos::GCRODRSolMgr: \"Num Blocks\" must be strictly positive.");
00456 
00457     // Update parameter in our list.
00458     params_->set("Num Blocks", numBlocks_);
00459   }
00460 
00461   // Check for the maximum number of blocks.
00462   if (params->isParameter("Num Recycled Blocks")) {
00463     recycledBlocks_ = params->get("Num Recycled Blocks",recycledBlocks_default_);
00464     TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
00465            "Belos::GCRODRSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
00466 
00467     TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
00468            "Belos::GCRODRSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
00469 
00470     // Update parameter in our list.
00471     params_->set("Num Recycled Blocks", recycledBlocks_);
00472   }
00473 
00474   // Check to see if the timer label changed.
00475   if (params->isParameter("Timer Label")) {
00476     string tempLabel = params->get("Timer Label", label_default_);
00477 
00478     // Update parameter in our list and solver timer
00479     if (tempLabel != label_) {
00480       label_ = tempLabel;
00481       params_->set("Timer Label", label_);
00482       string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00483       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00484     }
00485   }
00486 
00487   // Check if the orthogonalization changed.
00488   if (params->isParameter("Orthogonalization")) {
00489     string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00490     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00491       std::invalid_argument,
00492       "Belos::GCRODRSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00493     if (tempOrthoType != orthoType_) {
00494       orthoType_ = tempOrthoType;
00495       // Create orthogonalization manager
00496       if (orthoType_=="DGKS") {
00497   if (orthoKappa_ <= 0) {
00498     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00499   }
00500   else {
00501     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00502     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00503   }
00504       }
00505       else if (orthoType_=="ICGS") {
00506   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00507       } 
00508       else if (orthoType_=="IMGS") {
00509   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00510       } 
00511     }  
00512   }
00513 
00514   // Check which orthogonalization constant to use.
00515   if (params->isParameter("Orthogonalization Constant")) {
00516     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00517 
00518     // Update parameter in our list.
00519     params_->set("Orthogonalization Constant",orthoKappa_);
00520     if (orthoType_=="DGKS") {
00521       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00522   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00523       }
00524     } 
00525   }
00526 
00527   // Check for a change in verbosity level
00528   if (params->isParameter("Verbosity")) {
00529     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00530       verbosity_ = params->get("Verbosity", verbosity_default_);
00531     } else {
00532       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00533     }
00534 
00535     // Update parameter in our list.
00536     params_->set("Verbosity", verbosity_);
00537     if (printer_ != Teuchos::null)
00538       printer_->setVerbosity(verbosity_);
00539   }
00540 
00541   // output stream
00542   if (params->isParameter("Output Stream")) {
00543     outputStream_ = Teuchos::getParameter<Teuchos::RCP<ostream> >(*params,"Output Stream");
00544 
00545     // Update parameter in our list.
00546     params_->set("Output Stream", outputStream_);
00547     if (printer_ != Teuchos::null)
00548       printer_->setOStream( outputStream_ );
00549   }
00550 
00551   // frequency level
00552   if (verbosity_ & Belos::StatusTestDetails) {
00553     if (params->isParameter("Output Frequency")) {
00554       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00555     }
00556 
00557     // Update parameter in out list and output status test.
00558     params_->set("Output Frequency", outputFreq_);
00559     if (outputTest_ != Teuchos::null)
00560       outputTest_->setOutputFrequency( outputFreq_ );
00561   }
00562 
00563   // Create output manager if we need to.
00564   if (printer_ == Teuchos::null) {
00565     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00566   }  
00567   
00568   // Convergence
00569   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00570   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00571 
00572   // Check for convergence tolerance
00573   if (params->isParameter("Convergence Tolerance")) {
00574     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00575 
00576     // Update parameter in our list and residual tests.
00577     params_->set("Convergence Tolerance", convtol_);
00578     if (impConvTest_ != Teuchos::null)
00579       impConvTest_->setTolerance( convtol_ );
00580     if (expConvTest_ != Teuchos::null)
00581       expConvTest_->setTolerance( convtol_ );
00582   }
00583  
00584   // Check for a change in scaling, if so we need to build new residual tests.
00585   if (params->isParameter("Implicit Residual Scaling")) {
00586     string tempImpResScale = Teuchos::getParameter<string>( *params, "Implicit Residual Scaling" );
00587 
00588     // Only update the scaling if it's different.
00589     if (impResScale_ != tempImpResScale) {
00590       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00591       impResScale_ = tempImpResScale;
00592 
00593       // Update parameter in our list and residual tests
00594       params_->set("Implicit Residual Scaling", impResScale_);
00595       if (impConvTest_ != Teuchos::null) {
00596         try { 
00597           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00598         }
00599         catch (std::exception& e) { 
00600           // Delete the convergence test so it gets constructed again.
00601     impConvTest_ = Teuchos::null;
00602           convTest_ = Teuchos::null;
00603         }
00604       }
00605     }      
00606   }
00607   
00608   if (params->isParameter("Explicit Residual Scaling")) {
00609     string tempExpResScale = Teuchos::getParameter<string>( *params, "Explicit Residual Scaling" );
00610 
00611     // Only update the scaling if it's different.
00612     if (expResScale_ != tempExpResScale) {
00613       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00614       expResScale_ = tempExpResScale;
00615 
00616       // Update parameter in our list and residual tests
00617       params_->set("Explicit Residual Scaling", expResScale_);
00618       if (expConvTest_ != Teuchos::null) {
00619         try { 
00620           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00621         }
00622         catch (std::exception& e) {
00623           // Delete the convergence test so it gets constructed again.
00624     expConvTest_ = Teuchos::null;
00625           convTest_ = Teuchos::null;
00626         }
00627       }
00628     }      
00629   }
00630 
00631   // Create status tests if we need to.
00632 
00633   // Basic test checks maximum iterations and native residual.
00634   if (maxIterTest_ == Teuchos::null)
00635     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00636 
00637   // Implicit residual test, using the native residual to determine if convergence was achieved.
00638   if (impConvTest_ == Teuchos::null) {
00639     impConvTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_ ) );
00640     impConvTest_->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00641   }
00642 
00643   // Explicit residual test once the native residual is below the tolerance
00644   if (expConvTest_ == Teuchos::null) {
00645     expConvTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_ ) );
00646     expConvTest_->defineResForm( StatusTestResNorm_t::Explicit, Belos::TwoNorm );
00647     expConvTest_->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00648   }
00649 
00650   if (convTest_ == Teuchos::null) {
00651     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00652   }
00653 
00654   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00655   
00656   if (outputFreq_ > 0) {
00657     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00658                           sTest_, 
00659                   outputFreq_, 
00660                   Passed+Failed+Undefined ) ); 
00661   }
00662   else {
00663     outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00664                     sTest_, 1 ) );
00665   }
00666 
00667   // Create orthogonalization manager if we need to.
00668   if (ortho_ == Teuchos::null) {
00669     if (orthoType_=="DGKS") {
00670       if (orthoKappa_ <= 0) {
00671   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00672       }
00673       else {
00674   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00675   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00676       }
00677     }
00678     else if (orthoType_=="ICGS") {
00679       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00680     } 
00681     else if (orthoType_=="IMGS") {
00682       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00683     } 
00684     else {
00685       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00686        "Belos::GCRODRSolMgr(): Invalid orthogonalization type.");
00687     }  
00688   }
00689 
00690   // Create the timer if we need to.
00691   if (timerSolve_ == Teuchos::null) {
00692     string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00693     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00694   }
00695 
00696   // Inform the solver manager that the current parameters were set.
00697   isSet_ = true;
00698 }
00699 
00700     
00701 template<class ScalarType, class MV, class OP>
00702 Teuchos::RCP<const Teuchos::ParameterList>
00703 GCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const
00704 {
00705   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00706   if (is_null(validPL)) {
00707     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00708     // Set all the valid parameters and their default values.
00709     pl->set("Convergence Tolerance", convtol_default_,
00710       "The relative residual tolerance that needs to be achieved by the\n"
00711       "iterative solver in order for the linear system to be declared converged.");
00712     pl->set("Maximum Restarts", maxRestarts_default_,
00713       "The maximum number of cycles allowed for each\n"
00714       "set of RHS solved.");
00715     pl->set("Maximum Iterations", maxIters_default_,
00716       "The maximum number of iterations allowed for each\n"
00717       "set of RHS solved.");
00718     pl->set("Num Blocks", numBlocks_default_,
00719       "The maximum number of vectors allowed in the Krylov subspace\n"
00720       "for each set of RHS solved.");
00721     pl->set("Num Recycled Blocks", recycledBlocks_default_,
00722       "The maximum number of vectors in the recycled subspace." );
00723     pl->set("Verbosity", verbosity_default_,
00724       "What type(s) of solver information should be outputted\n"
00725       "to the output stream.");
00726     pl->set("Output Frequency", outputFreq_default_,
00727       "How often convergence information should be outputted\n"
00728       "to the output stream.");  
00729     pl->set("Output Stream", outputStream_default_,
00730       "A reference-counted pointer to the output stream where all\n"
00731       "solver output is sent.");
00732     pl->set("Implicit Residual Scaling", impResScale_default_,
00733       "The type of scaling used in the implicit residual convergence test.");
00734     pl->set("Explicit Residual Scaling", expResScale_default_,
00735       "The type of scaling used in the explicit residual convergence test.");
00736     pl->set("Timer Label", label_default_,
00737       "The string to use as a prefix for the timer labels.");
00738     //  pl->set("Restart Timers", restartTimers_);
00739     pl->set("Orthogonalization", orthoType_default_,
00740       "The type of orthogonalization to use: DGKS, ICGS, IMGS");
00741     pl->set("Orthogonalization Constant",orthoKappa_default_,
00742       "The constant used by DGKS orthogonalization to determine\n"
00743       "whether another step of classical Gram-Schmidt is necessary.");
00744     validPL = pl;
00745   }
00746   return validPL;
00747 }
00748 
00749   
00750 // solve()
00751 template<class ScalarType, class MV, class OP>
00752 ReturnType GCRODRSolMgr<ScalarType,MV,OP>::solve() {
00753 
00754   // Set the current parameters if they were not set before.
00755   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00756   // then didn't set any parameters using setParameters().
00757   if (!isSet_) { setParameters( params_ ); }
00758 
00759   Teuchos::BLAS<int,ScalarType> blas;
00760   Teuchos::LAPACK<int,ScalarType> lapack;
00761   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00762   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00763   
00764   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure,
00765                      "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
00766 
00767   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,
00768                      "Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00769 
00770   // Create indices for the linear systems to be solved.
00771   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00772   std::vector<int> currIdx(1);
00773   currIdx[0] = 0;
00774 
00775   // Inform the linear problem of the current linear system to solve.
00776   problem_->setLSIndex( currIdx );
00777 
00779   int keff = 0;
00780   if (U_ != Teuchos::null) {
00781     keff = MVT::GetNumberVecs( *U_ );
00782   }
00783 
00784   // Check the number of blocks and change them is necessary.
00785   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
00786   if (numBlocks_ > dim) {
00787     numBlocks_ = dim;
00788     printer_->stream(Warnings) << 
00789       "Warning! Requested Krylov subspace dimension is larger that operator dimension!" << endl <<
00790       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << endl;
00791     params_->set("Num Blocks", numBlocks_);
00792   } 
00793 
00794   // Assume convergence is achieved, then let any failed convergence set this to false.
00795   bool isConverged = true;  
00796 
00798   // Parameter list
00799   Teuchos::ParameterList plist;
00800   
00801   plist.set("Num Blocks",numBlocks_);
00802   plist.set("Recycled Blocks",recycledBlocks_);
00803   
00805   // GCRODR solver
00806   
00807   Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
00808   gcrodr_iter = Teuchos::rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00809   // Number of iterations required to generate initial recycle space (if needed)
00810   int prime_iterations = 0;
00811 
00812   // Enter solve() iterations
00813   {
00814     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00815     
00816     while ( numRHS2Solve > 0 ) {
00817       
00818       // Reset the status test.  
00819       outputTest_->reset();
00820       
00822       // Initialize recycled subspace for GCRODR
00823       
00824       // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
00825       if (keff > 0) {
00826   TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
00827          "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
00828 
00829   printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
00830   // Compute image of U_ under the new operator
00831   std::vector<int> index(keff);
00832   for (int i=0; i<keff; ++i) { index[i] = i; }
00833   Teuchos::RCP<MV> Ckeff = MVT::CloneView( *C_, index );
00834   Teuchos::RCP<MV> Ukeff = MVT::CloneView( *U_, index );
00835   problem_->apply( *Ukeff, *Ckeff );
00836 
00837   // Orthogonalize this block
00838   // Get a matrix to hold the orthonormalization coefficients.
00839   Teuchos::SerialDenseMatrix<int,ScalarType> R(keff,keff);
00840   int rank = ortho_->normalize(*Ckeff, Teuchos::rcp(&R,false));
00841   
00842   // Throw an error if we could not orthogonalize this block
00843   TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,
00844          "Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
00845   
00846   // U_ = U_*R^{-1} 
00847   // First, compute LU factorization of R
00848   int info = 0;
00849   std::vector<int> ipiv(R.numRows());
00850   lapack.GETRF(R.numRows(),R.numCols(),R.values(),R.numRows(),&ipiv[0],&info);
00851   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
00852          "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
00853   
00854   // Now, form inv(R)
00855   int lwork = R.numRows();
00856   std::vector<ScalarType> work(lwork);
00857   lapack.GETRI(R.numRows(),R.values(),R.numRows(),&ipiv[0],&work[0],lwork,&info);
00858   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
00859          "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
00860   
00861   Teuchos::RCP<MV> tempU = MVT::Clone( *U_, keff );
00862   MVT::MvTimesMatAddMv( one, *U_, R, zero, *tempU );
00863   U_ = tempU;
00864 
00865   // Compute C_'*r_
00866   Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
00867   problem_->computeCurrPrecResVec( &*r_ );
00868   MVT::MvTransMv( one, *C_, *r_, Ctr );
00869 
00870   // Update solution ( x += U_*C_'*r_ )
00871   MVT::MvTimesMatAddMv( one, *U_, Ctr, one, *problem_->getCurrLHSVec() );
00872   
00873   // Update residual norm ( r -= C_*C_'*r_ )  
00874   MVT::MvTimesMatAddMv( -one, *C_, Ctr, one, *r_ );
00875 
00876         // We recycled space from previous call
00877         prime_iterations = 0;
00878 
00879       }
00880       else {
00881   
00882   // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
00883   printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
00884   
00885   Teuchos::ParameterList primeList;
00886   
00887   // Tell the block solver that the block size is one.
00888   primeList.set("Block Size",1);
00889   primeList.set("Num Blocks",numBlocks_);
00890   primeList.set("Keep Hessenberg", true);
00891   primeList.set("Initialize Hessenberg", true);
00892   
00893   //  Create Gmres iteration object to perform one cycle of Gmres.
00894   Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
00895   gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
00896   
00897   // Create the first block in the current Krylov basis (residual).
00898   if (r_ == Teuchos::null)
00899     r_ = MVT::Clone( *(problem_->getRHS()), 1 );
00900   problem_->computeCurrPrecResVec( &*r_ );
00901   Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *r_ );
00902 
00903   // Get a matrix to hold the orthonormalization coefficients.
00904   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 = 
00905     rcp( new Teuchos::SerialDenseVector<int,ScalarType>(1) );
00906   
00907   // Orthonormalize the new V_0
00908   int rank = ortho_->normalize( *V_0, z_0 );
00909   TEST_FOR_EXCEPTION(rank != 1,GCRODRSolMgrOrthoFailure,
00910          "Belos::GCRODRSolMgr::solve(): Failed to compute initial block of orthonormal vectors for priming solve.");
00911   
00912   // Set the new state and initialize the solver.
00913   GmresIterationState<ScalarType,MV> newstate;
00914   newstate.V = V_0;
00915   newstate.z = z_0;
00916   newstate.curDim = 0;
00917   gmres_iter->initializeGmres(newstate);
00918   
00919   // Perform one cycle of Gmres iteration
00920   bool primeConverged = false;
00921   try {
00922     gmres_iter->iterate();
00923 
00924     // Check convergence first
00925     if ( convTest_->getStatus() == Passed ) {
00926       // we have convergence
00927       primeConverged = true;
00928     }
00929   }
00930   catch (GmresIterationOrthoFailure e) {
00931     // Try to recover the most recent least-squares solution
00932     gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
00933     
00934     // Check to see if the most recent least-squares solution yielded convergence.
00935     sTest_->checkStatus( &*gmres_iter );
00936     if (convTest_->getStatus() == Passed)
00937       primeConverged = true;
00938   }
00939   catch (std::exception e) {
00940     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
00941            << gmres_iter->getNumIters() << endl 
00942            << e.what() << endl;
00943     throw;
00944   }
00945 
00946         // Record number of iterations in generating initial recycle spacec
00947         prime_iterations = gmres_iter->getNumIters();
00948            
00949   // Update the linear problem.
00950   Teuchos::RCP<MV> update = gmres_iter->getCurrentUpdate();
00951   problem_->updateSolution( update, true );
00952 
00953   // Get the state.
00954   GmresIterationState<ScalarType,MV> oldState = gmres_iter->getState();
00955   int p = oldState.curDim;
00956 
00957   // Compute harmonic Ritz vectors 
00958   // NOTE:  The storage for the harmonic Ritz vectors (PP) is made one column larger 
00959         //        just in case we split a complex conjugate pair.
00960   // NOTE:  Generate a recycled subspace only if we have enough vectors.  If we converged
00961   //        too early, move on to the next linear system and try to generate a subspace again.
00962   if (recycledBlocks_ < p+1) {
00963     int info = 0;
00964     Teuchos::SerialDenseMatrix<int,ScalarType> PP( p, recycledBlocks_+1 );  
00965     keff = getHarmonicVecs1( p, *oldState.H, PP );
00966     // We can generate a subspace to recycle and we know its size, so intialize U_ and C_;
00967     if (U_ == Teuchos::null) {
00968       U_ = MVT::Clone( *problem_->getRHS(), keff );
00969     }
00970     else {
00971       if (MVT::GetNumberVecs( *U_ ) < keff) {
00972         U_ = MVT::Clone( *problem_->getRHS(), keff );
00973       }
00974     }
00975     if (C_ == Teuchos::null) {
00976       C_ = MVT::Clone( *problem_->getRHS(), keff );
00977     }
00978     else {
00979       if (MVT::GetNumberVecs( *C_ ) < keff) {
00980         C_ = MVT::Clone( *problem_->getRHS(), keff );
00981       }
00982     }   
00983 
00984     // Form U (the subspace to recycle)
00985           // U = oldState.V(:,1:p) * PP;
00986     std::vector<int> index( p );
00987           for (int i=0; i < p; ++i) { index[i] = i; }
00988     Teuchos::RCP<const MV> Vp = MVT::CloneView( *oldState.V, index );
00989           index.resize(keff);  // keff <= p
00990     Teuchos::RCP<MV> Up = MVT::CloneView( *U_, index );
00991     const Teuchos::SerialDenseMatrix<int,ScalarType> PPview( Teuchos::View, PP, p, keff );
00992           MVT::MvTimesMatAddMv( one, *Vp, PPview, zero, *Up );
00993     Vp = null;
00994 
00995     // Form orthonormalized C and adjust U so that C = A*U
00996     // [Q, R] = qr(H*P);
00997     const Teuchos::SerialDenseMatrix<int,ScalarType> Hview( Teuchos::View, *oldState.H, p+1, p );
00998     Teuchos::SerialDenseMatrix<int,ScalarType> HP( p+1, keff );
00999     HP.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Hview, PPview, zero );
01000 
01001     // Workspace size query for QR factorization of HP (the worksize will be placed in tau[0])
01002           int lwork = -1;
01003     std::vector<ScalarType> tau(keff);
01004           lapack.GEQRF(HP.numRows(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01005            &tau[0],lwork,&info);
01006     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01007            "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01008 
01009           lwork = (int)tau[0];
01010     std::vector<ScalarType> work(lwork);
01011           lapack.GEQRF(HP.numRows(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01012            &work[0],lwork,&info);
01013     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01014            "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01015 
01016           // Explicitly construct Q and R factors 
01017     // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01018     Teuchos::SerialDenseMatrix<int,ScalarType> R( Teuchos::Copy, HP, keff, keff );
01019           for(int i=1;i<keff;i++) { for(int j=0;j<i;j++) R(i,j) = zero; }
01020           lapack.ORGQR(HP.numRows(),HP.numCols(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01021            &work[0],lwork,&info);
01022     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01023            "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01024 
01025     // Now compute C = V(:,1:p+1) * Q 
01026     index.resize( p+1 );
01027           for (int i=0; i < p+1; ++i) { index[i] = i; }
01028     Vp = MVT::CloneView( *oldState.V, index );
01029           index.resize(keff);  // keff <= p
01030           for (int i=0; i < keff; ++i) { index[i] = i; }
01031     Teuchos::RCP<MV> Cp = MVT::CloneView( *C_, index );
01032           MVT::MvTimesMatAddMv( one, *Vp, HP, zero, *Cp );
01033     Vp = null;
01034           Cp = null;
01035         
01036     // Finally, compute U_ = U_*R^{-1}  
01037     // First, compute LU factorization of R
01038     std::vector<int> ipiv(R.numRows());
01039     lapack.GETRF(R.numRows(),R.numCols(),R.values(),R.numRows(),&ipiv[0],&info);
01040     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01041            "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01042     
01043     // Now, form inv(R)
01044     lwork = R.numRows();
01045     work.resize(lwork);
01046     lapack.GETRI(R.numRows(),R.values(),R.numRows(),&ipiv[0],&work[0],lwork,&info);
01047     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01048            "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01049     
01050     Teuchos::RCP<MV> tempU = MVT::Clone( *U_, keff );
01051     MVT::MvTimesMatAddMv( one, *Up, R, zero, *tempU );
01052     U_ = tempU;
01053 
01054     printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01055   }
01056 
01057         // Return to outer loop if the priming solve converged, set the next linear system.
01058   if (primeConverged) {
01059     // Inform the linear problem that we are finished with this block linear system.
01060     problem_->setCurrLS();
01061 
01062           // Update indices for the linear systems to be solved.
01063           numRHS2Solve -= 1;
01064           if ( numRHS2Solve > 0 ) {
01065             currIdx[0]++;
01066  
01067             // Set the next indices.
01068             problem_->setLSIndex( currIdx );
01069           }
01070           else {
01071             currIdx.resize( numRHS2Solve );
01072           }
01073 
01074     continue;
01075   }
01076       } // if (keff > 0) ...
01077       
01078       // Prepare for the Gmres iterations with the recycled subspace.
01079 
01080       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01081       gcrodr_iter->setSize( keff, numBlocks_ );
01082       
01083       // Reset the number of iterations.
01084       gcrodr_iter->resetNumIters(prime_iterations);
01085 
01086       // Reset the number of calls that the status test output knows about.
01087       outputTest_->resetNumCalls();
01088 
01089       // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
01090       if (r_ == Teuchos::null)
01091   r_ = MVT::Clone( *(problem_->getRHS()), 1 );
01092       problem_->computeCurrPrecResVec( &*r_ );
01093 
01094       // Get a matrix to hold the orthonormalization coefficients.
01095       Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > z_0 = 
01096         rcp( new Teuchos::SerialDenseVector<int,ScalarType>(1) );
01097       
01098       // Orthonormalize the new r_
01099       int rank = ortho_->normalize( *r_, z_0 );
01100       TEST_FOR_EXCEPTION(rank != 1,GCRODRSolMgrOrthoFailure,
01101        "Belos::GCRODRSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
01102       
01103       // Set the new state and initialize the solver.
01104       GCRODRIterState<ScalarType,MV> newstate;
01105       newstate.V = r_;
01106       newstate.z = z_0;
01107       newstate.U = U_;
01108       newstate.C = C_;
01109       newstate.curDim = 0;
01110       gcrodr_iter->initialize(newstate);
01111       int numRestarts = 0;
01112       while(1) {
01113   
01114   // tell gcrodr_iter to iterate
01115   try {
01116           printf("********** Calling iterate...\n");
01117     gcrodr_iter->iterate();
01118     
01120     //
01121     // check convergence first
01122     //
01124     if ( convTest_->getStatus() == Passed ) {
01125       // we have convergence
01126       break;  // break from while(1){gcrodr_iter->iterate()}
01127     }
01129     //
01130     // check for maximum iterations
01131     //
01133     else if ( maxIterTest_->getStatus() == Passed ) {
01134       // we don't have convergence
01135       isConverged = false;
01136       break;  // break from while(1){gcrodr_iter->iterate()}
01137     }
01139     //
01140     // check for restarting, i.e. the subspace is full
01141     //
01143     else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
01144 
01145             // Update the recycled subspace even if we have hit the maximum number of restarts.
01146  
01147       // Get the state.
01148       GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
01149             int p = oldState.curDim;     
01150 
01151             // Update the linear problem.
01152             Teuchos::RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01153             problem_->updateSolution( update, true );
01154 
01155             // Take the norm of the recycled vectors
01156             std::vector<MagnitudeType> d(keff);
01157             MVT::MvNorm( *U_, d );
01158             for (int i=0; i<keff; ++i) {
01159         d[i] = one / d[i];
01160       }
01161       MVT::MvScale( *U_, d );
01162           
01163             // Construct the full block upper Hessenberg matrix
01164             Teuchos::SerialDenseMatrix<int,ScalarType> H2(p+keff+1, p+keff);
01165 
01166             // Insert D into the leading keff block of H2
01167             for (int i=0; i<keff; ++i) {
01168               H2(i,i) = d[i];
01169             }
01170             
01171       // Insert B into the upper right-hand keff by p block of H2
01172       Teuchos::SerialDenseMatrix<int,ScalarType> H2temp(Teuchos::View, H2, keff, p, 0, keff);
01173             H2temp.assign(*oldState.B);
01174 
01175             // Insert H into the lower p+1 by p block of H2
01176             Teuchos::SerialDenseMatrix<int,ScalarType> H2temp2(Teuchos::View, H2, p+1, p, keff, keff);
01177             H2temp2.assign(*oldState.H);
01178 
01179             // Compute the harmoic Ritz pairs for the generalized eigenproblem
01180       Teuchos::SerialDenseMatrix<int,ScalarType> PP( p+keff, recycledBlocks_+1 );  
01181       int keff_new = getHarmonicVecs2( keff, p, H2, oldState.V, PP );
01182       printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff_new << std::endl << std::endl;
01183 
01184       // Code to form new U, C
01185       // U = [U V(:,1:p)] * P; (in two steps)
01186       
01187       // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01188       Teuchos::RCP<MV> tempU = MVT::Clone( *U_, keff_new );
01189       Teuchos::SerialDenseMatrix<int,ScalarType> tempPP( Teuchos::View, PP, keff, keff_new );
01190       MVT::MvTimesMatAddMv( one, *U_, tempPP, zero, *tempU );
01191 
01192       // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01193       std::vector<int> index(p);
01194       for (int i=0; i < p; i++) { index[i] = i; }
01195       Teuchos::RCP<const MV> Vp = MVT::CloneView( *oldState.V, index );
01196       Teuchos::SerialDenseMatrix<int,ScalarType> tempPP2( Teuchos::View, PP, p, keff_new, keff );
01197       MVT::MvTimesMatAddMv( one, *Vp, tempPP2, one, *tempU );
01198      
01199       // Form HP = H*P
01200             Teuchos::SerialDenseMatrix<int,ScalarType> tempPP3( Teuchos::View, PP, p+keff, keff_new );
01201       Teuchos::SerialDenseMatrix<int,ScalarType> HP(H2.numRows(),tempPP3.numCols());
01202       HP.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,H2,tempPP3,zero);
01203 
01204       // Workspace size query for QR factorization of HP (the worksize will be placed in tau[0])
01205       int info = 0, lwork = -1;
01206       std::vector<ScalarType> tau(keff_new);
01207       lapack.GEQRF(HP.numRows(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01208        &tau[0],lwork,&info);
01209       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01210              "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01211       
01212       lwork = (int)tau[0];
01213       std::vector<ScalarType> work(lwork);
01214       lapack.GEQRF(HP.numRows(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01215        &work[0],lwork,&info);
01216       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01217              "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01218 
01219       // Explicitly construct Q and R factors 
01220       // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01221       Teuchos::SerialDenseMatrix<int,ScalarType> R( Teuchos::Copy, HP, keff_new, keff_new );
01222             for(int i=1;i<keff_new;i++) { for(int j=0;j<i;j++) R(i,j) = zero; }
01223       lapack.ORGQR(HP.numRows(),HP.numCols(),HP.numCols(),HP.values(),HP.numRows(),&tau[0],
01224        &work[0],lwork,&info);
01225       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01226              "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01227 
01228       // Form orthonormalized C and adjust U accordingly so that C = A*U      
01229       // C = [C V] * Q;
01230 
01231       // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01232       index.resize(keff_new);
01233       for (int i=0; i < keff_new; i++) { index[i] = i; }
01234       Teuchos::RCP<MV> tempC = MVT::Clone( *C_, keff_new );
01235       Teuchos::SerialDenseMatrix<int,ScalarType> tempQ( Teuchos::View, HP, keff, keff_new );
01236       MVT::MvTimesMatAddMv( one, *C_, tempQ, zero, *tempC );
01237 
01238       // Now compute C = V(:,1:p+1) * Q 
01239       index.resize( p+1 );
01240       for (int i=0; i < p+1; ++i) { index[i] = i; }
01241       Vp = MVT::CloneView( *oldState.V, index );
01242       Teuchos::SerialDenseMatrix<int,ScalarType> tempQ2( Teuchos::View, HP, p+1, keff_new, keff );
01243       MVT::MvTimesMatAddMv( one, *Vp, tempQ2, one, *tempC );
01244             C_ = tempC;     
01245 
01246       // Finally, compute U_ = U_*R^{-1}  
01247       // First, compute LU factorization of R
01248       std::vector<int> ipiv(R.numRows());
01249       lapack.GETRF(R.numRows(),R.numCols(),R.values(),R.numRows(),&ipiv[0],&info);
01250       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01251              "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01252       
01253       // Now, form inv(R)
01254       lwork = R.numRows();
01255       work.resize(lwork);
01256       lapack.GETRI(R.numRows(),R.values(),R.numRows(),&ipiv[0],&work[0],lwork,&info);
01257       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01258              "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01259     
01260       if (keff != keff_new) {
01261         U_ = MVT::Clone( *problem_->getRHS(), keff_new );
01262       }
01263       MVT::MvTimesMatAddMv( one, *tempU, R, zero, *U_ );
01264 
01265             // NOTE:  If we have hit the maximum number of restarts then QUIT! 
01266       if ( numRestarts >= maxRestarts_ ) {
01267         isConverged = false;
01268         break; // break from while(1){gcrodr_iter->iterate()}
01269       }
01270       numRestarts++;
01271       
01272       printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01273 
01274       // Compute the restart vector.
01275       if (r_ == Teuchos::null)
01276         r_ = MVT::Clone( *(problem_->getRHS()), 1 );
01277       problem_->computeCurrPrecResVec( &*r_ );
01278       
01279       // Orthonormalize the new r_, z_0 already exists
01280       rank = ortho_->normalize( *r_, z_0 );
01281       TEST_FOR_EXCEPTION(rank != 1,GCRODRSolMgrOrthoFailure,
01282              "Belos::GCRODRSolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
01283       
01284       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01285       keff = keff_new;
01286       gcrodr_iter->setSize( keff, numBlocks_ );
01287 
01288       // Set the new state and initialize the solver.
01289       GCRODRIterState<ScalarType,MV> restartState;
01290       restartState.V = r_;
01291       restartState.z = z_0;
01292       restartState.U = U_;
01293       restartState.C = C_;
01294       restartState.curDim = 0;
01295       gcrodr_iter->initialize(restartState);
01296       
01297     } // end of restarting
01298     
01300     //
01301     // we returned from iterate(), but none of our status tests Passed.
01302     // something is wrong, and it is probably our fault.
01303     //
01305     
01306     else {
01307       TEST_FOR_EXCEPTION(true,std::logic_error,
01308              "Belos::GCRODRSolMgr::solve(): Invalid return from GCRODRIter::iterate().");
01309     }
01310   }
01311         catch (GCRODRIterOrthoFailure e) {
01312     // Try to recover the most recent least-squares solution
01313     gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
01314     
01315     // Check to see if the most recent least-squares solution yielded convergence.
01316     sTest_->checkStatus( &*gcrodr_iter );
01317     if (convTest_->getStatus() != Passed)
01318       isConverged = false;
01319     break;
01320         }
01321         catch (std::exception e) {
01322     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01323                              << gcrodr_iter->getNumIters() << endl 
01324            << e.what() << endl;
01325           throw;
01326   }
01327       }
01328      
01329       // Compute the current solution.
01330       // Update the linear problem.
01331       Teuchos::RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01332       problem_->updateSolution( update, true );
01333 
01334       // Inform the linear problem that we are finished with this block linear system.
01335       problem_->setCurrLS();
01336 
01337       // Update indices for the linear systems to be solved.
01338       numRHS2Solve -= 1;
01339       if ( numRHS2Solve > 0 ) {
01340   currIdx[0]++;
01341 
01342         // Set the next indices.
01343         problem_->setLSIndex( currIdx );
01344       }
01345       else {
01346         currIdx.resize( numRHS2Solve );
01347       }
01348  
01349     }// while ( numRHS2Solve > 0 )
01350     
01351   }
01352   
01353   // print final summary
01354   sTest_->print( printer_->stream(FinalSummary) );
01355   
01356   // print timing information
01357   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01358  
01359   // get iteration information for this solve
01360   numIters_ = maxIterTest_->getNumIters();
01361  
01362   if (!isConverged) {
01363     return Unconverged; // return from GCRODRSolMgr::solve() 
01364   }
01365   return Converged; // return from GCRODRSolMgr::solve() 
01366 }
01367 
01368 
01369 //  Compute the harmonic eigenpairs of the projected, dense system.
01370 template<class ScalarType, class MV, class OP>
01371 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs1(int m, 
01372                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01373                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP)
01374 {
01375   int i, j;
01376   bool xtraVec = false;
01377   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01378   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01379 
01380   // The LAPACK interface
01381   Teuchos::LAPACK<int,ScalarType> lapack;
01382   
01383   // Real and imaginary eigenvalue components
01384   std::vector<MagnitudeType> wr(m), wi(m);
01385 
01386   // Real and imaginary (right) eigenvectors
01387   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m);
01388 
01389   // Magnitude of harmonic Ritz values
01390   std::vector<MagnitudeType> w(m);
01391 
01392   // Sorted order of harmonic Ritz values, also used for DGEEV
01393   std::vector<int> iperm(m);
01394 
01395   // Size of workspace and workspace for DGEEV
01396   int lwork = 4*m;
01397   std::vector<ScalarType> work(lwork);
01398 
01399   // Output info
01400   int info = 0;
01401 
01402   // Solve linear system:  H_m^{-H}*e_m 
01403   Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
01404   Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
01405   e_m[m-1] = one;
01406   lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
01407   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01408          "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01409 
01410   // Compute H_m + d*H_m^{-H}*e_m*e_m^H
01411   ScalarType d = HH(m, m-1) * HH(m, m-1);
01412   Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( HH );
01413   for( i=0; i<m; ++i ) {
01414     harmHH(i, m-1) += d * e_m[i];
01415   }
01416 
01417   // Revise to do query for optimal workspace first
01418   // Create simple storage for the left eigenvectors, which we don't care about.
01419   const int ldvl = m;
01420   ScalarType* vl = 0;
01421   lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
01422         vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01423   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01424          "Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01425 
01426   // Construct magnitude of each harmonic Ritz value
01427   for( i=0; i<m; ++i ) {
01428     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01429   }
01430 
01431   // Construct magnitude of each harmonic Ritz value
01432   this->sort(w, m, iperm);
01433 
01434   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01435   if (wi[iperm[recycledBlocks_-1]] != zero) {
01436     int countImag = 0;
01437     for ( i=0; i<recycledBlocks_; ++i ) {
01438       if (wi[iperm[i]] != zero)
01439   countImag++;
01440     }
01441     // Check to see if this count is even or odd:
01442     if (countImag % 2)
01443       xtraVec = true;
01444   }
01445 
01446   // Select recycledBlocks_ smallest eigenvectors
01447   for( i=0; i<recycledBlocks_; ++i ) {
01448     for( j=0; j<m; j++ ) {
01449       PP(j,i) = vr(j,iperm[i]);
01450     }
01451   }
01452   
01453   if (xtraVec) { // we need to store one more vector
01454     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
01455       for( j=0; j<m; ++j ) {   // so get the "imag" component
01456   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
01457       }
01458     }
01459     else { //  I picked the "imag" component
01460       for( j=0; j<m; ++j ) {   // so get the "real" component
01461   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
01462       }
01463     }
01464   }
01465 
01466   // Return whether we needed to store an additional vector
01467   if (xtraVec) {
01468     return recycledBlocks_+1;
01469   }
01470   return recycledBlocks_;
01471 }
01472 
01473 //  Compute the harmonic eigenpairs of the projected, dense system.
01474 template<class ScalarType, class MV, class OP>
01475 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs2(int keff, int m, 
01476                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01477                  const Teuchos::RCP<const MV>& VV,
01478                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP)
01479 {
01480   int i, j;
01481   int m2 = HH.numCols();
01482   bool xtraVec = false;
01483   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01484   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01485   std::vector<int> index;
01486   
01487   // The LAPACK interface
01488   Teuchos::LAPACK<int,ScalarType> lapack;
01489   
01490   // Real and imaginary eigenvalue components
01491   std::vector<ScalarType> wr(m2), wi(m2);
01492 
01493   // Magnitude of harmonic Ritz values
01494   std::vector<MagnitudeType> w(m2);
01495 
01496   // Real and imaginary (right) eigenvectors
01497   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2);
01498 
01499   // Sorted order of harmonic Ritz values
01500   std::vector<int> iperm(m2);
01501 
01502   // Form matrices for generalized eigenproblem
01503   
01504   // B = H2' * H2;
01505   Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2);
01506   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
01507   
01508   // A_tmp = | C'*U        0 |
01509   //         | V_{m+1}'*U  I |
01510   Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keff+m+1, keff+m );
01511 
01512   // A_tmp(1:keff,1:keff) = C' * U;
01513   Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keff, keff );
01514   MVT::MvTransMv( one, *C_, *U_, A11 );
01515 
01516   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
01517   Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keff, keff );
01518   index.resize(m+1);
01519   for (i=0; i < m+1; i++) { index[i] = i; }
01520   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
01521   MVT::MvTransMv( one, *Vp, *U_, A21 );
01522 
01523   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
01524   for( i=keff; i<keff+m; i++ ) {
01525     A_tmp(i,i) = one;
01526   }
01527 
01528   // A = H2' * A_tmp;
01529   Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
01530   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
01531 
01532   // Compute k smallest harmonic Ritz pairs
01533   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
01534   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
01535   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
01536   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
01537   // MLP : 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting (for now)
01538   char balanc='P', jobvl='N', jobvr='V', sense='N'; 
01539   int ld = A.numRows();
01540   int lwork = 6*ld;
01541   int ldvl = ld, ldvr = ld;
01542   int info = 0,ilo = 0,ihi = 0;
01543   ScalarType abnrm = zero, bbnrm = zero;
01544   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
01545   std::vector<ScalarType> beta(ld);
01546   std::vector<ScalarType> work(lwork);
01547   std::vector<MagnitudeType> lscale(ld), rscale(ld);  
01548   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
01549   std::vector<int> iwork(ld+6);
01550   int *bwork = 0; // If sense == 'N', bwork is never referenced
01551   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0], 
01552                &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0], 
01553                &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
01554   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,
01555          "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
01556   
01557   // Construct magnitude of each harmonic Ritz value
01558   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
01559   for( i=0; i<ld; i++ ) {
01560     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
01561   }
01562 
01563   // Construct magnitude of each harmonic Ritz value
01564   this->sort(w,ld,iperm);
01565 
01566   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01567   if (wi[iperm[ld-recycledBlocks_]] != zero) {
01568     int countImag = 0;
01569     for ( i=ld-recycledBlocks_; i<ld; i++ ) {
01570       if (wi[iperm[i]] != zero)
01571   countImag++;
01572     }
01573     // Check to see if this count is even or odd:
01574     if (countImag % 2)
01575       xtraVec = true;
01576   }
01577   
01578   // Select recycledBlocks_ smallest eigenvectors
01579   for( i=0; i<recycledBlocks_; i++ ) {
01580     for( j=0; j<ld; j++ ) {
01581       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
01582     }
01583   }
01584   
01585   if (xtraVec) { // we need to store one more vector
01586     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
01587       for( j=0; j<ld; j++ ) {   // so get the "imag" component
01588   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
01589       }
01590     }
01591     else { // I picked the "imag" component
01592       for( j=0; j<ld; j++ ) {   // so get the "real" component
01593   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
01594       }
01595     }
01596   }
01597   
01598   // Return whether we needed to store an additional vector
01599   if (xtraVec) {
01600     return recycledBlocks_+1;
01601   }
01602   return recycledBlocks_;
01603 }
01604 
01605 
01606 // This method sorts list of n floating-point numbers and return permutation vector
01607 template<class ScalarType, class MV, class OP>
01608 void GCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
01609 {
01610   int l, r, j, i, flag;
01611   int    RR2;
01612   double dRR, dK;
01613   
01614   // Initialize the permutation vector.
01615   for(j=0;j<n;j++)
01616     iperm[j] = j;
01617   
01618   if (n <= 1) return;
01619   
01620   l    = n / 2 + 1;
01621   r    = n - 1;
01622   l    = l - 1;
01623   dRR  = dlist[l - 1];
01624   dK   = dlist[l - 1];
01625   
01626   RR2 = iperm[l - 1];
01627   while (r != 0) {
01628     j = l;
01629     flag = 1;
01630     
01631     while (flag == 1) {
01632       i = j;
01633       j = j + j;
01634       
01635       if (j > r + 1)
01636   flag = 0;
01637       else {
01638   if (j < r + 1)
01639     if (dlist[j] > dlist[j - 1]) j = j + 1;
01640   
01641   if (dlist[j - 1] > dK) {
01642     dlist[i - 1] = dlist[j - 1];
01643     iperm[i - 1] = iperm[j - 1];
01644   }
01645   else {
01646     flag = 0;
01647   }
01648       }
01649     }
01650     dlist[i - 1] = dRR;
01651     iperm[i - 1] = RR2;
01652     
01653     if (l == 1) {
01654       dRR  = dlist [r];
01655       RR2 = iperm[r];
01656       dK = dlist[r];
01657       dlist[r] = dlist[0];
01658       iperm[r] = iperm[0];
01659       r = r - 1;
01660     }
01661     else {
01662       l   = l - 1;
01663       dRR  = dlist[l - 1];
01664       RR2  = iperm[l - 1];
01665       dK   = dlist[l - 1];
01666     }
01667   }
01668   dlist[0] = dRR;
01669   iperm[0] = RR2; 
01670 }
01671 
01672 
01673 //  This method requires the solver manager to return a string that describes itself.
01674 template<class ScalarType, class MV, class OP>
01675 std::string GCRODRSolMgr<ScalarType,MV,OP>::description() const
01676 {
01677   std::ostringstream oss;
01678   oss << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01679   oss << "{";
01680   oss << "Ortho Type='"<<orthoType_;
01681   oss << ", Num Blocks=" <<numBlocks_<< ", Max Restarts=" << maxRestarts_;
01682   oss << "}";
01683   return oss.str();
01684 }
01685   
01686 } // end Belos namespace
01687 
01688 #endif /* BELOS_GCRODR_SOLMGR_HPP */

Generated on Wed May 12 21:45:51 2010 for Belos by  doxygen 1.4.7