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