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