Belos Version of the Day
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_GCRODR_SOLMGR_HPP
00043 #define BELOS_GCRODR_SOLMGR_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosOrthoManagerFactory.hpp"
00052 
00053 #include "BelosLinearProblem.hpp"
00054 #include "BelosSolverManager.hpp"
00055 
00056 #include "BelosGCRODRIter.hpp"
00057 #include "BelosBlockFGmresIter.hpp"
00058 #include "BelosStatusTestMaxIters.hpp"
00059 #include "BelosStatusTestGenResNorm.hpp"
00060 #include "BelosStatusTestCombo.hpp"
00061 #include "BelosStatusTestOutputFactory.hpp"
00062 #include "BelosOutputManager.hpp"
00063 #include "Teuchos_BLAS.hpp"
00064 #include "Teuchos_LAPACK.hpp"
00065 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00066 #include "Teuchos_TimeMonitor.hpp"
00067 #endif // BELOS_TEUCHOS_TIME_MONITOR
00068 
00105 namespace Belos {
00106   
00108 
00109   
00116   class GCRODRSolMgrLinearProblemFailure : public BelosError {
00117     public:
00118       GCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
00119   };
00120   
00127   class GCRODRSolMgrOrthoFailure : public BelosError {
00128     public:
00129       GCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00130   };
00131   
00138   class GCRODRSolMgrLAPACKFailure : public BelosError {
00139     public:
00140       GCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
00141   };
00142 
00149   class GCRODRSolMgrRecyclingFailure : public BelosError {
00150     public:
00151       GCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
00152   };
00153 
00155 
00156   template<class ScalarType, class MV, class OP>
00157   class GCRODRSolMgr : public SolverManager<ScalarType,MV,OP> {
00158     
00159   private:
00160     typedef MultiVecTraits<ScalarType,MV> MVT;
00161     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00162     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00163     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00164     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00165     typedef OrthoManagerFactory<ScalarType, MV, OP> ortho_factory_type;
00166 
00167     
00168   public:
00170 
00171    
00177      GCRODRSolMgr();
00178  
00231     GCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00232       const Teuchos::RCP<Teuchos::ParameterList> &pl);
00233     
00235     virtual ~GCRODRSolMgr() {};
00237     
00239 
00240     
00243     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00244       return *problem_;
00245     }
00246 
00249     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00250 
00253     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { 
00254       return params_; 
00255     }
00256  
00262     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00263       return Teuchos::tuple(timerSolve_);
00264     }
00265 
00271     MagnitudeType achievedTol() const {
00272       return achievedTol_;
00273     }
00274   
00276     int getNumIters() const {
00277       return numIters_;
00278     }
00279  
00282     bool isLOADetected() const { return false; }
00283  
00285     
00287 
00288    
00290     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { 
00291       problem_ = problem; 
00292     }
00293    
00295     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00296     
00298    
00300 
00301 
00305     void reset( const ResetType type ) { 
00306       if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
00307       else if (type & Belos::RecycleSubspace) keff = 0;
00308     }
00310  
00312 
00313     
00330     ReturnType solve();
00331     
00333     
00336     
00338     std::string description() const;
00339     
00341     
00342   private:
00343 
00344     // Called by all constructors; Contains init instructions common to all constructors
00345     void init();
00346 
00347     // Initialize solver state storage
00348     void initializeStateStorage();
00349 
00350     // Compute updated recycle space given existing recycle space and newly generated Krylov space
00351     void buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
00352 
00353     //  Computes harmonic eigenpairs of projected matrix created during the priming solve.
00354     //  HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
00355     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00356     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00357     int getHarmonicVecs1(int m, 
00358        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00359        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00360 
00361     //  Computes harmonic eigenpairs of projected matrix created during one cycle.
00362     //  HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
00363     //  VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
00364     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00365     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00366     int getHarmonicVecs2(int keff, int m, 
00367        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00368        const Teuchos::RCP<const MV>& VV,
00369        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00370 
00371     // Sort list of n floating-point numbers and return permutation vector
00372     void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00373 
00374     // Lapack interface
00375     Teuchos::LAPACK<int,ScalarType> lapack;
00376 
00377     // Linear problem.
00378     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00379     
00380     // Output manager.
00381     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00382     Teuchos::RCP<std::ostream> outputStream_;
00383 
00384     // Status test.
00385     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00386     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00387     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00388     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00389     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00390 
00394     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00395     
00396     // Current parameter list.
00397     Teuchos::RCP<Teuchos::ParameterList> params_;
00398 
00399     // Default solver values.
00400     static const MagnitudeType convTol_default_;
00401     static const MagnitudeType orthoKappa_default_;
00402     static const int maxRestarts_default_;
00403     static const int maxIters_default_;
00404     static const int numBlocks_default_;
00405     static const int blockSize_default_;
00406     static const int recycledBlocks_default_;
00407     static const int verbosity_default_;
00408     static const int outputStyle_default_;
00409     static const int outputFreq_default_;
00410     static const std::string impResScale_default_; 
00411     static const std::string expResScale_default_; 
00412     static const std::string label_default_;
00413     static const std::string orthoType_default_;
00414     static const Teuchos::RCP<std::ostream> outputStream_default_;
00415 
00416     // Current solver values.
00417     MagnitudeType convTol_, orthoKappa_, achievedTol_;
00418     int maxRestarts_, maxIters_, numIters_;
00419     int verbosity_, outputStyle_, outputFreq_;
00420     std::string orthoType_; 
00421     std::string impResScale_, expResScale_;
00422 
00424     // Solver State Storage
00426     // 
00427     // The number of blocks and recycle blocks (m and k, respectively)
00428     int numBlocks_, recycledBlocks_;
00429     // Current size of recycled subspace 
00430     int keff;
00431     //
00432     // Residual vector
00433     Teuchos::RCP<MV> r_;
00434     //
00435     // Search space
00436     Teuchos::RCP<MV> V_;
00437     //
00438     // Recycled subspace and its image
00439     Teuchos::RCP<MV> U_, C_;
00440     //
00441     // Updated recycle space and its image
00442     Teuchos::RCP<MV> U1_, C1_;
00443     //
00444     // Storage used in constructing harmonic Ritz values/vectors
00445     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
00446     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00447     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00448     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
00449     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
00450     std::vector<ScalarType> tau_;
00451     std::vector<ScalarType> work_;
00452     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00453     std::vector<int> ipiv_;
00455 
00456     // Timers.
00457     std::string label_;
00458     Teuchos::RCP<Teuchos::Time> timerSolve_;
00459 
00460     // Internal state variables.
00461     bool isSet_;
00462 
00463     // Have we generated or regenerated a recycle space yet this solve?
00464     bool builtRecycleSpace_;
00465   };
00466 
00467 
00468 // Default solver values.
00469 template<class ScalarType, class MV, class OP>
00470 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::convTol_default_ = 1e-8;
00471 
00472 template<class ScalarType, class MV, class OP>
00473 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = 0.0;
00474 
00475 template<class ScalarType, class MV, class OP>
00476 const int GCRODRSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 100;
00477 
00478 template<class ScalarType, class MV, class OP>
00479 const int GCRODRSolMgr<ScalarType,MV,OP>::maxIters_default_ = 5000;
00480 
00481 template<class ScalarType, class MV, class OP>
00482 const int GCRODRSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 50;
00483 
00484 template<class ScalarType, class MV, class OP>
00485 const int GCRODRSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00486 
00487 template<class ScalarType, class MV, class OP>
00488 const int GCRODRSolMgr<ScalarType,MV,OP>::recycledBlocks_default_ = 5;
00489 
00490 template<class ScalarType, class MV, class OP>
00491 const int GCRODRSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00492 
00493 template<class ScalarType, class MV, class OP>
00494 const int GCRODRSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00495 
00496 template<class ScalarType, class MV, class OP>
00497 const int GCRODRSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00498 
00499 template<class ScalarType, class MV, class OP>
00500 const std::string GCRODRSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00501 
00502 template<class ScalarType, class MV, class OP>
00503 const std::string GCRODRSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00504 
00505 template<class ScalarType, class MV, class OP>
00506 const std::string GCRODRSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00507 
00508 template<class ScalarType, class MV, class OP>
00509 const std::string GCRODRSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00510 
00511 template<class ScalarType, class MV, class OP>
00512 const Teuchos::RCP<std::ostream> GCRODRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00513 
00514 
00515 // Empty Constructor
00516 template<class ScalarType, class MV, class OP>
00517 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr() {
00518   init();
00519 }
00520 
00521 
00522 // Basic Constructor
00523 template<class ScalarType, class MV, class OP>
00524 GCRODRSolMgr<ScalarType,MV,OP>::
00525 GCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00526        const Teuchos::RCP<Teuchos::ParameterList> &pl ) 
00527 {
00528   // Initialize local pointers to null, and initialize local variables
00529   // to default values.
00530   init();
00531 
00532   TEUCHOS_TEST_FOR_EXCEPTION(problem == Teuchos::null, std::invalid_argument, 
00533          "Belos::GCRODRSolMgr constructor: The solver manager's "
00534          "constructor needs the linear problem argument 'problem' "
00535          "to be non-null.");
00536   problem_ = problem;
00537 
00538   // Set the parameters using the list that was passed in.  If null,
00539   // we defer initialization until a non-null list is set (by the
00540   // client calling setParameters(), or by calling solve() -- in
00541   // either case, a null parameter list indicates that default
00542   // parameters should be used).
00543   if (! is_null (pl))
00544     setParameters (pl);
00545 }
00546 
00547 // Common instructions executed in all constructors
00548 template<class ScalarType, class MV, class OP>
00549 void GCRODRSolMgr<ScalarType,MV,OP>::init() {
00550   outputStream_ = outputStream_default_;
00551   convTol_ = convTol_default_;
00552   orthoKappa_ = orthoKappa_default_;
00553   maxRestarts_ = maxRestarts_default_;
00554   maxIters_ = maxIters_default_;
00555   numBlocks_ = numBlocks_default_;
00556   recycledBlocks_ = recycledBlocks_default_;
00557   verbosity_ = verbosity_default_;
00558   outputStyle_ = outputStyle_default_;
00559   outputFreq_ = outputFreq_default_;
00560   orthoType_ = orthoType_default_;
00561   impResScale_ = impResScale_default_;
00562   expResScale_ = expResScale_default_;
00563   label_ = label_default_;
00564   isSet_ = false;
00565   builtRecycleSpace_ = false;
00566   keff = 0;
00567   r_ = Teuchos::null;
00568   V_ = Teuchos::null;
00569   U_ = Teuchos::null; 
00570   C_ = Teuchos::null;
00571   U1_ = Teuchos::null;
00572   C1_ = Teuchos::null;
00573   PP_ = Teuchos::null;
00574   HP_ = Teuchos::null;
00575   H2_ = Teuchos::null;
00576   R_ = Teuchos::null;
00577   H_ = Teuchos::null;
00578   B_ = Teuchos::null;
00579 }
00580 
00581 template<class ScalarType, class MV, class OP>
00582 void 
00583 GCRODRSolMgr<ScalarType,MV,OP>::
00584 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00585 {
00586   using Teuchos::isParameterType;
00587   using Teuchos::getParameter;
00588   using Teuchos::null;
00589   using Teuchos::ParameterList;
00590   using Teuchos::parameterList;
00591   using Teuchos::RCP;
00592   using Teuchos::rcp;
00593   using Teuchos::rcp_dynamic_cast;
00594   using Teuchos::rcpFromRef;
00595   using Teuchos::Exceptions::InvalidParameter;
00596   using Teuchos::Exceptions::InvalidParameterName;
00597   using Teuchos::Exceptions::InvalidParameterType;
00598 
00599   // The default parameter list contains all parameters that
00600   // GCRODRSolMgr understands, and none that it doesn't understand.
00601   RCP<const ParameterList> defaultParams = getValidParameters();
00602 
00603   // Create the internal parameter list if one doesn't already exist.
00604   //
00605   // (mfh 28 Feb 2011, 10 Mar 2011) At the time this code was written,
00606   // ParameterList did not have validators or validateParameters().
00607   // This is why the code below carefully validates the parameters one
00608   // by one and fills in defaults.  This code could be made a lot
00609   // shorter by using validators.  To do so, we would have to define
00610   // appropriate validators for all the parameters.  (This would more
00611   // or less just move all that validation code out of this routine
00612   // into to getValidParameters().)
00613   //
00614   // For an analogous reason, GCRODRSolMgr defines default parameter
00615   // values as class data, as well as in the default ParameterList.
00616   // This redundancy could be removed by defining the default
00617   // parameter values only in the default ParameterList (which
00618   // documents each parameter as well -- handy!).
00619   if (params_.is_null()) {
00620     params_ = parameterList (*defaultParams);
00621   } else {
00622     // A common case for setParameters() is for it to be called at the
00623     // beginning of the solve() routine.  This follows the Belos
00624     // pattern of delaying initialization until the last possible
00625     // moment (when the user asks Belos to perform the solve).  In
00626     // this common case, we save ourselves a deep copy of the input
00627     // parameter list.
00628     if (params_ != params) {
00629       // Make a deep copy of the input parameter list.  This allows
00630       // the caller to modify or change params later, without
00631       // affecting the behavior of this solver.  This solver will then
00632       // only change its internal parameters if setParameters() is
00633       // called again.
00634       params_ = parameterList (*params);
00635     }
00636 
00637     // Fill in any missing parameters and their default values.  Also,
00638     // throw an exception if the parameter list has any misspelled or
00639     // "extra" parameters.  If you don't like this behavior, you'll
00640     // want to replace the line of code below with your desired
00641     // validation scheme.  Note that Teuchos currently only implements
00642     // two options: 
00643     //
00644     // 1. validateParameters() requires that params_ has all the
00645     //    parameters that the default list has, and none that it
00646     //    doesn't have.
00647     //
00648     // 2. validateParametersAndSetDefaults() fills in missing
00649     //    parameters in params_ using the default list, but requires
00650     //    that any parameter provided in params_ is also in the
00651     //    default list.
00652     //
00653     // Here is an easy way to ignore any "extra" or misspelled
00654     // parameters: Make a deep copy of the default list, fill in any
00655     // "missing" parameters from the _input_ list, and then validate
00656     // the input list using the deep copy of the default list.  We
00657     // show this in code:
00658     //
00659     // RCP<ParameterList> defaultCopy = parameterList (*getValidParameters ());
00660     // defaultCopy->validateParametersAndSetDefaults (params);
00661     // params->validateParametersAndSetDefaults (defaultCopy);
00662     //
00663     // This method is not entirely robust, because the input list may
00664     // have incorrect validators set for existing parameters in the
00665     // default list.  This would then cause "validation" of the
00666     // default list to throw an exception.  As a result, we've chosen
00667     // for now to be intolerant of misspellings and "extra" parameters
00668     // in the input list.
00669     params_->validateParametersAndSetDefaults (*defaultParams);
00670   }
00671 
00672   // Check for maximum number of restarts.
00673   if (params->isParameter ("Maximum Restarts")) {
00674     maxRestarts_ = params->get("Maximum Restarts", maxRestarts_default_);
00675 
00676     // Update parameter in our list.
00677     params_->set ("Maximum Restarts", maxRestarts_);
00678   }
00679 
00680   // Check for maximum number of iterations
00681   if (params->isParameter ("Maximum Iterations")) {
00682     maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
00683 
00684     // Update parameter in our list and in status test.
00685     params_->set ("Maximum Iterations", maxIters_);
00686     if (! maxIterTest_.is_null())
00687       maxIterTest_->setMaxIters (maxIters_);
00688   }
00689 
00690   // Check for the maximum number of blocks.
00691   if (params->isParameter ("Num Blocks")) {
00692     numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
00693     TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00694            "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
00695            "be strictly positive, but you specified a value of "
00696            << numBlocks_ << ".");
00697     // Update parameter in our list.
00698     params_->set ("Num Blocks", numBlocks_);
00699   }
00700 
00701   // Check for the maximum number of blocks.
00702   if (params->isParameter ("Num Recycled Blocks")) {
00703     recycledBlocks_ = params->get ("Num Recycled Blocks", 
00704            recycledBlocks_default_);
00705     TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
00706            "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
00707            "parameter must be strictly positive, but you specified "
00708            "a value of " << recycledBlocks_ << ".");
00709     TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
00710            "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
00711            "parameter must be less than the \"Num Blocks\" "
00712            "parameter, but you specified \"Num Recycled Blocks\" "
00713            "= " << recycledBlocks_ << " and \"Num Blocks\" = " 
00714            << numBlocks_ << ".");
00715     // Update parameter in our list.
00716     params_->set("Num Recycled Blocks", recycledBlocks_);
00717   }
00718 
00719   // Check to see if the timer label changed.  If it did, update it in
00720   // the parameter list, and create a new timer with that label (if
00721   // Belos was compiled with timers enabled).
00722   if (params->isParameter ("Timer Label")) {
00723     std::string tempLabel = params->get ("Timer Label", label_default_);
00724 
00725     // Update parameter in our list and solver timer
00726     if (tempLabel != label_) {
00727       label_ = tempLabel;
00728       params_->set ("Timer Label", label_);
00729       std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00730 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00731       timerSolve_ = Teuchos::TimeMonitor::getNewTimer (solveLabel);
00732 #endif
00733     }
00734   }
00735 
00736   // Check for a change in verbosity level
00737   if (params->isParameter ("Verbosity")) {
00738     if (isParameterType<int> (*params, "Verbosity")) {
00739       verbosity_ = params->get ("Verbosity", verbosity_default_);
00740     } else {
00741       verbosity_ = (int) getParameter<Belos::MsgType> (*params, "Verbosity");
00742     }
00743     // Update parameter in our list.
00744     params_->set ("Verbosity", verbosity_);
00745     // If the output manager (printer_) is null, then we will
00746     // instantiate it later with the correct verbosity.
00747     if (! printer_.is_null())
00748       printer_->setVerbosity (verbosity_);
00749   }
00750 
00751   // Check for a change in output style
00752   if (params->isParameter ("Output Style")) {
00753     if (isParameterType<int> (*params, "Output Style")) {
00754       outputStyle_ = params->get ("Output Style", outputStyle_default_);
00755     } else {
00756       outputStyle_ = (int) getParameter<OutputType> (*params, "Output Style");
00757     }
00758 
00759     // Update parameter in our list.
00760     params_->set ("Output Style", outputStyle_);
00761     // We will (re)instantiate the output status test afresh below.
00762     outputTest_ = null;
00763   }
00764 
00765   // Get the output stream for the output manager.
00766   //
00767   // While storing the output stream in the parameter list (either as
00768   // an RCP or as a nonconst reference) is convenient and safe for
00769   // programming, it makes it impossible to serialize the parameter
00770   // list, read it back in from the serialized representation, and get
00771   // the same output stream as before.  This is because output streams
00772   // may be arbitrary constructed objects.
00773   //
00774   // In case the user tried reading in the parameter list from a
00775   // serialized representation and the output stream can't be read
00776   // back in, we set the output stream to point to std::cout.  This
00777   // ensures reasonable behavior.
00778   if (params->isParameter ("Output Stream")) {
00779     try {
00780       outputStream_ = getParameter<RCP<std::ostream> > (*params, "Output Stream");
00781     } catch (InvalidParameter&) {
00782       outputStream_ = rcpFromRef (std::cout);
00783     }
00784     // We assume that a null output stream indicates that the user
00785     // doesn't want to print anything, so we replace it with a "black
00786     // hole" stream that prints nothing sent to it.  (We can't use a
00787     // null output stream, since the output manager always sends
00788     // things it wants to print to the output stream.)
00789     if (outputStream_.is_null()) {
00790       outputStream_ = rcp (new Teuchos::oblackholestream);
00791     }
00792     // Update parameter in our list.
00793     params_->set ("Output Stream", outputStream_);
00794     // If the output manager (printer_) is null, then we will
00795     // instantiate it later with the correct output stream.
00796     if (! printer_.is_null()) {
00797       printer_->setOStream (outputStream_);
00798     }
00799   }
00800 
00801   // frequency level
00802   if (verbosity_ & Belos::StatusTestDetails) {
00803     if (params->isParameter ("Output Frequency")) {
00804       outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
00805     }
00806 
00807     // Update parameter in out list and output status test.
00808     params_->set("Output Frequency", outputFreq_);
00809     if (! outputTest_.is_null())
00810       outputTest_->setOutputFrequency (outputFreq_);
00811   }
00812 
00813   // Create output manager if we need to, using the verbosity level
00814   // and output stream that we fetched above.  We do this here because
00815   // instantiating an OrthoManager using OrthoManagerFactory requires
00816   // a valid OutputManager.
00817   if (printer_.is_null()) {
00818     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00819   }  
00820 
00821   // Get the orthogonalization manager name ("Orthogonalization").
00822   //
00823   // Getting default values for the orthogonalization manager
00824   // parameters ("Orthogonalization Parameters") requires knowing the
00825   // orthogonalization manager name.  Save it for later, and also
00826   // record whether it's different than before.
00827   OrthoManagerFactory<ScalarType, MV, OP> factory;
00828   bool changedOrthoType = false;
00829   if (params->isParameter ("Orthogonalization")) {
00830     const std::string& tempOrthoType = 
00831       params->get ("Orthogonalization", orthoType_default_);
00832     // Ensure that the specified orthogonalization type is valid.
00833     if (! factory.isValidName (tempOrthoType)) {
00834       std::ostringstream os;
00835       os << "Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 
00836    << tempOrthoType << "\".  The following are valid options "
00837    << "for the \"Orthogonalization\" name parameter: ";
00838       factory.printValidNames (os);
00839       throw std::invalid_argument (os.str());
00840     }
00841     if (tempOrthoType != orthoType_) {
00842       changedOrthoType = true;
00843       orthoType_ = tempOrthoType;
00844       // Update parameter in our list.
00845       params_->set ("Orthogonalization", orthoType_);
00846     }
00847   }
00848 
00849   // Get any parameters for the orthogonalization ("Orthogonalization
00850   // Parameters").  If not supplied, the orthogonalization manager
00851   // factory will supply default values.
00852   //
00853   // NOTE (mfh 12 Jan 2011) For the sake of backwards compatibility,
00854   // if params has an "Orthogonalization Constant" parameter and the
00855   // DGKS orthogonalization manager is to be used, the value of this
00856   // parameter will override DGKS's "depTol" parameter.
00857   //
00858   // Users must supply the orthogonalization manager parameters as a
00859   // sublist (supplying it as an RCP<ParameterList> would make the
00860   // resulting parameter list not serializable).
00861   RCP<ParameterList> orthoParams;
00862   { // The nonmember function sublist() returns an RCP<ParameterList>,
00863     // which is what we want here.
00864     using Teuchos::sublist;
00865     // Abbreviation to avoid typos.
00866     const std::string paramName ("Orthogonalization Parameters");
00867 
00868     try {
00869       orthoParams = sublist (params_, paramName, true);
00870     } catch (InvalidParameter&) {
00871       // We didn't get the parameter list from params, so get a
00872       // default parameter list from the OrthoManagerFactory.  Modify
00873       // params_ so that it has the default parameter list, and set
00874       // orthoParams to ensure it's a sublist of params_ (and not just
00875       // a copy of one).
00876       params_->set (paramName, factory.getDefaultParameters (orthoType_));
00877       orthoParams = sublist (params_, paramName, true);
00878     }
00879   }
00880   TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error, 
00881            "Failed to get orthogonalization parameters.  "
00882            "Please report this bug to the Belos developers.");
00883 
00884   // Instantiate a new MatOrthoManager subclass instance if necessary.
00885   // If not necessary, then tell the existing instance about the new
00886   // parameters.
00887   if (ortho_.is_null() || changedOrthoType) {
00888     // We definitely need to make a new MatOrthoManager, since either
00889     // we haven't made one yet, or we've changed orthogonalization
00890     // methods.  Creating the orthogonalization manager requires that
00891     // the OutputManager (printer_) already be initialized.
00892     ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_, 
00893             label_, orthoParams);
00894   } else {
00895     // If the MatOrthoManager implements the ParameterListAcceptor
00896     // mix-in interface, we can propagate changes to its parameters
00897     // without reinstantiating the MatOrthoManager.
00898     //
00899     // We recommend that all MatOrthoManager subclasses implement
00900     // Teuchos::ParameterListAcceptor, but do not (yet) require this.
00901     typedef Teuchos::ParameterListAcceptor PLA;
00902     RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
00903     if (pla.is_null()) {
00904       // Oops, it's not a ParameterListAcceptor.  We have to
00905       // reinstantiate the MatOrthoManager in order to pass in the
00906       // possibly new parameters.
00907       ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_,
00908               label_, orthoParams);
00909     } else {
00910       pla->setParameterList (orthoParams);
00911     }
00912   }
00913 
00914   // The DGKS orthogonalization accepts a "Orthogonalization Constant"
00915   // parameter (also called kappa in the code, but not in the
00916   // parameter list).  If its value is provided in the given parameter
00917   // list, and its value is positive, use it.  Ignore negative values.
00918   //
00919   // NOTE (mfh 12 Jan 2011) This overrides the "depTol" parameter that
00920   // may have been specified in "Orthogonalization Parameters".  We
00921   // retain this behavior for backwards compatibility.
00922   bool gotValidOrthoKappa = false;
00923   if (params->isParameter ("Orthogonalization Constant")) {
00924     const MagnitudeType orthoKappa = 
00925       params->get ("Orthogonalization Constant", orthoKappa_default_);
00926     if (orthoKappa > 0) {
00927       orthoKappa_ = orthoKappa;
00928       gotValidOrthoKappa = true;
00929       // Update parameter in our list.
00930       params_->set("Orthogonalization Constant", orthoKappa_);
00931       // Only DGKS currently accepts this parameter.
00932       if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
00933   typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
00934   // This cast should always succeed; it's a bug
00935   // otherwise.  (If the cast fails, then orthoType_
00936   // doesn't correspond to the OrthoManager subclass
00937   // instance that we think we have, so we initialized the
00938   // wrong subclass somehow.)
00939   rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
00940       }
00941     }
00942   }
00943   
00944   // Convergence
00945   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00946   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00947 
00948   // Check for convergence tolerance
00949   if (params->isParameter("Convergence Tolerance")) {
00950     convTol_ = params->get ("Convergence Tolerance", convTol_default_);
00951 
00952     // Update parameter in our list and residual tests.
00953     params_->set ("Convergence Tolerance", convTol_);
00954     if (! impConvTest_.is_null())
00955       impConvTest_->setTolerance (convTol_);
00956     if (! expConvTest_.is_null())
00957       expConvTest_->setTolerance (convTol_);
00958   }
00959  
00960   // Check for a change in scaling, if so we need to build new residual tests.
00961   if (params->isParameter ("Implicit Residual Scaling")) {
00962     std::string tempImpResScale = 
00963       getParameter<std::string> (*params, "Implicit Residual Scaling");
00964 
00965     // Only update the scaling if it's different.
00966     if (impResScale_ != tempImpResScale) {
00967       ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
00968       impResScale_ = tempImpResScale;
00969 
00970       // Update parameter in our list and residual tests
00971       params_->set("Implicit Residual Scaling", impResScale_);
00972       // NOTE (mfh 28 Feb 2011) StatusTestImpResNorm only lets you
00973       // call defineScaleForm() once.  The code below attempts to call
00974       // defineScaleForm(); if the scale form has already been
00975       // defined, it constructs a new StatusTestImpResNorm instance.
00976       // StatusTestImpResNorm really should not expose the
00977       // defineScaleForm() method, since it's serving an
00978       // initialization purpose; all initialization should happen in
00979       // the constructor whenever possible.  In that case, the code
00980       // below could be simplified into a single (re)instantiation.
00981       if (! impConvTest_.is_null()) {
00982         try { 
00983           impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
00984         }
00985         catch (StatusTestError&) {
00986           // Delete the convergence test so it gets constructed again.
00987     impConvTest_ = null;
00988           convTest_ = null;
00989         }
00990       }
00991     }      
00992   }
00993   
00994   if (params->isParameter("Explicit Residual Scaling")) {
00995     std::string tempExpResScale = 
00996       getParameter<std::string> (*params, "Explicit Residual Scaling");
00997 
00998     // Only update the scaling if it's different.
00999     if (expResScale_ != tempExpResScale) {
01000       ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
01001       expResScale_ = tempExpResScale;
01002 
01003       // Update parameter in our list and residual tests
01004       params_->set("Explicit Residual Scaling", expResScale_);
01005       // NOTE (mfh 28 Feb 2011) See note above on the (re)construction
01006       // of StatusTestImpResNorm.
01007       if (! expConvTest_.is_null()) {
01008         try { 
01009           expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
01010         }
01011         catch (StatusTestError&) {
01012           // Delete the convergence test so it gets constructed again.
01013     expConvTest_ = null;
01014           convTest_ = null;
01015         }
01016       }
01017     }      
01018   }
01019   //
01020   // Create iteration stopping criteria ("status tests") if we need
01021   // to, by combining three different stopping criteria.
01022   //
01023   // First, construct maximum-number-of-iterations stopping criterion.
01024   if (maxIterTest_.is_null())
01025     maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
01026 
01027   // Implicit residual test, using the native residual to determine if
01028   // convergence was achieved.
01029   if (impConvTest_.is_null()) {
01030     impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01031     impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), 
01032            Belos::TwoNorm);
01033   }
01034 
01035   // Explicit residual test once the native residual is below the tolerance
01036   if (expConvTest_.is_null()) {
01037     expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01038     expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
01039     expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), 
01040            Belos::TwoNorm);
01041   }
01042   // Convergence test first tests the implicit residual, then the
01043   // explicit residual if the implicit residual test passes.
01044   if (convTest_.is_null()) {
01045     convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ, 
01046               impConvTest_, 
01047               expConvTest_));
01048   }
01049   // Construct the complete iteration stopping criterion:
01050   //
01051   // "Stop iterating if the maximum number of iterations has been
01052   // reached, or if the convergence test passes."
01053   sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, 
01054                maxIterTest_, 
01055                convTest_));
01056   // Create the status test output class.
01057   // This class manages and formats the output from the status test.
01058   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
01059   outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_, 
01060            Passed+Failed+Undefined);
01061 
01062   // Set the solver string for the output test
01063   std::string solverDesc = " GCRODR ";
01064   outputTest_->setSolverDesc( solverDesc );
01065 
01066   // Create the timer if we need to.
01067   if (timerSolve_.is_null()) {
01068     std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
01069 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01070     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
01071 #endif
01072   }
01073 
01074   // Inform the solver manager that the current parameters were set.
01075   isSet_ = true;
01076 }
01077 
01078     
01079 template<class ScalarType, class MV, class OP>
01080 Teuchos::RCP<const Teuchos::ParameterList> 
01081 GCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const 
01082 {
01083   using Teuchos::ParameterList;
01084   using Teuchos::parameterList;
01085   using Teuchos::RCP;
01086 
01087   static RCP<const ParameterList> validPL;
01088   if (is_null(validPL)) {
01089     RCP<ParameterList> pl = parameterList ();
01090 
01091     // Set all the valid parameters and their default values.
01092     pl->set("Convergence Tolerance", convTol_default_,
01093       "The relative residual tolerance that needs to be achieved by the\n"
01094       "iterative solver in order for the linear system to be declared converged.");
01095     pl->set("Maximum Restarts", maxRestarts_default_,
01096       "The maximum number of cycles allowed for each\n"
01097       "set of RHS solved.");
01098     pl->set("Maximum Iterations", maxIters_default_,
01099       "The maximum number of iterations allowed for each\n"
01100       "set of RHS solved.");
01101     // mfh 25 Oct 2010: "Block Size" must be 1 because GCRODR is
01102     // currently not a block method: i.e., it does not work on
01103     // multiple right-hand sides at once.
01104     pl->set("Block Size", blockSize_default_,
01105       "Block Size Parameter -- currently must be 1 for GCRODR");
01106     pl->set("Num Blocks", numBlocks_default_,
01107       "The maximum number of vectors allowed in the Krylov subspace\n"
01108       "for each set of RHS solved.");
01109     pl->set("Num Recycled Blocks", recycledBlocks_default_,
01110       "The maximum number of vectors in the recycled subspace." );
01111     pl->set("Verbosity", verbosity_default_,
01112       "What type(s) of solver information should be outputted\n"
01113       "to the output stream.");
01114     pl->set("Output Style", outputStyle_default_,
01115       "What style is used for the solver information outputted\n"
01116       "to the output stream.");
01117     pl->set("Output Frequency", outputFreq_default_,
01118       "How often convergence information should be outputted\n"
01119       "to the output stream.");  
01120     pl->set("Output Stream", outputStream_default_,
01121       "A reference-counted pointer to the output stream where all\n"
01122       "solver output is sent.");
01123     pl->set("Implicit Residual Scaling", impResScale_default_,
01124       "The type of scaling used in the implicit residual convergence test.");
01125     pl->set("Explicit Residual Scaling", expResScale_default_,
01126       "The type of scaling used in the explicit residual convergence test.");
01127     pl->set("Timer Label", label_default_,
01128       "The string to use as a prefix for the timer labels.");
01129     //  pl->set("Restart Timers", restartTimers_);
01130     {
01131       OrthoManagerFactory<ScalarType, MV, OP> factory;
01132       pl->set("Orthogonalization", orthoType_default_,
01133         "The type of orthogonalization to use.  Valid options: " + 
01134         factory.validNamesString());
01135       RCP<const ParameterList> orthoParams = 
01136   factory.getDefaultParameters (orthoType_default_);
01137       pl->set ("Orthogonalization Parameters", *orthoParams, 
01138          "Parameters specific to the type of orthogonalization used.");
01139     }
01140     pl->set("Orthogonalization Constant", orthoKappa_default_,
01141       "When using DGKS orthogonalization: the \"depTol\" constant, used "
01142       "to determine whether another step of classical Gram-Schmidt is "
01143       "necessary.  Otherwise ignored.");
01144     validPL = pl;
01145   }
01146   return validPL;
01147 }
01148 
01149 // initializeStateStorage
01150 template<class ScalarType, class MV, class OP>
01151 void GCRODRSolMgr<ScalarType,MV,OP>::initializeStateStorage() {
01152 
01153     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01154 
01155     // Check if there is any multivector to clone from.
01156     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
01157     if (rhsMV == Teuchos::null) {
01158       // Nothing to do
01159       return;
01160     }
01161     else {
01162 
01163       // Initialize the state storage
01164       TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
01165                          "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
01166 
01167       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01168       if (U_ == Teuchos::null) {
01169         U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01170       }
01171       else {
01172         // Generate U_ by cloning itself ONLY if more space is needed.
01173         if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
01174           Teuchos::RCP<const MV> tmp = U_;
01175           U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01176         }
01177       }
01178 
01179       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01180       if (C_ == Teuchos::null) {
01181         C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01182       }
01183       else {
01184         // Generate C_ by cloning itself ONLY if more space is needed.
01185         if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
01186           Teuchos::RCP<const MV> tmp = C_;
01187           C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01188         }
01189       }
01190 
01191       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01192       if (V_ == Teuchos::null) {
01193         V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
01194       }
01195       else {
01196         // Generate V_ by cloning itself ONLY if more space is needed.
01197         if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
01198           Teuchos::RCP<const MV> tmp = V_;
01199           V_ = MVT::Clone( *tmp, numBlocks_+1 );
01200         }
01201       }
01202 
01203       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01204       if (U1_ == Teuchos::null) {
01205         U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01206       }
01207       else {
01208         // Generate U1_ by cloning itself ONLY if more space is needed.
01209         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01210           Teuchos::RCP<const MV> tmp = U1_;
01211           U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01212         }
01213       }
01214 
01215       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01216       if (C1_ == Teuchos::null) {
01217         C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01218       }
01219       else {
01220         // Generate C1_ by cloning itself ONLY if more space is needed.
01221         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01222           Teuchos::RCP<const MV> tmp = C1_;
01223           C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01224         }
01225       }
01226 
01227       // Generate r_ only if it doesn't exist
01228       if (r_ == Teuchos::null)
01229         r_ = MVT::Clone( *rhsMV, 1 );
01230 
01231       // Size of tau_ will change during computation, so just be sure it starts with appropriate size
01232       tau_.resize(recycledBlocks_+1);
01233 
01234       // Size of work_ will change during computation, so just be sure it starts with appropriate size
01235       work_.resize(recycledBlocks_+1);
01236 
01237       // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
01238       ipiv_.resize(recycledBlocks_+1);
01239 
01240       // Generate H2_ only if it doesn't exist, otherwise resize it.
01241       if (H2_ == Teuchos::null)
01242         H2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
01243       else {
01244         if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
01245           H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
01246       }
01247       H2_->putScalar(zero);
01248 
01249       // Generate R_ only if it doesn't exist, otherwise resize it.
01250       if (R_ == Teuchos::null)
01251         R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
01252       else {
01253         if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
01254           R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
01255       }
01256       R_->putScalar(zero);
01257 
01258       // Generate PP_ only if it doesn't exist, otherwise resize it.
01259       if (PP_ == Teuchos::null)
01260         PP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
01261       else {
01262         if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
01263           PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
01264       }
01265 
01266       // Generate HP_ only if it doesn't exist, otherwise resize it.
01267       if (HP_ == Teuchos::null)
01268         HP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
01269       else {
01270         if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
01271           HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
01272       }
01273 
01274     } // end else
01275 }
01276 
01277   
01278 // solve()
01279 template<class ScalarType, class MV, class OP>
01280 ReturnType GCRODRSolMgr<ScalarType,MV,OP>::solve() {
01281   using Teuchos::RCP;
01282   using Teuchos::rcp;
01283 
01284   // Set the current parameters if they were not set before.
01285   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
01286   // then didn't set any parameters using setParameters().
01287   if (!isSet_) { setParameters( params_ ); }
01288 
01289   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01290   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01291   std::vector<int> index(numBlocks_+1);
01292   
01293   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
01294 
01295   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
01296 
01297   // Create indices for the linear systems to be solved.
01298   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01299   std::vector<int> currIdx(1);
01300   currIdx[0] = 0;
01301 
01302   // Inform the linear problem of the current linear system to solve.
01303   problem_->setLSIndex( currIdx );
01304 
01305   // Check the number of blocks and change them is necessary.
01306   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
01307   if (numBlocks_ > dim) {
01308     numBlocks_ = dim;
01309     printer_->stream(Warnings) << 
01310       "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
01311       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
01312     params_->set("Num Blocks", numBlocks_);
01313   } 
01314 
01315   // Assume convergence is achieved, then let any failed convergence set this to false.
01316   bool isConverged = true;  
01317 
01318   // Initialize storage for all state variables
01319   initializeStateStorage();
01320 
01322   // Parameter list
01323   Teuchos::ParameterList plist;
01324   
01325   plist.set("Num Blocks",numBlocks_);
01326   plist.set("Recycled Blocks",recycledBlocks_);
01327   
01329   // GCRODR solver
01330   
01331   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
01332   gcrodr_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01333   // Number of iterations required to generate initial recycle space (if needed)
01334   int prime_iterations = 0;
01335 
01336   // Enter solve() iterations
01337   {
01338 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01339     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01340 #endif
01341     
01342     while ( numRHS2Solve > 0 ) {
01343 
01344       // Set flag indicating recycle space has not been generated this solve
01345       builtRecycleSpace_ = false;
01346       
01347       // Reset the status test.  
01348       outputTest_->reset();
01349       
01351       // Initialize recycled subspace for GCRODR
01352       
01353       // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
01354       if (keff > 0) {
01355   TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
01356          "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
01357 
01358   printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
01359   // Compute image of U_ under the new operator
01360   index.resize(keff);
01361   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01362   RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01363   RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01364   problem_->apply( *Utmp, *Ctmp );
01365 
01366   RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01367 
01368   // Orthogonalize this block
01369   // Get a matrix to hold the orthonormalization coefficients.
01370         Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff ); 
01371   int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,false));
01372   // Throw an error if we could not orthogonalize this block
01373   TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
01374   
01375   // U_ = U_*R^{-1} 
01376   // First, compute LU factorization of R
01377   int info = 0;
01378         ipiv_.resize(Rtmp.numRows());
01379   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01380   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01381   
01382   // Now, form inv(R)
01383   int lwork = Rtmp.numRows();
01384         work_.resize(lwork);
01385   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01386   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01387   
01388         // U_ = U1_; (via a swap)
01389   MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
01390         std::swap(U_, U1_);
01391 
01392         // Must reinitialize after swap
01393   index.resize(keff);
01394   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01395   Ctmp  = MVT::CloneViewNonConst( *C_, index );
01396   Utmp  = MVT::CloneView( *U_, index );
01397 
01398   // Compute C_'*r_
01399   Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
01400   problem_->computeCurrPrecResVec( &*r_ );
01401   MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
01402 
01403   // Update solution ( x += U_*C_'*r_ )
01404         RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
01405         MVT::MvInit( *update, 0.0 );
01406         MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
01407         problem_->updateSolution( update, true );
01408   
01409   // Update residual norm ( r -= C_*C_'*r_ )  
01410   MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
01411 
01412         // We recycled space from previous call
01413         prime_iterations = 0;
01414 
01415       }
01416       else {
01417   
01418   // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
01419   printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
01420   
01421   Teuchos::ParameterList primeList;
01422   
01423   // Tell the block solver that the block size is one.
01424         primeList.set("Num Blocks",numBlocks_);
01425         primeList.set("Recycled Blocks",0);
01426   
01427   //  Create GCRODR iterator object to perform one cycle of GMRES.
01428   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
01429   gcrodr_prime_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
01430   
01431   // Create the first block in the current Krylov basis (residual).
01432   problem_->computeCurrPrecResVec( &*r_ );
01433         index.resize( 1 ); index[0] = 0;
01434   RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01435   MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01436 
01437   // Set the new state and initialize the solver.
01438   GCRODRIterState<ScalarType,MV> newstate;
01439         index.resize( numBlocks_+1 );
01440         for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01441         newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01442   newstate.U = Teuchos::null;
01443   newstate.C = Teuchos::null;
01444         newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
01445   newstate.B = Teuchos::null;
01446   newstate.curDim = 0;
01447   gcrodr_prime_iter->initialize(newstate);
01448   
01449   // Perform one cycle of GMRES
01450   bool primeConverged = false;
01451   try {
01452     gcrodr_prime_iter->iterate();
01453 
01454     // Check convergence first
01455     if ( convTest_->getStatus() == Passed ) {
01456       // we have convergence
01457       primeConverged = true;
01458     }
01459   }
01460   catch (const GCRODRIterOrthoFailure &e) {
01461     // Try to recover the most recent least-squares solution
01462     gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
01463     
01464     // Check to see if the most recent least-squares solution yielded convergence.
01465     sTest_->checkStatus( &*gcrodr_prime_iter );
01466     if (convTest_->getStatus() == Passed)
01467       primeConverged = true;
01468   }
01469   catch (const std::exception &e) {
01470     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01471            << gcrodr_prime_iter->getNumIters() << std::endl 
01472            << e.what() << std::endl;
01473     throw;
01474   }
01475         // Record number of iterations in generating initial recycle spacec
01476         prime_iterations = gcrodr_prime_iter->getNumIters();
01477            
01478   // Update the linear problem.
01479   RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
01480   problem_->updateSolution( update, true );
01481 
01482   // Get the state.
01483   newstate = gcrodr_prime_iter->getState();
01484   int p = newstate.curDim;
01485 
01486   // Compute harmonic Ritz vectors 
01487   // NOTE:  The storage for the harmonic Ritz vectors (PP) is made one column larger 
01488         //        just in case we split a complex conjugate pair.
01489   // NOTE:  Generate a recycled subspace only if we have enough vectors.  If we converged
01490   //        too early, move on to the next linear system and try to generate a subspace again.
01491   if (recycledBlocks_ < p+1) {
01492     int info = 0;
01493           RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01494           // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
01495     keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
01496           // Hereafter, only keff columns of PP are needed
01497           PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
01498           // Now get views into C, U, V
01499           index.resize(keff);
01500           for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01501           RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01502           RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01503           RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01504     index.resize(p);
01505           for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
01506     RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01507 
01508     // Form U (the subspace to recycle)
01509           // U = newstate.V(:,1:p) * PP;
01510           MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01511 
01512     // Form orthonormalized C and adjust U so that C = A*U
01513 
01514     // First, compute [Q, R] = qr(H*P);
01515 
01516           // Step #1: Form HP = H*P
01517     Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
01518     Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
01519     HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
01520 
01521     // Step #1.5: Perform workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01522           int lwork = -1;
01523           tau_.resize(keff);
01524           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01525     TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01526 
01527           // Step #2: Compute QR factorization of HP
01528           lwork = (int)work_[0];
01529           work_.resize(lwork);
01530           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01531     TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01532 
01533           // Step #3: Explicitly construct Q and R factors 
01534     // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01535           Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
01536           for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01537           lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01538     TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01539 
01540           // Now we have [Q,R] = qr(H*P)
01541 
01542     // Now compute C = V(:,1:p+1) * Q 
01543     index.resize( p+1 );
01544           for (int ii=0; ii < (p+1); ++ii) { index[ii] = ii; }
01545     Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
01546           MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01547         
01548     // Finally, compute U = U*R^{-1}. 
01549           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as 
01550           // backsolve capabilities don't exist in the Belos::MultiVec class
01551 
01552     // Step #1: First, compute LU factorization of R
01553     ipiv_.resize(Rtmp.numRows());
01554     lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01555     TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01556     
01557     // Step #2: Form inv(R)
01558     lwork = Rtmp.numRows();
01559     work_.resize(lwork);
01560     lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01561     TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01562     
01563           // Step #3: Let U = U * R^{-1}
01564     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01565 
01566     printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01567 
01568   }  // if (recycledBlocks_ < p+1)
01569 
01570         // Return to outer loop if the priming solve converged, set the next linear system.
01571   if (primeConverged) {
01572     // Inform the linear problem that we are finished with this block linear system.
01573     problem_->setCurrLS();
01574 
01575           // Update indices for the linear systems to be solved.
01576           numRHS2Solve -= 1;
01577           if ( numRHS2Solve > 0 ) {
01578             currIdx[0]++;
01579  
01580             // Set the next indices.
01581             problem_->setLSIndex( currIdx );
01582           }
01583           else {
01584             currIdx.resize( numRHS2Solve );
01585           }
01586 
01587     continue;
01588   }
01589       } // if (keff > 0) ...
01590       
01591       // Prepare for the Gmres iterations with the recycled subspace.
01592 
01593       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01594       gcrodr_iter->setSize( keff, numBlocks_ );
01595       
01596       // Reset the number of iterations.
01597       gcrodr_iter->resetNumIters(prime_iterations);
01598 
01599       // Reset the number of calls that the status test output knows about.
01600       outputTest_->resetNumCalls();
01601 
01602       // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
01603       problem_->computeCurrPrecResVec( &*r_ );
01604       index.resize( 1 ); index[0] = 0;
01605       RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01606       MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01607 
01608       // Set the new state and initialize the solver.
01609       GCRODRIterState<ScalarType,MV> newstate;
01610       index.resize( numBlocks_+1 );
01611       for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01612       newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01613       index.resize( keff );
01614       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01615       newstate.C  = MVT::CloneViewNonConst( *C_,  index );
01616       newstate.U  = MVT::CloneViewNonConst( *U_,  index );
01617       newstate.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_,         keff, numBlocks_,    0, keff ) );
01618       newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01619       newstate.curDim = 0;
01620       gcrodr_iter->initialize(newstate);
01621 
01622       // variables needed for inner loop
01623       int numRestarts = 0;
01624       while(1) {
01625   
01626   // tell gcrodr_iter to iterate
01627   try {
01628     gcrodr_iter->iterate();
01629     
01631     //
01632     // check convergence first
01633     //
01635     if ( convTest_->getStatus() == Passed ) {
01636       // we have convergence
01637       break;  // break from while(1){gcrodr_iter->iterate()}
01638     }
01640     //
01641     // check for maximum iterations
01642     //
01644     else if ( maxIterTest_->getStatus() == Passed ) {
01645       // we don't have convergence
01646       isConverged = false;
01647       break;  // break from while(1){gcrodr_iter->iterate()}
01648     }
01650     //
01651     // check for restarting, i.e. the subspace is full
01652     //
01654     else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
01655 
01656             // Update the recycled subspace even if we have hit the maximum number of restarts.
01657  
01658             // Update the linear problem.
01659             RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01660             problem_->updateSolution( update, true );
01661 
01662             buildRecycleSpace2(gcrodr_iter);
01663 
01664             printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01665 
01666             // NOTE:  If we have hit the maximum number of restarts then quit
01667             if ( numRestarts >= maxRestarts_ ) {
01668               isConverged = false;
01669               break; // break from while(1){gcrodr_iter->iterate()}
01670             }
01671             numRestarts++;
01672 
01673             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01674 
01675             // Create the restart vector (first block in the current Krylov basis)
01676             problem_->computeCurrPrecResVec( &*r_ );
01677             index.resize( 1 ); index[0] = 0;
01678             RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01679             MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01680 
01681             // Set the new state and initialize the solver.
01682             GCRODRIterState<ScalarType,MV> restartState;
01683             index.resize( numBlocks_+1 );
01684             for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01685             restartState.V  = MVT::CloneViewNonConst( *V_,  index );
01686             index.resize( keff );
01687             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01688             restartState.U  = MVT::CloneViewNonConst( *U_,  index );
01689             restartState.C  = MVT::CloneViewNonConst( *C_,  index );
01690             restartState.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
01691             restartState.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01692             restartState.curDim = 0;
01693             gcrodr_iter->initialize(restartState);
01694 
01695 
01696     } // end of restarting
01697     
01699     //
01700     // we returned from iterate(), but none of our status tests Passed.
01701     // something is wrong, and it is probably our fault.
01702     //
01704     
01705     else {
01706       TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::GCRODRSolMgr::solve(): Invalid return from GCRODRIter::iterate().");
01707     }
01708   }
01709         catch (const GCRODRIterOrthoFailure &e) {
01710     // Try to recover the most recent least-squares solution
01711     gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
01712     
01713     // Check to see if the most recent least-squares solution yielded convergence.
01714     sTest_->checkStatus( &*gcrodr_iter );
01715     if (convTest_->getStatus() != Passed)
01716       isConverged = false;
01717     break;
01718         }
01719         catch (const std::exception &e) {
01720     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01721                              << gcrodr_iter->getNumIters() << std::endl 
01722            << e.what() << std::endl;
01723           throw;
01724   }
01725       }
01726      
01727       // Compute the current solution.
01728       // Update the linear problem.
01729       RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01730       problem_->updateSolution( update, true );
01731 
01732       // Inform the linear problem that we are finished with this block linear system.
01733       problem_->setCurrLS();
01734 
01735       // Update indices for the linear systems to be solved.
01736       numRHS2Solve -= 1;
01737       if ( numRHS2Solve > 0 ) {
01738   currIdx[0]++;
01739 
01740         // Set the next indices.
01741         problem_->setLSIndex( currIdx );
01742       }
01743       else {
01744         currIdx.resize( numRHS2Solve );
01745       }
01746 
01747       // If we didn't build a recycle space this solve but ran at least k iterations,
01748       // force build of new recycle space
01749 
01750       if (!builtRecycleSpace_) {
01751         buildRecycleSpace2(gcrodr_iter);
01752         printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01753       }
01754 
01755     }// while ( numRHS2Solve > 0 )
01756     
01757   }
01758   
01759   // print final summary
01760   sTest_->print( printer_->stream(FinalSummary) );
01761   
01762   // print timing information
01763 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01764   // Calling summarize() can be expensive, so don't call unless the
01765   // user wants to print out timing details.  summarize() will do all
01766   // the work even if it's passed a "black hole" output stream.
01767   if (verbosity_ & TimingDetails)
01768     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01769 #endif
01770  
01771   // get iteration information for this solve
01772   numIters_ = maxIterTest_->getNumIters();
01773 
01774   // Save the convergence test value ("achieved tolerance") for this
01775   // solve.  This solver (unlike BlockGmresSolMgr) always has two
01776   // residual norm status tests: an explicit and an implicit test.
01777   // The master convergence test convTest_ is a SEQ combo of the
01778   // implicit resp. explicit tests.  If the implicit test never
01779   // passes, then the explicit test won't ever be executed.  This
01780   // manifests as expConvTest_->getTestValue()->size() < 1.  We deal
01781   // with this case by using the values returned by
01782   // impConvTest_->getTestValue().
01783   {
01784     const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
01785     if (pTestValues == NULL || pTestValues->size() < 1) {
01786       pTestValues = impConvTest_->getTestValue();
01787     }
01788     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
01789       "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
01790       "method returned NULL.  Please report this bug to the Belos developers.");
01791     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
01792       "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
01793       "method returned a vector of length zero.  Please report this bug to the "
01794       "Belos developers.");
01795 
01796     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01797     // achieved tolerances for all vectors in the current solve(), or
01798     // just for the vectors from the last deflation?
01799     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01800   }
01801  
01802   if (!isConverged) {
01803     return Unconverged; // return from GCRODRSolMgr::solve() 
01804   }
01805   return Converged; // return from GCRODRSolMgr::solve() 
01806 }
01807 
01808 //  Given existing recycle space and Krylov space, build new recycle space
01809 template<class ScalarType, class MV, class OP>
01810 void GCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter) {
01811 
01812   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01813   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01814 
01815   std::vector<MagnitudeType> d(keff);
01816   std::vector<int> index(numBlocks_+1);
01817 
01818   // Get the state
01819   GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
01820   int p = oldState.curDim;
01821 
01822   // insufficient new information to update recycle space
01823   if (p<1) return;
01824 
01825   // Take the norm of the recycled vectors
01826   {
01827     index.resize(keff);
01828     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01829     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01830     d.resize(keff);
01831     MVT::MvNorm( *Utmp, d );
01832     for (int i=0; i<keff; ++i) {
01833       d[i] = one / d[i];
01834     }
01835     MVT::MvScale( *Utmp, d );
01836   }
01837 
01838   // Get view into current "full" upper Hessnburg matrix
01839   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
01840 
01841   // Insert D into the leading keff x keff  block of H2
01842   for (int i=0; i<keff; ++i) {
01843     (*H2tmp)(i,i) = d[i];
01844   }
01845 
01846   // Compute the harmoic Ritz pairs for the generalized eigenproblem
01847   // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
01848   int keff_new;
01849   {
01850     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
01851     keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, PPtmp );
01852   }
01853 
01854   // Code to form new U, C
01855   // U = [U V(:,1:p)] * P; (in two steps)
01856 
01857   // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01858   Teuchos::RCP<MV> U1tmp;
01859   {
01860     index.resize( keff );
01861     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01862     Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
01863     index.resize( keff_new );
01864     for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01865     U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
01866     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
01867     MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
01868   }
01869 
01870   // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01871   {
01872     index.resize(p);
01873     for (int ii=0; ii < p; ii++) { index[ii] = ii; }
01874     Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01875     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
01876     MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
01877   }
01878 
01879   // Form HP = H*P
01880   Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
01881   {
01882     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
01883     HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
01884   }
01885 
01886   // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01887   int info = 0, lwork = -1;
01888   tau_.resize(keff_new);
01889   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01890   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01891 
01892   lwork = (int)work_[0];
01893   work_.resize(lwork);
01894   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01895   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01896 
01897   // Explicitly construct Q and R factors
01898   // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01899   Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
01900   for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
01901   lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01902   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01903 
01904   // Form orthonormalized C and adjust U accordingly so that C = A*U
01905   // C = [C V] * Q;
01906 
01907   // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01908   {
01909     Teuchos::RCP<MV> C1tmp;
01910     {
01911       index.resize(keff);
01912       for (int i=0; i < keff; i++) { index[i] = i; }
01913       Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_,  index );
01914       index.resize(keff_new);
01915       for (int i=0; i < keff_new; i++) { index[i] = i; }
01916       C1tmp  = MVT::CloneViewNonConst( *C1_,  index );
01917       Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
01918       MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
01919     }
01920     // Now compute C += V(:,1:p+1) * Q
01921     {
01922       index.resize( p+1 );
01923       for (int i=0; i < p+1; ++i) { index[i] = i; }
01924       Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01925       Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
01926       MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
01927     }
01928   }
01929 
01930   // C_ = C1_; (via a swap)
01931   std::swap(C_, C1_);
01932 
01933   // Finally, compute U_ = U_*R^{-1}
01934   // First, compute LU factorization of R
01935   ipiv_.resize(Rtmp.numRows());
01936   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01937   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01938 
01939   // Now, form inv(R)
01940   lwork = Rtmp.numRows();
01941   work_.resize(lwork);
01942   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01943   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01944 
01945   {
01946     index.resize(keff_new);
01947     for (int i=0; i < keff_new; i++) { index[i] = i; }
01948     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_,  index );
01949     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01950   }
01951 
01952   // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01953   if (keff != keff_new) {
01954     keff = keff_new;
01955     gcrodr_iter->setSize( keff, numBlocks_ );
01956     // Important to zero this out before next cyle
01957     Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01958     b1.putScalar(zero);
01959   }
01960 
01961 }
01962 
01963 
01964 //  Compute the harmonic eigenpairs of the projected, dense system.
01965 template<class ScalarType, class MV, class OP>
01966 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs1(int m, 
01967                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01968                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01969   int i, j;
01970   bool xtraVec = false;
01971   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01972   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01973 
01974   // Real and imaginary eigenvalue components
01975   std::vector<MagnitudeType> wr(m), wi(m);
01976 
01977   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01978   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,false);
01979 
01980   // Magnitude of harmonic Ritz values
01981   std::vector<MagnitudeType> w(m);
01982 
01983   // Sorted order of harmonic Ritz values, also used for DGEEV
01984   std::vector<int> iperm(m);
01985 
01986   // Size of workspace and workspace for DGEEV
01987   int lwork = 4*m;
01988   std::vector<ScalarType> work(lwork);
01989 
01990   // Output info
01991   int info = 0;
01992 
01993   // Set flag indicating recycle space has been generated this solve
01994   builtRecycleSpace_ = true;
01995 
01996   // Solve linear system:  H_m^{-H}*e_m 
01997   Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
01998   Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
01999   e_m[m-1] = one;
02000   lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
02001   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
02002 
02003   // Compute H_m + d*H_m^{-H}*e_m*e_m^H
02004   ScalarType d = HH(m, m-1) * HH(m, m-1);
02005   Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( HH );
02006   for( i=0; i<m; ++i )
02007     harmHH(i, m-1) += d * e_m[i];
02008 
02009   // Revise to do query for optimal workspace first
02010   // Create simple storage for the left eigenvectors, which we don't care about.
02011   const int ldvl = m;
02012   ScalarType* vl = 0;
02013   lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
02014         vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
02015   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
02016 
02017   // Construct magnitude of each harmonic Ritz value
02018   for( i=0; i<m; ++i )
02019     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
02020 
02021   // Construct magnitude of each harmonic Ritz value
02022   this->sort(w, m, iperm);
02023 
02024   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
02025   if (wi[iperm[recycledBlocks_-1]] != zero) {
02026     int countImag = 0;
02027     for ( i=0; i<recycledBlocks_; ++i ) {
02028       if (wi[iperm[i]] != zero)
02029   countImag++;
02030     }
02031     // Check to see if this count is even or odd:
02032     if (countImag % 2)
02033       xtraVec = true;
02034   }
02035 
02036   // Select recycledBlocks_ smallest eigenvectors
02037   for( i=0; i<recycledBlocks_; ++i ) {
02038     for( j=0; j<m; j++ ) {
02039       PP(j,i) = vr(j,iperm[i]);
02040     }
02041   }
02042   
02043   if (xtraVec) { // we need to store one more vector
02044     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
02045       for( j=0; j<m; ++j ) {   // so get the "imag" component
02046   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
02047       }
02048     }
02049     else { //  I picked the "imag" component
02050       for( j=0; j<m; ++j ) {   // so get the "real" component
02051   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
02052       }
02053     }
02054   }
02055 
02056   // Return whether we needed to store an additional vector
02057   if (xtraVec) {
02058     return recycledBlocks_+1;
02059   }
02060   return recycledBlocks_;
02061 }
02062 
02063 //  Compute the harmonic eigenpairs of the projected, dense system.
02064 template<class ScalarType, class MV, class OP>
02065 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs2(int keff, int m, 
02066                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
02067                  const Teuchos::RCP<const MV>& VV,
02068                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
02069   int i, j;
02070   int m2 = HH.numCols();
02071   bool xtraVec = false;
02072   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
02073   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
02074   std::vector<int> index;
02075   
02076   // Real and imaginary eigenvalue components
02077   std::vector<ScalarType> wr(m2), wi(m2);
02078 
02079   // Magnitude of harmonic Ritz values
02080   std::vector<MagnitudeType> w(m2);
02081 
02082   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
02083   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,false);
02084 
02085   // Sorted order of harmonic Ritz values
02086   std::vector<int> iperm(m2);
02087 
02088   // Set flag indicating recycle space has been generated this solve
02089   builtRecycleSpace_ = true;
02090 
02091   // Form matrices for generalized eigenproblem
02092   
02093   // B = H2' * H2; Don't zero out matrix when constructing
02094   Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,false);
02095   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
02096   
02097   // A_tmp = | C'*U        0 |
02098   //         | V_{m+1}'*U  I |
02099   Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keff+m+1, keff+m );
02100 
02101   // A_tmp(1:keff,1:keff) = C' * U;
02102   index.resize(keff);
02103   for (int i=0; i<keff; ++i) { index[i] = i; }
02104   Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_, index );
02105   Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
02106   Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keff, keff );
02107   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
02108 
02109   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
02110   Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keff, keff );
02111   index.resize(m+1);
02112   for (i=0; i < m+1; i++) { index[i] = i; }
02113   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
02114   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
02115 
02116   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
02117   for( i=keff; i<keff+m; i++ ) {
02118     A_tmp(i,i) = one;
02119   }
02120 
02121   // A = H2' * A_tmp;
02122   Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
02123   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
02124 
02125   // Compute k smallest harmonic Ritz pairs
02126   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
02127   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
02128   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
02129   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
02130   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
02131   char balanc='P', jobvl='N', jobvr='V', sense='N'; 
02132   int ld = A.numRows();
02133   int lwork = 6*ld;
02134   int ldvl = ld, ldvr = ld;
02135   int info = 0,ilo = 0,ihi = 0;
02136   ScalarType abnrm = zero, bbnrm = zero;
02137   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
02138   std::vector<ScalarType> beta(ld);
02139   std::vector<ScalarType> work(lwork);
02140   std::vector<MagnitudeType> lscale(ld), rscale(ld);  
02141   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
02142   std::vector<int> iwork(ld+6);
02143   int *bwork = 0; // If sense == 'N', bwork is never referenced
02144   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0], 
02145                &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0], 
02146                &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
02147   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
02148   
02149   // Construct magnitude of each harmonic Ritz value
02150   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
02151   for( i=0; i<ld; i++ ) {
02152     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
02153   }
02154 
02155   // Construct magnitude of each harmonic Ritz value
02156   this->sort(w,ld,iperm);
02157 
02158   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
02159   if (wi[iperm[ld-recycledBlocks_]] != zero) {
02160     int countImag = 0;
02161     for ( i=ld-recycledBlocks_; i<ld; i++ ) {
02162       if (wi[iperm[i]] != zero)
02163   countImag++;
02164     }
02165     // Check to see if this count is even or odd:
02166     if (countImag % 2)
02167       xtraVec = true;
02168   }
02169   
02170   // Select recycledBlocks_ smallest eigenvectors
02171   for( i=0; i<recycledBlocks_; i++ ) {
02172     for( j=0; j<ld; j++ ) {
02173       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
02174     }
02175   }
02176   
02177   if (xtraVec) { // we need to store one more vector
02178     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
02179       for( j=0; j<ld; j++ ) {   // so get the "imag" component
02180   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
02181       }
02182     }
02183     else { // I picked the "imag" component
02184       for( j=0; j<ld; j++ ) {   // so get the "real" component
02185   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
02186       }
02187     }
02188   }
02189   
02190   // Return whether we needed to store an additional vector
02191   if (xtraVec) {
02192     return recycledBlocks_+1;
02193   }
02194   return recycledBlocks_;
02195 }
02196 
02197 
02198 // This method sorts list of n floating-point numbers and return permutation vector
02199 template<class ScalarType, class MV, class OP>
02200 void GCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
02201   int l, r, j, i, flag;
02202   int    RR2;
02203   double dRR, dK;
02204   
02205   // Initialize the permutation vector.
02206   for(j=0;j<n;j++)
02207     iperm[j] = j;
02208   
02209   if (n <= 1) return;
02210   
02211   l    = n / 2 + 1;
02212   r    = n - 1;
02213   l    = l - 1;
02214   dRR  = dlist[l - 1];
02215   dK   = dlist[l - 1];
02216   
02217   RR2 = iperm[l - 1];
02218   while (r != 0) {
02219     j = l;
02220     flag = 1;
02221     
02222     while (flag == 1) {
02223       i = j;
02224       j = j + j;
02225       
02226       if (j > r + 1)
02227   flag = 0;
02228       else {
02229   if (j < r + 1)
02230     if (dlist[j] > dlist[j - 1]) j = j + 1;
02231   
02232   if (dlist[j - 1] > dK) {
02233     dlist[i - 1] = dlist[j - 1];
02234     iperm[i - 1] = iperm[j - 1];
02235   }
02236   else {
02237     flag = 0;
02238   }
02239       }
02240     }
02241     dlist[i - 1] = dRR;
02242     iperm[i - 1] = RR2;
02243     
02244     if (l == 1) {
02245       dRR  = dlist [r];
02246       RR2 = iperm[r];
02247       dK = dlist[r];
02248       dlist[r] = dlist[0];
02249       iperm[r] = iperm[0];
02250       r = r - 1;
02251     }
02252     else {
02253       l   = l - 1;
02254       dRR  = dlist[l - 1];
02255       RR2  = iperm[l - 1];
02256       dK   = dlist[l - 1];
02257     }
02258   }
02259   dlist[0] = dRR;
02260   iperm[0] = RR2; 
02261 }
02262 
02263 
02264 //  This method requires the solver manager to return a string that describes itself.
02265 template<class ScalarType, class MV, class OP>
02266 std::string GCRODRSolMgr<ScalarType,MV,OP>::description() const {
02267   std::ostringstream oss;
02268   oss << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
02269   oss << "{";
02270   oss << "Ortho Type='"<<orthoType_;
02271   oss << ", Num Blocks=" <<numBlocks_;
02272   oss << ", Num Recycle Blocks=" << recycledBlocks_;
02273   oss << ", Max Restarts=" << maxRestarts_;
02274   oss << "}";
02275   return oss.str();
02276 }
02277 
02278 
02279 } // end Belos namespace
02280 
02281 #endif /* BELOS_GCRODR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines