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