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