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

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7