Belos Version of the Day
BelosBlockGCRODRSolMgr.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 
00045 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
00046 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP
00047 
00048 #include "BelosConfigDefs.hpp"
00049 #include "BelosTypes.hpp"
00050 #include "BelosOrthoManagerFactory.hpp"
00051 #include "BelosLinearProblem.hpp"
00052 #include "BelosSolverManager.hpp"
00053 #include "BelosGmresIteration.hpp"
00054 #include "BelosBlockGCRODRIter.hpp"
00055 #include "BelosBlockGmresIter.hpp"
00056 #include "BelosBlockFGmresIter.hpp"
00057 #include "BelosStatusTestMaxIters.hpp"
00058 #include "BelosStatusTestGenResNorm.hpp"
00059 #include "BelosStatusTestCombo.hpp"
00060 #include "BelosStatusTestOutputFactory.hpp"
00061 #include "BelosOutputManager.hpp"
00062 #include "Teuchos_BLAS.hpp"
00063 #include "Teuchos_LAPACK.hpp"
00064 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00065 #include "Teuchos_TimeMonitor.hpp"
00066 #endif // BELOS_TEUCHOS_TIME_MONITOR
00067 
00068 // MLP: Add links to examples here
00069 
00070 namespace Belos{
00071 
00073 
00074 
00079   class BlockGCRODRSolMgrLinearProblemFailure : public BelosError {
00080   public:
00081     BlockGCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
00082   };
00083 
00088   class BlockGCRODRSolMgrOrthoFailure : public BelosError {
00089   public:
00090     BlockGCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00091   };
00092 
00097   class BlockGCRODRSolMgrLAPACKFailure : public BelosError {
00098   public:
00099     BlockGCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
00100   };
00101 
00106   class BlockGCRODRSolMgrRecyclingFailure : public BelosError {
00107   public:
00108     BlockGCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
00109   };
00110 
00112 
00124 
00125 template<class ScalarType, class MV, class OP>
00126 class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP> {
00127 private:
00128 
00129   typedef MultiVecTraits<ScalarType,MV> MVT;
00130   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00131   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00132   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00133   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00134   typedef Teuchos::ScalarTraits<MagnitudeType> SMT;
00135   typedef OrthoManagerFactory<ScalarType, MV, OP> ortho_factory_type;
00136   typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
00137   typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
00138 
00139 public:
00141 
00142 
00148   BlockGCRODRSolMgr();
00149     
00175   BlockGCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00176          const Teuchos::RCP<Teuchos::ParameterList> &pl);
00177 
00179   virtual ~BlockGCRODRSolMgr() {};
00181   
00184 
00186   std::string description() const;      
00187 
00189 
00190 
00192 
00193     
00195   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00196     return *problem_;
00197   }
00198 
00200   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00201 
00203   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00204 
00206   int getNumIters() const { return numIters_; }
00207 
00209   bool isLOADetected() const { return loaDetected_; }
00210 
00212 
00214 
00215 
00217   void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) {
00218     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00219       "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
00220 
00221     // Check state of problem object before proceeding
00222     if (! problem->isProblemSet()) {
00223       const bool success = problem->setProblem();
00224       TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
00225         "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed.  This likely means that the "
00226         "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B.  Please set these items in the LinearProblem and try again.");
00227     }
00228 
00229     problem_ = problem;
00230   }
00231 
00233   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00234 
00236 
00238 
00239     
00246   void reset (const ResetType type) {
00247     if ((type & Belos::Problem) && ! problem_.is_null ())
00248       problem_->setProblem();
00249     else if (type & Belos::RecycleSubspace)
00250       keff = 0;
00251   }
00252 
00254 
00256 
00257 
00273   // \returns ::ReturnType specifying:
00276   ReturnType solve();
00277 
00279 
00280 private: 
00281 
00282   // Called by all constructors; Contains init instructions common to all constructors
00283   void init();
00284 
00285   // Initialize solver state storage
00286   void initializeStateStorage();
00287 
00288   // Recycling Methods
00289   // Appending Function name by:
00290   //  "Kryl" indicates it is specialized for building a recycle space after an 
00291   //         initial run of Block GMRES which generates an initial unaugmented block Krylov subspace
00292   //  "AugKryl" indicates  it is specialized for building a recycle space from the augmented Krylov subspace
00293 
00294   // Functions which control the building of a recycle space
00295   void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter);
00296   void buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
00297 
00298   // Recycling with Harmonic Ritz Vectors
00299   // Computes harmonic eigenpairs of projected matrix created during the priming solve.
00300   // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00301 
00302   // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension (numBlocks+1)*blockSize x numBlocks.
00303   // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00304   int getHarmonicVecsKryl (int m, const SDM& HH, SDM& PP);
00305 
00306   // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+(numBlocks+1)*blockSize x keff+numBlocksm.
00307   // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) (numBlocks+1)*blockSize vectors.
00308   // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00309   int getHarmonicVecsAugKryl (int keff, 
00310             int m,
00311             const SDM& HH,
00312             const Teuchos::RCP<const MV>& VV,
00313             SDM& PP);
00314 
00315   // Sort list of n floating-point numbers and return permutation vector
00316   void sort (std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00317 
00318   // Lapack interface
00319   Teuchos::LAPACK<int,ScalarType> lapack;
00320 
00322   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00323 
00324   //Output Manager
00325   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00326   Teuchos::RCP<std::ostream> outputStream_;
00327 
00328   //Status Test
00329   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00330   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00331   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00332   Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00333   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00334 
00336   ortho_factory_type orthoFactory_;
00337 
00343   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00344 
00346   Teuchos::RCP<Teuchos::ParameterList> params_;
00347 
00358   mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00359 
00360   //Default Solver Values
00361   static const bool adaptiveBlockSize_default_;
00362   static const std::string recycleMethod_default_;
00363 
00364   //Current Solver Values
00365   MagnitudeType convTol_, orthoKappa_;
00366   int blockSize_, maxRestarts_, maxIters_, numIters_;
00367   int verbosity_, outputStyle_, outputFreq_;
00368   bool adaptiveBlockSize_;
00369   std::string orthoType_, recycleMethod_;
00370   std::string impResScale_, expResScale_;
00371   std::string label_;
00372 
00374   // Solver State Storage
00376   //
00377   // The number of blocks and recycle blocks (m and k, respectively)
00378   int numBlocks_, recycledBlocks_;
00379   // Current size of recycled subspace
00380   int keff;
00381   //
00382   // Residual Vector
00383   Teuchos::RCP<MV> R_;
00384   //
00385   // Search Space
00386   Teuchos::RCP<MV> V_;
00387   //
00388   // Recycle subspace and its image under action of the operator
00389   Teuchos::RCP<MV> U_, C_;
00390   //
00391   // Updated recycle Space and its image under action of the operator
00392   Teuchos::RCP<MV> U1_, C1_;
00393   //
00394   // Storage used in constructing recycle space
00395   Teuchos::RCP<SDM > G_;
00396   Teuchos::RCP<SDM > H_;
00397   Teuchos::RCP<SDM > B_;
00398   Teuchos::RCP<SDM > PP_;
00399   Teuchos::RCP<SDM > HP_;
00400   std::vector<ScalarType> tau_;
00401   std::vector<ScalarType> work_;
00402   Teuchos::RCP<SDM > F_;
00403   std::vector<int> ipiv_;
00404 
00406   Teuchos::RCP<Teuchos::Time> timerSolve_;
00407 
00409   bool isSet_;
00410 
00412   bool loaDetected_;
00413 
00415   bool builtRecycleSpace_; 
00416 };
00417 
00418   //
00419   // Set default solver values
00420   //
00421   template<class ScalarType, class MV, class OP>
00422   const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ = true;
00423 
00424   template<class ScalarType, class MV, class OP>
00425   const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ = "harmvecs";
00426 
00427   // 
00428   // Method definitions 
00429   //
00430 
00431   template<class ScalarType, class MV, class OP>
00432   BlockGCRODRSolMgr<ScalarType,MV,OP>::BlockGCRODRSolMgr() {
00433     init();
00434   }
00435 
00436   //Basic Constructor
00437   template<class ScalarType, class MV, class OP>
00438   BlockGCRODRSolMgr<ScalarType,MV,OP>::
00439   BlockGCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00440         const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
00441     // Initialize local pointers to null, and initialize local
00442     // variables to default values.
00443     init ();
00444 
00445     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00446       "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
00447       "the linear problem argument 'problem' to be nonnull.");
00448 
00449     problem_ = problem;
00450 
00451     // Set the parameters using the list that was passed in.  If null, we defer initialization until a non-null list is set (by the
00452     // client calling setParameters(), or by calling solve() -- in either case, a null parameter list indicates that default parameters should be used).
00453     if (! pl.is_null())
00454       setParameters (pl);
00455   }
00456 
00457   template<class ScalarType, class MV, class OP>
00458   void BlockGCRODRSolMgr<ScalarType,MV,OP>::init() {
00459     adaptiveBlockSize_ = adaptiveBlockSize_default_;
00460     recycleMethod_ = recycleMethod_default_;
00461     isSet_ = false;
00462     loaDetected_ = false;
00463     builtRecycleSpace_ = false;
00464     keff = 0;//Effective Size of Recycle Space
00465     //The following are all RCP smart pointers to indicated matrices/vectors.
00466     //Some MATLAB notation used in comments.
00467     R_ = Teuchos::null;//The Block Residual
00468     V_ = Teuchos::null;//Block Arnoldi Vectors
00469     U_ = Teuchos::null;//Recycle Space
00470     C_ = Teuchos::null;//Image of U Under Action of Operator
00471     U1_ = Teuchos::null;//Newly Computed Recycle Space
00472     C1_ = Teuchos::null;//Image of Newly Computed Recycle Space
00473     PP_ = Teuchos::null;//Coordinates of New Recyc. Vectors in Augmented Space
00474     HP_ = Teuchos::null;//H_*PP_ or G_*PP_
00475     G_ = Teuchos::null;//G_ such that A*[U V(:,1:numBlocks_*blockSize_)] = [C V_]*G_
00476     F_ = Teuchos::null;//Upper Triangular portion of QR factorization for HP_
00477     H_ = Teuchos::null;//H_ such that A*V(:,1:numBlocks_*blockSize_) = V_*H_ + C_*C_'*V_
00478     B_ = Teuchos::null;//B_ = C_'*V_
00479     
00480     //THIS BLOCK OF CODE IS COMMENTED OUT AND PLACED ELSEWHERE IN THE CODE
00481 /*
00482     //WE TEMPORARILY INITIALIZE STATUS TESTS HERE FOR TESTING PURPOSES, BUT THEY SHOULD BE 
00483     //INITIALIZED SOMEWHERE ELSE, LIKE THE setParameters() FUNCTION
00484 
00485     //INSTANTIATE AND INITIALIZE TEST OBJECTS AS NEEDED
00486     if (maxIterTest_.is_null()){
00487       maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00488     }
00489     //maxIterTest_->setMaxIters (maxIters_);
00490 
00491     //INSTANTIATE THE PRINTER
00492     if (printer_.is_null()) {
00493       printer_ = Teuchos::rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00494     }
00495 
00496     if (ortho_.is_null()) // || changedOrthoType) %%%In other codes, this is also triggered if orthogonalization type changed {
00497       // Create orthogonalization manager.  This requires that the OutputManager (printer_) already be initialized.
00498       Teuchos::RCP<const Teuchos::ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
00499       ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, orthoParams);
00500     }
00501 
00502     // Convenience typedefs
00503     typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00504     typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00505 
00506     if (impConvTest_.is_null()) {
00507       impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
00508       impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), Belos::TwoNorm);
00509       impConvTest_->setShowMaxResNormOnly( true );
00510     }
00511 
00512     if (expConvTest_.is_null()) {
00513       expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
00514       expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
00515       expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), Belos::TwoNorm);
00516       expConvTest_->setShowMaxResNormOnly( true );
00517     }
00518 
00519     if (convTest_.is_null()) {
00520       convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,  impConvTest_, expConvTest_));
00521     }
00522 
00523     sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
00524 
00525     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00526     outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
00527     Passed+Failed+Undefined);
00528 
00529  */
00530 
00531   }
00532 
00533   //  This method requires the solver manager to return a string that describes itself.
00534   template<class ScalarType, class MV, class OP>
00535   std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::description() const {
00536     std::ostringstream oss;
00537     oss << "Belos::BlockGCRODRSolMgr<" << SCT::name() << ", ...>";
00538     oss << "{";
00539     oss << "Ortho Type='"<<orthoType_ ;
00540     oss << ", Num Blocks=" <<numBlocks_;
00541     oss << ", Block Size=" <<blockSize_;
00542     oss << ", Num Recycle Blocks=" << recycledBlocks_;
00543     oss << ", Max Restarts=" << maxRestarts_;
00544     oss << "}";
00545     return oss.str();
00546   }
00547    
00548    template<class ScalarType, class MV, class OP>
00549    Teuchos::RCP<const Teuchos::ParameterList> 
00550    BlockGCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const {
00551      using Teuchos::ParameterList;
00552      using Teuchos::parameterList;
00553      using Teuchos::RCP;
00554      using Teuchos::rcpFromRef;
00555 
00556      if (defaultParams_.is_null()) {
00557        RCP<ParameterList> pl = parameterList ();
00558 
00559        const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
00560        const int maxRestarts = 1000;
00561        const int maxIters = 5000;
00562        const int blockSize = 2;
00563        const int numBlocks = 100;
00564        const int numRecycledBlocks = 25;
00565        const int verbosity = Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::StatusTestDetails;
00566        const Belos::OutputType outputStyle = Belos::General;
00567        const int outputFreq = 1;
00568        RCP<std::ostream> outputStream = rcpFromRef (std::cout);
00569        const std::string impResScale ("Norm of Preconditioned Initial Residual");
00570        const std::string expResScale ("Norm of Initial Residual");
00571        const std::string timerLabel ("Belos");
00572        const std::string orthoType ("DGKS");
00573        RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
00574        //const MagnitudeType orthoKappa = SCT::magnitude (SCT::eps());
00575        //
00576        // mfh 31 Dec 2011: Negative values of orthoKappa are ignored.
00577        // Setting this to a negative value by default ensures that
00578        // this parameter is only _not_ ignored if the user
00579        // specifically sets a valid value.
00580        const MagnitudeType orthoKappa = -SCT::one();
00581 
00582        // Set all the valid parameters and their default values.
00583        pl->set ("Convergence Tolerance", convTol,
00584     "The tolerance that the solver needs to achieve in order for "
00585     "the linear system(s) to be declared converged.  The meaning "
00586     "of this tolerance depends on the convergence test details.");
00587        pl->set("Maximum Restarts", maxRestarts,
00588          "The maximum number of restart cycles (not counting the first) "
00589          "allowed for each set of right-hand sides solved.");
00590        pl->set("Maximum Iterations", maxIters,
00591          "The maximum number of iterations allowed for each set of "
00592          "right-hand sides solved.");
00593        pl->set("Block Size", blockSize,
00594          "How many linear systems to solve at once.");
00595        pl->set("Num Blocks", numBlocks,
00596          "The maximum number of blocks allowed in the Krylov subspace "
00597          "for each set of right-hand sides solved.  This includes the "
00598          "initial block.  Total Krylov basis storage, not counting the "
00599          "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
00600        pl->set("Num Recycled Blocks", numRecycledBlocks,
00601          "The maximum number of vectors in the recycled subspace." );
00602        pl->set("Verbosity", verbosity,
00603          "What type(s) of solver information should be written "
00604          "to the output stream.");
00605        pl->set("Output Style", outputStyle,
00606          "What style is used for the solver information to write "
00607          "to the output stream.");
00608        pl->set("Output Frequency", outputFreq,
00609          "How often convergence information should be written "
00610          "to the output stream.");
00611        pl->set("Output Stream", outputStream,
00612          "A reference-counted pointer to the output stream where all "
00613          "solver output is sent.");
00614        pl->set("Implicit Residual Scaling", impResScale,
00615          "The type of scaling used in the implicit residual convergence test.");
00616        pl->set("Explicit Residual Scaling", expResScale,
00617          "The type of scaling used in the explicit residual convergence test.");
00618        pl->set("Timer Label", timerLabel,
00619          "The string to use as a prefix for the timer labels.");
00620        //  pl->set("Restart Timers", restartTimers_);
00621        pl->set("Orthogonalization", orthoType,
00622          "The orthogonalization method to use.  Valid options: " +
00623          orthoFactory_.validNamesString());
00624        pl->set ("Orthogonalization Parameters", *orthoParams,
00625     "Sublist of parameters specific to the orthogonalization method to use.");
00626        pl->set("Orthogonalization Constant", orthoKappa,
00627          "When using DGKS orthogonalization: the \"depTol\" constant, used "
00628          "to determine whether another step of classical Gram-Schmidt is "
00629          "necessary.  Otherwise ignored.  Nonpositive values are ignored.");
00630        defaultParams_ = pl;
00631      }
00632      return defaultParams_;
00633    }
00634    
00635    template<class ScalarType, class MV, class OP>
00636    void
00637    BlockGCRODRSolMgr<ScalarType,MV,OP>::
00638    setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params) {
00639      using Teuchos::isParameterType;
00640      using Teuchos::getParameter;
00641      using Teuchos::null;
00642      using Teuchos::ParameterList;
00643      using Teuchos::parameterList;
00644      using Teuchos::RCP;
00645      using Teuchos::rcp;
00646      using Teuchos::rcp_dynamic_cast;
00647      using Teuchos::rcpFromRef;
00648      using Teuchos::Exceptions::InvalidParameter;
00649      using Teuchos::Exceptions::InvalidParameterName;
00650      using Teuchos::Exceptions::InvalidParameterType;
00651 
00652      // The default parameter list contains all parameters that BlockGCRODRSolMgr understands, and none that it doesn't
00653      RCP<const ParameterList> defaultParams = getValidParameters();
00654 
00655      if (params.is_null()) {
00656        // Users really should supply a non-null input ParameterList,
00657        // so that we can fill it in.  However, if they really did give
00658        // us a null list, be nice and use default parameters.  In this
00659        // case, we don't have to do any validation.
00660        params_ = parameterList (*defaultParams);
00661      } 
00662      else {
00663        // Validate the user's parameter list and fill in any missing parameters with default values.
00664        //
00665        // Currently, validation is quite strict.  The following line
00666        // will throw an exception for misspelled or extra parameters.
00667        // However, it will _not_ throw an exception if there are
00668        // missing parameters: these will be filled in with their
00669        // default values.  There is additional discussion of other
00670        // validation strategies in the comments of this function for
00671        // Belos::GCRODRSolMgr.
00672        params->validateParametersAndSetDefaults (*defaultParams);
00673        // No side effects on the solver until after validation passes.
00674        params_ = params;
00675      }
00676 
00677      // 
00678      // Validate and set parameters relating to the Krylov subspace
00679      // dimensions and the number of blocks.  
00680      //
00681      // We try to read and validate all these parameters' values
00682      // before setting them in the solver.  This is because the
00683      // validity of each may depend on the others' values.  Our goal
00684      // is to avoid incorrect side effects, so we don't set the values
00685      // in the solver until we know they are right.
00686      //
00687      {
00688        const int maxRestarts = params_->get<int> ("Maximum Restarts");
00689        TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
00690          "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
00691          "must be nonnegative, but you specified a negative value of " 
00692          << maxRestarts << ".");
00693 
00694        const int maxIters = params_->get<int> ("Maximum Iterations");
00695        TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
00696          "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
00697          "must be positive, but you specified a nonpositive value of " 
00698          << maxIters << ".");
00699 
00700        const int numBlocks = params_->get<int> ("Num Blocks");
00701        TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
00702          "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
00703          "positive, but you specified a nonpositive value of " << numBlocks 
00704          << ".");
00705 
00706        const int blockSize = params_->get<int> ("Block Size");
00707        TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
00708          "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
00709          "positive, but you specified a nonpositive value of " << blockSize 
00710          << ".");
00711 
00712        const int recycledBlocks = params_->get<int> ("Num Recycled Blocks");
00713        TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
00714          "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
00715          "be positive, but you specified a nonpositive value of " 
00716          << recycledBlocks << ".");
00717        TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks, 
00718          std::invalid_argument, "Belos::BlockGCRODRSolMgr: The \"Num Recycled "
00719          "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
00720          "but you specified \"Num Recycled Blocks\" = " << recycledBlocks 
00721          << " and \"Num Blocks\" = " << numBlocks << ".");
00722   
00723        // Now that we've validated the various dimension-related
00724        // parameters, we can set them in the solvers.
00725        maxRestarts_ = maxRestarts;
00726        maxIters_ = maxIters;
00727        numBlocks_ = numBlocks;
00728        blockSize_ = blockSize;
00729        recycledBlocks_ = recycledBlocks;
00730      }
00731 
00732      // Check to see if the timer label changed.  If it did, update it in
00733      // the parameter list, and create a new timer with that label (if
00734      // Belos was compiled with timers enabled).
00735      {
00736        std::string tempLabel = params_->get<std::string> ("Timer Label");
00737        const bool labelChanged = (tempLabel != label_);
00738 
00739 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00740        std::string solveLabel = label_ + ": BlockGCRODRSolMgr total solve time";
00741        if (timerSolve_.is_null()) {
00742    // We haven't created a timer yet.
00743    timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00744        } else if (labelChanged) {
00745    // We created a timer before with a different label.  In that
00746    // case, clear the old timer and create a new timer with the
00747    // new label.  Clearing old timers prevents them from piling
00748    // up, since they are stored statically, possibly without
00749    // checking for duplicates.
00750    Teuchos::TimeMonitor::clearCounter (solveLabel);
00751    timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00752        }
00753 #endif // BELOS_TEUCHOS_TIME_MONITOR
00754      }
00755 
00756      // Check for a change in verbosity level.  Just in case, we allow
00757      // the type of this parameter to be either int or MsgType (an
00758      // enum).
00759      if (params_->isParameter ("Verbosity")) {
00760        if (isParameterType<int> (*params_, "Verbosity")) {
00761    verbosity_ = params_->get<int> ("Verbosity");
00762        } 
00763        else {
00764    verbosity_ = (int) getParameter<MsgType> (*params_, "Verbosity");
00765        }
00766      }
00767 
00768      // Check for a change in output style.  Just in case, we allow
00769      // the type of this parameter to be either int or OutputType (an
00770      // enum).
00771      if (params_->isParameter ("Output Style")) {
00772        if (isParameterType<int> (*params_, "Output Style")) {
00773    outputStyle_ = params_->get<int> ("Output Style");
00774        } 
00775        else {
00776    outputStyle_ = (int) getParameter<OutputType> (*params_, "Output Style");
00777        }
00778 
00779        // Currently, the only way we can update the output style is to
00780        // create a new output test.  This is because we must first
00781        // instantiate a new StatusTestOutputFactory with the new
00782        // style, and then use the factory to make a new output test.
00783        // Thus, we set outputTest_ to null for now, so we remember
00784        // later to create a new one.
00785        outputTest_ = null;
00786      }
00787   
00788      // Get the output stream for the output manager.
00789      //
00790      // It has been pointed out (mfh 28 Feb 2011 in GCRODRSolMgr code)
00791      // that including an RCP<std::ostream> parameter makes it
00792      // impossible to serialize the parameter list.  If you serialize
00793      // the list and read it back in from the serialized
00794      // representation, you might not even get a valid
00795      // RCP<std::ostream>, let alone the same output stream that you
00796      // serialized.
00797      //
00798      // However, existing Belos users depend on setting the output
00799      // stream in the parameter list.  We retain this behavior for
00800      // backwards compatibility.
00801      //
00802      // In case the output stream can't be read back in, we default to
00803      // stdout (std::cout), just to ensure reasonable behavior.
00804      if (params_->isParameter ("Output Stream")) {
00805        try {
00806    outputStream_ = getParameter<RCP<std::ostream> > (*params_, "Output Stream");
00807        } 
00808        catch (InvalidParameter&) {
00809    outputStream_ = rcpFromRef (std::cout);
00810        }
00811        // We assume that a null output stream indicates that the user
00812        // doesn't want to print anything, so we replace it with a
00813        // "black hole" stream that prints nothing sent to it.  (We
00814        // can't use outputStream_ = Teuchos::null, since the output
00815        // manager assumes that operator<< works on its output stream.)
00816        if (outputStream_.is_null()) {
00817    outputStream_ = rcp (new Teuchos::oblackholestream);
00818        }
00819      }
00820 
00821      outputFreq_ = params_->get<int> ("Output Frequency");
00822      if (verbosity_ & Belos::StatusTestDetails) {
00823        // Update parameter in our output status test.
00824        if (! outputTest_.is_null()) {
00825    outputTest_->setOutputFrequency (outputFreq_);
00826        }
00827      }
00828 
00829      // Create output manager if we need to, using the verbosity level
00830      // and output stream that we fetched above.  We do this here because
00831      // instantiating an OrthoManager using OrthoManagerFactory requires
00832      // a valid OutputManager.
00833      if (printer_.is_null()) {
00834        printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00835      } else {
00836        printer_->setVerbosity (verbosity_);
00837        printer_->setOStream (outputStream_);
00838      }
00839 
00840      // Get the orthogonalization manager name ("Orthogonalization").
00841      //
00842      // Getting default values for the orthogonalization manager
00843      // parameters ("Orthogonalization Parameters") requires knowing the
00844      // orthogonalization manager name.  Save it for later, and also
00845      // record whether it's different than before.
00846      bool changedOrthoType = false;
00847      if (params_->isParameter ("Orthogonalization")) {
00848        const std::string& tempOrthoType = 
00849    params_->get<std::string> ("Orthogonalization");
00850        // Ensure that the specified orthogonalization type is valid.
00851        if (! orthoFactory_.isValidName (tempOrthoType)) {
00852    std::ostringstream os;
00853    os << "Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
00854       << tempOrthoType << "\".  The following are valid options "
00855       << "for the \"Orthogonalization\" name parameter: ";
00856    orthoFactory_.printValidNames (os);
00857    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00858        }
00859        if (tempOrthoType != orthoType_) {
00860    changedOrthoType = true;
00861    orthoType_ = tempOrthoType;
00862        }
00863      }
00864 
00865      // Get any parameters for the orthogonalization
00866      // ("Orthogonalization Parameters").  If not supplied, the
00867      // orthogonalization manager factory will supply default values.
00868      //
00869      // NOTE (mfh 12 Jan 2011, summer 2011, 31 Dec 2011) For backwards
00870      // compatibility with other Belos GMRES-type solvers, if params
00871      // has an "Orthogonalization Constant" parameter and the DGKS
00872      // orthogonalization manager is to be used, the value of this
00873      // parameter will override DGKS's "depTol" parameter.
00874      //
00875      // Users must supply the orthogonalization manager parameters as
00876      // a sublist (supplying it as an RCP<ParameterList> would make
00877      // the resulting parameter list not serializable).
00878      RCP<ParameterList> orthoParams = sublist (params, "Orthogonalization Parameters", true);
00879      TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error, 
00880         "Failed to get orthogonalization parameters.  "
00881         "Please report this bug to the Belos developers.");
00882 
00883      // Check if the desired orthogonalization method changed, or if
00884      // the orthogonalization manager has not yet been instantiated.
00885      // If either is the case, instantiate a new MatOrthoManager
00886      // subclass instance corresponding to the desired
00887      // orthogonalization method.  We've already fetched the
00888      // orthogonalization method name (orthoType_) and its parameters
00889      // (orthoParams) above.
00890      //
00891      // As suggested (by mfh 12 Jan 2011 in Belos::GCRODRSolMgr) In
00892      // order to ensure parameter changes get propagated to the
00893      // orthomanager we simply reinstantiate the OrthoManager every
00894      // time, whether or not the orthogonalization method name or
00895      // parameters have changed.  This is not efficient for some
00896      // orthogonalization managers that keep a lot of internal state.
00897      // A more general way to fix this inefficiency would be to
00898      //
00899      // 1. make each orthogonalization implement 
00900      //    Teuchos::ParameterListAcceptor, and
00901      //
00902      // 2. make each orthogonalization implement a way to set the
00903      //    OutputManager instance and timer label.
00904      //
00905      // Create orthogonalization manager.  This requires that the
00906      // OutputManager (printer_) already be initialized.
00907      ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
00908              label_, orthoParams);
00909 
00910      // The DGKS orthogonalization accepts a "Orthogonalization
00911      // Constant" parameter (also called kappa in the code, but not in
00912      // the parameter list).  If its value is provided in the given
00913      // parameter list and its value is positive, use it.  Ignore
00914      // negative values.
00915      //
00916      // The "Orthogonalization Constant" parameter is retained only
00917      // for backwards compatibility.
00918      bool gotValidOrthoKappa = false;
00919      if (params_->isParameter ("Orthogonalization Constant")) {
00920        const MagnitudeType orthoKappa =
00921    params_->get<MagnitudeType> ("Orthogonalization Constant");
00922        if (orthoKappa > 0) {
00923    orthoKappa_ = orthoKappa;
00924    gotValidOrthoKappa = true;
00925    // Only DGKS currently accepts this parameter.
00926    if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
00927      typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
00928      // This cast should always succeed; it's a bug otherwise.
00929      // (If the cast fails, then orthoType_ doesn't correspond
00930      // to the OrthoManager subclass instance that we think we
00931      // have, so we initialized the wrong subclass somehow.)
00932      rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
00933    }
00934        }
00935      }
00936 
00937      // Convergence
00938      typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00939      typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00940 
00941      // Check for convergence tolerance
00942      convTol_ = params_->get<MagnitudeType> ("Convergence Tolerance");
00943      if (! impConvTest_.is_null()) {
00944        impConvTest_->setTolerance (convTol_);
00945      }
00946      if (! expConvTest_.is_null()) {
00947        expConvTest_->setTolerance (convTol_);
00948      }
00949 
00950      // Check for a change in scaling, if so we need to build new residual tests.
00951      if (params_->isParameter ("Implicit Residual Scaling")) {
00952        std::string tempImpResScale =
00953    getParameter<std::string> (*params_, "Implicit Residual Scaling");
00954 
00955        // Only update the scaling if it's different.
00956        if (impResScale_ != tempImpResScale) {
00957    ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
00958    impResScale_ = tempImpResScale;
00959 
00960    if (! impConvTest_.is_null()) {
00961      try {
00962        impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
00963      }
00964      catch (StatusTestError&) {
00965        // Delete the convergence test so it gets constructed again.
00966        impConvTest_ = null;
00967        convTest_ = null;
00968      }
00969    }
00970        }
00971      }
00972 
00973      if (params_->isParameter("Explicit Residual Scaling")) {
00974        std::string tempExpResScale =
00975    getParameter<std::string> (*params_, "Explicit Residual Scaling");
00976 
00977        // Only update the scaling if it's different.
00978        if (expResScale_ != tempExpResScale) {
00979    ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
00980    expResScale_ = tempExpResScale;
00981 
00982    if (! expConvTest_.is_null()) {
00983      try {
00984        expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
00985      }
00986      catch (StatusTestError&) {
00987        // Delete the convergence test so it gets constructed again.
00988        expConvTest_ = null;
00989        convTest_ = null;
00990      }
00991    }
00992        }
00993      }
00994 
00995      //
00996      // Create iteration stopping criteria ("status tests") if we need
00997      // to, by combining three different stopping criteria.
00998      //
00999      // First, construct maximum-number-of-iterations stopping criterion.
01000      if (maxIterTest_.is_null()) {
01001        maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
01002      } else {
01003        maxIterTest_->setMaxIters (maxIters_);
01004      }
01005 
01006      // Implicit residual test, using the native residual to determine if
01007      // convergence was achieved.
01008      if (impConvTest_.is_null()) {
01009        impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01010        impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
01011               Belos::TwoNorm);
01012      }
01013 
01014      // Explicit residual test once the native residual is below the tolerance
01015      if (expConvTest_.is_null()) {
01016        expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01017        expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
01018        expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
01019               Belos::TwoNorm);
01020      }
01021      // Convergence test first tests the implicit residual, then the
01022      // explicit residual if the implicit residual test passes.
01023      if (convTest_.is_null()) {
01024        convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
01025                  impConvTest_,
01026                  expConvTest_));
01027      }
01028      // Construct the complete iteration stopping criterion:
01029      //
01030      // "Stop iterating if the maximum number of iterations has been
01031      // reached, or if the convergence test passes."
01032      sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, 
01033             maxIterTest_, convTest_));
01034      // Create the status test output class.
01035      // This class manages and formats the output from the status test.
01036      StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
01037      outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
01038               Passed+Failed+Undefined);
01039 
01040      // Set the solver string for the output test
01041      std::string solverDesc = "Block GCRODR ";
01042      outputTest_->setSolverDesc (solverDesc);
01043 
01044      // Inform the solver manager that the current parameters were set.
01045      isSet_ = true;
01046    }
01047    
01048   // initializeStateStorage.
01049   template<class ScalarType, class MV, class OP>
01050   void 
01051   BlockGCRODRSolMgr<ScalarType,MV,OP>::initializeStateStorage()
01052   {
01053 
01054     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01055 
01056     // Check if there is any multivector to clone from.
01057     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
01058  
01059     //The Dimension of the Krylov Subspace 
01060     int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
01061 
01062     //Number of columns in [U_ V_(:,1:KrylSpaDim)]
01063     int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;// + 1 is for possible extra recycle vector
01064  
01065     //Number of columns in [C_ V_]
01066     int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
01067 
01068     //TEMPORARY SKELETON DEFINITION OF THIS FUNCTION TO GET THINGS WORKING
01069     //NOT EVERYTHING IS INITIALIZE CORRECTLY YET.
01070 
01071     //INITIALIZE RECYCLE SPACE VARIABLES HERE
01072 
01073     //WE DO NOT ALLOCATE V HERE IN THIS SOLVERMANAGER. If no recycle space exists, then block_gmres_iter
01074     //will allocated V for us.  If a recycle space already exists, then we will allocate V after updating the
01075     //recycle space for the current problem.
01076     // If the block Krylov subspace has not been initialized before, generate it using the RHS from lp_.
01077     /*if (V_ == Teuchos::null) {
01078       V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
01079       }
01080       else{
01081       // Generate V_ by cloning itself ONLY if more space is needed.
01082       if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
01083       Teuchos::RCP<const MV> tmp = V_;
01084       V_ = MVT::Clone( *tmp, numBlocks_+1 );
01085       }
01086       }*/
01087 
01088     //INTITIALIZE SPACE FOR THE NEWLY COMPUTED RECYCLE SPACE VARIABLES HERE
01089 
01090     if (U_ == Teuchos::null) {
01091       U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01092     }
01093     else {
01094       // Generate U_ by cloning itself ONLY if more space is needed.
01095       if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
01096   Teuchos::RCP<const MV> tmp = U_;
01097   U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01098       }
01099     }
01100 
01101     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01102     if (C_ == Teuchos::null) {
01103       C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01104     }
01105     else {
01106       // Generate C_ by cloning itself ONLY if more space is needed.
01107       if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
01108   Teuchos::RCP<const MV> tmp = C_;
01109   C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01110       }
01111     }
01112 
01113     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01114     if (U1_ == Teuchos::null) {
01115       U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01116     }
01117     else {
01118       // Generate U1_ by cloning itself ONLY if more space is needed.
01119       if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01120   Teuchos::RCP<const MV> tmp = U1_;
01121   U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01122       }
01123     }
01124 
01125     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01126     if (C1_ == Teuchos::null) {
01127       C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01128     }
01129     else {
01130       // Generate C1_ by cloning itself ONLY if more space is needed.
01131       if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01132   Teuchos::RCP<const MV> tmp = C1_;
01133   C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01134       }
01135     }
01136 
01137     // Generate R_ only if it doesn't exist
01138     if (R_ == Teuchos::null){
01139       R_ = MVT::Clone( *rhsMV, blockSize_ );
01140     }
01141 
01142     //INITIALIZE SOME WORK VARIABLES
01143         
01144     // Generate G_ only if it doesn't exist, otherwise resize it.
01145     if (G_ == Teuchos::null){
01146       G_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
01147     }
01148     else{
01149       if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )       
01150   {
01151     G_->reshape( augSpaImgDim, augSpaDim );
01152   }
01153       G_->putScalar(zero);
01154     }
01155 
01156     // Generate H_ only if it doesn't exist by pointing it to a view of G_.
01157     if (H_ == Teuchos::null){
01158       H_ = Teuchos::rcp (new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
01159     }
01160 
01161     // Generate F_ only if it doesn't exist, otherwise resize it.
01162     if (F_ == Teuchos::null){
01163       F_ = Teuchos::rcp( new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
01164     }
01165     else {
01166       if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
01167   F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
01168       }
01169     }
01170     F_->putScalar(zero);
01171 
01172     // Generate PP_ only if it doesn't exist, otherwise resize it.
01173     if (PP_ == Teuchos::null){
01174       PP_ = Teuchos::rcp( new SDM( augSpaImgDim, recycledBlocks_+1 ) );
01175     }
01176     else{
01177       if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
01178   PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
01179       }
01180     }
01181 
01182     // Generate HP_ only if it doesn't exist, otherwise resize it.
01183     if (HP_ == Teuchos::null)
01184       HP_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
01185     else{
01186       if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
01187   HP_->reshape( augSpaImgDim, augSpaDim );
01188       }
01189     }
01190 
01191     // Size of tau_ will change during computation, so just be sure it starts with appropriate size
01192     tau_.resize(recycledBlocks_+1);
01193 
01194     // Size of work_ will change during computation, so just be sure it starts with appropriate size
01195     work_.resize(recycledBlocks_+1);
01196 
01197     // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
01198     ipiv_.resize(recycledBlocks_+1);
01199 
01200   }
01201 
01202 template<class ScalarType, class MV, class OP>
01203 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
01204 
01205   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01206   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01207 
01208 
01209   int p = block_gmres_iter->getState().curDim; //Dimension of the Krylov space generated
01210   std::vector<int> index(keff);//we use this to index certain columns of U, C, and V to 
01211   //get views into pieces of these matrices.
01212 
01213   //GET CORRECT PIECE OF MATRIX H APPROPRIATE TO SIZE OF KRYLOV SUBSPACE
01214   SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
01215   if(recycledBlocks_ >= p + blockSize_){//keep whole block Krylov subspace
01216     //IF THIS HAS HAPPENED, THIS MEANS WE CONVERGED DURING THIS CYCLE
01217     //THEREFORE, WE DO NOT CONSTRUCT C = A*U;
01218     index.resize(p);
01219     for (int ii=0; ii < p; ++ii) index[ii] = ii;
01220     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01221     MVT::SetBlock(*V_, index, *Utmp);
01222     keff = p;
01223   }
01224   else{ //use a subspace selection method to get recycle space
01225     int info = 0; 
01226     Teuchos::RCP<SDM > PPtmp = rcp (new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01227     if(recycleMethod_ == "harmvecs"){
01228       keff = getHarmonicVecsKryl(p, HH, *PPtmp);
01229       printer_->stream(Debug) << "keff = " << keff << std::endl;
01230     }
01231 // Hereafter, only keff columns of PP are needed
01232 PPtmp = rcp (new SDM ( Teuchos::View, *PP_, p, keff ) );
01233 // Now get views into C, U, V
01234 index.resize(keff);
01235 for (int ii=0; ii<keff; ++ii) index[ii] = ii;
01236 Teuchos::RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01237 Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01238 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01239 index.resize(p);
01240 for (int ii=0; ii < p; ++ii) index[ii] = ii;
01241 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01242 
01243 // Form U (the subspace to recycle)
01244 // U = newstate.V(:,1:p) * PP;
01245 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01246 
01247 // Step #1: Form HP = H*P
01248 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
01249 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
01250 // Step #1.5: Perform workspace size query for QR factorization of HP
01251 int lwork = -1;
01252 tau_.resize(keff);
01253 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01254 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01255 
01256 // Step #2: Compute QR factorization of HP
01257 lwork = (int)work_[0];
01258 work_.resize(lwork);
01259 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01260 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,  "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01261 
01262 // Step #3: Explicitly construct Q and R factors
01263 // The upper triangular part of HP is copied into R and HP becomes Q.
01264 SDM Rtmp( Teuchos::View, *F_, keff, keff );
01265 for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01266 lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01267 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01268     // Now we have [Q,R] = qr(H*P)
01269 
01270             // Now compute C = V(:,1:p+blockSize_) * Q
01271                 index.resize( p+blockSize_ );
01272             for (int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
01273             Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+blockSize_ vectors now; needed p above)
01274             MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01275 
01276     // Finally, compute U = U*R^{-1}.
01277           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
01278                 // backsolve capabilities don't exist in the Belos::MultiVec class
01279 
01280 // Step #1: First, compute LU factorization of R
01281 ipiv_.resize(Rtmp.numRows());
01282 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01283 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01284 // Step #2: Form inv(R)
01285 lwork = Rtmp.numRows();
01286 work_.resize(lwork);
01287 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01288 //TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01289 // Step #3: Let U = U * R^{-1}
01290 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01291 
01292 } //end else from if(recycledBlocks_ >= p + 1)
01293 return;
01294 } // end buildRecycleSpaceKryl defnition
01295 
01296 template<class ScalarType, class MV, class OP>
01297 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
01298   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01299   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01300 
01301   std::vector<MagnitudeType> d(keff);
01302   std::vector<int> index(numBlocks_+1);
01303 
01304   // Get the state
01305   BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
01306   int p = oldState.curDim;
01307 
01308   // insufficient new information to update recycle space
01309   if (p<1) return;
01310 
01311   if(recycledBlocks_ >= keff + p){ //we add new Krylov vectors to existing recycle space
01312      // If this has happened, we have converged during this cycle. Therefore, do not construct C = A*C;
01313      index.resize(p);
01314      for (int ii=0; ii < p; ++ii) { index[ii] = keff+ii; } //get a view after current reycle vectors
01315      Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01316      for (int ii=0; ii < p; ++ii) { index[ii] =ii; }
01317      MVT::SetBlock(*V_, index, *Utmp);
01318      keff += p;
01319   }
01320 
01321   // Take the norm of the recycled vectors
01322   {
01323     index.resize(keff);
01324     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01325     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01326     d.resize(keff);
01327     MVT::MvNorm( *Utmp, d );
01328     for (int i=0; i<keff; ++i)
01329       d[i] = one / d[i];
01330     MVT::MvScale( *Utmp, d );
01331   }
01332 
01333   // Get view into current "full" upper Hessnburg matrix
01334   // note that p describes the dimension of the iter+1 block Krylov space so we have to adjust how we use it to index Gtmp
01335   Teuchos::RCP<SDM> Gtmp = Teuchos::rcp( new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
01336 
01337   // Insert D into the leading keff x keff  block of G 
01338   for (int i=0; i<keff; ++i)
01339     (*Gtmp)(i,i) = d[i];
01340 
01341   // Compute the harmoic Ritz pairs for the generalized eigenproblem
01342   // getHarmonicVecsKryl assumes PP has recycledBlocks_+1 columns available
01343   // See previous block of comments for why we subtract p-blockSize_
01344   int keff_new;
01345   {
01346     SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
01347     keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
01348   }
01349 
01350   // Code to form new U, C
01351   // U = [U V(:,1:p)] * P; (in two steps)
01352   
01353   // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01354   Teuchos::RCP<MV> U1tmp;
01355   {
01356     index.resize( keff );
01357     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01358     Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
01359     index.resize( keff_new );
01360     for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01361     U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
01362     SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
01363     MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
01364   }
01365 
01366   // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01367   {
01368     index.resize(p-blockSize_);
01369     for (int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
01370     Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01371     SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
01372     MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
01373   }
01374 
01375   // Form GP = G*P
01376   SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
01377   {
01378     SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
01379     HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
01380   }
01381 
01382   // Workspace size query for QR factorization HP= QF (the worksize will be placed in work_[0])
01383   int info = 0, lwork = -1;
01384   tau_.resize(keff_new);
01385   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01386   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01387 
01388   lwork = (int)work_[0];
01389   work_.resize(lwork);
01390   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01391   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01392 
01393   // Explicitly construct Q and F factors
01394   // NOTE:  The upper triangular part of HP is copied into F and HP becomes Q.
01395   SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
01396   for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
01397   lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01398   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01399 
01400   // Form orthonormalized C and adjust U accordingly so that C = A*U
01401   // C = [C V] * Q;
01402 
01403   // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01404   {
01405     Teuchos::RCP<MV> C1tmp;
01406     {
01407       index.resize(keff);
01408       for (int i=0; i < keff; i++) { index[i] = i; }
01409       Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_,  index );
01410       index.resize(keff_new);
01411       for (int i=0; i < keff_new; i++) { index[i] = i; }
01412       C1tmp  = MVT::CloneViewNonConst( *C1_,  index );
01413       SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
01414       MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
01415     }
01416     // Now compute C += V(:,1:p+1) * Q
01417     {
01418       index.resize( p );
01419       for (int i=0; i < p; ++i) { index[i] = i; }
01420       Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01421       SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
01422       MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
01423     }
01424   }
01425 
01426   // C_ = C1_; (via a swap)
01427   std::swap(C_, C1_);
01428 
01429   // Finally, compute U_ = U_*R^{-1}
01430   // First, compute LU factorization of R
01431   ipiv_.resize(Ftmp.numRows());
01432   lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
01433   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01434 
01435   // Now, form inv(R)
01436   lwork = Ftmp.numRows();
01437   work_.resize(lwork);
01438   lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01439   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01440 
01441   {
01442     index.resize(keff_new);
01443     for (int i=0; i < keff_new; i++) { index[i] = i; }
01444     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_,  index );
01445     MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
01446   }
01447 
01448   // Set the current number of recycled blocks and subspace 
01449   // dimension with the Block GCRO-DR iteration.
01450   if (keff != keff_new) {
01451     keff = keff_new;
01452     block_gcrodr_iter->setSize( keff, numBlocks_ );
01453     // Important to zero this out before next cyle
01454     SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01455     b1.putScalar(zero);
01456   }
01457 
01458 } //end buildRecycleSpaceAugKryl definition
01459 
01460 template<class ScalarType, class MV, class OP>
01461 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(int keff, int m, const SDM& GG, const Teuchos::RCP<const MV>& VV, SDM& PP){
01462   int i, j;
01463   int m2 = GG.numCols();
01464   bool xtraVec = false;
01465   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01466   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01467   std::vector<int> index;
01468 
01469   // Real and imaginary eigenvalue components
01470   std::vector<ScalarType> wr(m2), wi(m2);
01471 
01472   // Magnitude of harmonic Ritz values
01473   std::vector<MagnitudeType> w(m2);
01474 
01475   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01476   SDM vr(m2,m2,false);
01477 
01478   // Sorted order of harmonic Ritz values
01479   std::vector<int> iperm(m2);
01480 
01481   // Set flag indicating recycle space has been generated this solve
01482   builtRecycleSpace_ = true;
01483 
01484   // Form matrices for generalized eigenproblem
01485 
01486   // B = G' * G; Don't zero out matrix when constructing
01487   SDM B(m2,m2,false);
01488   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
01489 
01490   // A_tmp = | C'*U        0 |
01491   //         | V_{m+1}'*U  I |
01492   SDM A_tmp( keff+m+blockSize_, keff+m );
01493 
01494 
01495   // A_tmp(1:keff,1:keff) = C' * U;
01496   index.resize(keff);
01497   for (int i=0; i<keff; ++i) { index[i] = i; }
01498   Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_, index );
01499   Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01500   SDM A11( Teuchos::View, A_tmp, keff, keff );
01501   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
01502 
01503   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
01504   SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
01505   index.resize(m+blockSize_);
01506   for (i=0; i < m+blockSize_; i++) { index[i] = i; }
01507   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
01508   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
01509 
01510   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
01511   for( i=keff; i<keff+m; i++)
01512     A_tmp(i,i) = one;
01513 
01514   // A = G' * A_tmp;
01515   SDM A( m2, A_tmp.numCols() );
01516   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
01517 
01518   // Compute k smallest harmonic Ritz pairs
01519   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
01520   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
01521   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
01522   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
01523   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
01524   char balanc='P', jobvl='N', jobvr='V', sense='N';
01525   int ld = A.numRows();
01526   int lwork = 6*ld;
01527   int ldvl = ld, ldvr = ld;
01528   int info = 0,ilo = 0,ihi = 0;
01529   ScalarType abnrm = zero, bbnrm = zero;
01530   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
01531   std::vector<ScalarType> beta(ld);
01532   std::vector<ScalarType> work(lwork);
01533   std::vector<MagnitudeType> lscale(ld), rscale(ld);
01534   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
01535   std::vector<int> iwork(ld+6);
01536   int *bwork = 0; // If sense == 'N', bwork is never referenced
01537   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
01538   &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
01539   &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
01540   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
01541 
01542   // Construct magnitude of each harmonic Ritz value
01543   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
01544   for( i=0; i<ld; i++ ) // Construct magnitude of each harmonic Ritz value
01545     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
01546 
01547   this->sort(w,ld,iperm);
01548 
01549   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01550   if (wi[iperm[ld-recycledBlocks_]] != zero) {
01551     int countImag = 0;
01552     for ( i=ld-recycledBlocks_; i<ld; i++ )
01553       if (wi[iperm[i]] != zero) countImag++;
01554     // Check to see if this count is even or odd:
01555     if (countImag % 2) xtraVec = true;
01556   }
01557 
01558   // Select recycledBlocks_ smallest eigenvectors
01559 
01560   for( i=0; i<recycledBlocks_; i++ )
01561     for( j=0; j<ld; j++ )
01562       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
01563 
01564   if (xtraVec) { // we need to store one more vector
01565     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
01566       for( j=0; j<ld; j++ )                  // so get the "imag" component
01567         PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
01568     }
01569     else {                  // I picked the "imag" component
01570       for( j=0; j<ld; j++ ) // so get the "real" component
01571         PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
01572     }
01573   }
01574 
01575   // Return whether we needed to store an additional vector
01576   if (xtraVec)
01577     return recycledBlocks_+1;
01578   else
01579     return recycledBlocks_;
01580   
01581 } //end getHarmonicVecsAugKryl definition
01582 
01583 template<class ScalarType, class MV, class OP>
01584 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(int m, const SDM& HH, SDM& PP){
01585   bool xtraVec = false;
01586   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01587   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01588 
01589   // Real and imaginary eigenvalue components
01590   std::vector<MagnitudeType> wr(m), wi(m);
01591 
01592   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01593   SDM vr(m,m,false);
01594 
01595   // Magnitude of harmonic Ritz values
01596   std::vector<MagnitudeType> w(m);
01597 
01598   // Sorted order of harmonic Ritz values, also used for DGEEV
01599   std::vector<int> iperm(m);
01600 
01601   // Size of workspace and workspace for DGEEV
01602   int lwork = 4*m;
01603   std::vector<ScalarType> work(lwork);
01604 
01605   // Output info
01606   int info = 0;
01607 
01608   // Set flag indicating recycle space has been generated this solve
01609   builtRecycleSpace_ = true;
01610 
01611   // Solve linear system:  H_m^{-H}*E_m where E_m is the last blockSize_ columns of the identity matrix
01612   SDM HHt( HH, Teuchos::TRANS );
01613   Teuchos::RCP<SDM> harmRitzMatrix = rcp( new SDM( m, blockSize_));
01614 
01615   //Initialize harmRitzMatrix as E_m
01616   for(int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
01617 
01618   //compute harmRitzMatrix <- H_m^{-H}*E_m
01619   lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
01620 
01621   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01622   // Compute H_m + H_m^{-H}*E_m*H_lbl^{H}*H_lbl  
01623   // H_lbl is bottom-right block of H_, which is a blockSize_ x blockSize_ matrix
01624 
01625   Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
01626   Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl, Teuchos::TRANS );
01627   
01628   { //So that HTemp will fall out of scope
01629 
01630     // HH_lbl_t <- H_lbl_t*H_lbl
01631     Teuchos::RCP<SDM> Htemp = Teuchos::null;
01632     Htemp = Teuchos::rcp(new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
01633     Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
01634     H_lbl_t.assign(*Htemp);
01635     // harmRitzMatrix <- harmRitzMatrix*HH_lbl_t
01636     Htemp = Teuchos::rcp(new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
01637     Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
01638     harmRitzMatrix -> assign(*Htemp);
01639 
01640     // We need to add harmRitzMatrix to the last blockSize_ columns of HH and store the result harmRitzMatrix
01641     int harmColIndex, HHColIndex;
01642     Htemp = Teuchos::rcp(new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
01643     for(int i = 0; i<blockSize_; i++){
01644       harmColIndex = harmRitzMatrix -> numCols() - i -1;
01645       HHColIndex = m-i-1;
01646       for(int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
01647     } 
01648     harmRitzMatrix = Htemp;
01649   }
01650 
01651   // Revise to do query for optimal workspace first
01652   // Create simple storage for the left eigenvectors, which we don't care about.
01653 
01654   const int ldvl = m;
01655   ScalarType* vl = 0;
01656   lapack.GEEV('N', 'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
01657               vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01658   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01659 
01660   // Construct magnitude of each harmonic Ritz value
01661   for( int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01662 
01663   this->sort(w, m, iperm);
01664 
01665   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01666   if (wi[iperm[recycledBlocks_-1]] != zero) {
01667     int countImag = 0;
01668     for (int i=0; i<recycledBlocks_; ++i )
01669       if (wi[iperm[i]] != zero) countImag++;
01670     // Check to see if this count is even or odd:
01671     if (countImag % 2) xtraVec = true;
01672   }
01673 
01674   // Select recycledBlocks_ smallest eigenvectors
01675   for( int i=0; i<recycledBlocks_; ++i ) 
01676     for(int j=0; j<m; j++ )
01677       PP(j,i) = vr(j,iperm[i]);
01678 
01679   if (xtraVec) { // we need to store one more vector
01680     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
01681       for(int j=0; j<m; ++j )               // so get the "imag" component
01682         PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
01683     }
01684     else{                         // I picked the "imag" component
01685       for(int j=0; j<m; ++j )     // so get the "real" component
01686         PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
01687     }
01688   }
01689 
01690   // Return whether we needed to store an additional vector
01691   if (xtraVec) {
01692     printer_->stream(Debug) << "Recycled " << recycledBlocks_+1 << " vectors" << std::endl;
01693     return recycledBlocks_+1;
01694   }
01695   printer_->stream(Debug) << "Recycled " << recycledBlocks_ << " vectors" << std::endl;
01696   return recycledBlocks_;
01697 
01698 } //end getHarmonicVecsKryl
01699 
01700 
01701 // This method sorts list of n floating-point numbers and return permutation vector
01702 template<class ScalarType, class MV, class OP>
01703 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
01704   int l, r, j, i, flag;
01705   int    RR2;
01706   ScalarType dRR, dK;
01707 
01708   // Initialize the permutation vector.
01709   for(j=0;j<n;j++)
01710     iperm[j] = j;
01711 
01712   if (n <= 1) return;
01713 
01714   l    = n / 2 + 1;
01715   r    = n - 1;
01716   l    = l - 1;
01717   dRR  = dlist[l - 1];
01718   dK   = dlist[l - 1];
01719 
01720   RR2 = iperm[l - 1];
01721   while (r != 0) {
01722     j = l;
01723     flag = 1;
01724     while (flag == 1) {
01725       i = j;
01726       j = j + j;
01727       if (j > r + 1)
01728         flag = 0;
01729       else {
01730         if (j < r + 1)
01731           if (dlist[j] > dlist[j - 1]) j = j + 1;
01732         if (dlist[j - 1] > dK) {
01733           dlist[i - 1] = dlist[j - 1];
01734           iperm[i - 1] = iperm[j - 1];
01735         }
01736         else {
01737           flag = 0;
01738         }
01739       }
01740     }
01741     dlist[i - 1] = dRR;
01742     iperm[i - 1] = RR2;
01743     if (l == 1) {
01744       dRR  = dlist [r];
01745       RR2 = iperm[r];
01746       dK = dlist[r];
01747       dlist[r] = dlist[0];
01748       iperm[r] = iperm[0];
01749       r = r - 1;
01750     }
01751     else {
01752       l   = l - 1;
01753       dRR  = dlist[l - 1];
01754       RR2  = iperm[l - 1];
01755       dK   = dlist[l - 1];
01756     }
01757   } 
01758   dlist[0] = dRR;
01759   iperm[0] = RR2;
01760 } //end sort() definition
01761 
01762 template<class ScalarType, class MV, class OP>
01763 ReturnType BlockGCRODRSolMgr<ScalarType,MV,OP>::solve() {
01764   using Teuchos::RCP;
01765   using Teuchos::rcp;
01766   using Teuchos::rcp_const_cast;
01767 
01768   // MLP: NEED TO ADD CHECK IF PARAMETERS ARE SET LATER
01769      
01770   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01771   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01772   std::vector<int> index(numBlocks_+1);
01773 
01774   // MLP: EXCEPTION TESTING SHOULD GO HERE
01775 
01776   //THE FOLLOWING BLOCK OF CODE IS TO INFORM THE PROBLEM HOW MANY RIGHT HAND
01777   //SIDES WILL BE SOLVED.  IF THE CHOSEN BLOCK SIZE IS THE SAME AS THE NUMBER
01778   //OF RIGHT HAND SIDES IN THE PROBLEM THEN WE PASS THE PROBLEM 
01779   //INDICES FOR ALL RIGHT HAND SIDES.  IF BLOCK SIZE IS GREATER THAN
01780   //THE NUMBER OF RIGHT HAND SIDES AND THE USER HAS ADAPTIVE BLOCK SIZING  
01781   //TURNED OFF, THEN THE PROBLEM WILL GENERATE RANDOM RIGHT HAND SIDES SO WE HAVE AS MANY
01782   //RIGHT HAND SIDES AS BLOCK SIZE INDICATES.  THIS MAY NEED TO BE CHANGED
01783   //LATER TO ACCOMADATE SMARTER CHOOSING OF FICTITIOUS RIGHT HAND SIDES.
01784   //THIS CODE IS DIRECTLY LIFTED FROM BelosBlockGmresSolMgr.hpp.
01785      
01786   // Create indices for the linear systems to be solved.
01787   int startPtr = 0;
01788   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01789   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01790 
01791   //currIdx holds indices to all right-hand sides to be solved
01792   //or -1's to indicate that random right-hand sides should be generated
01793   std::vector<int> currIdx;
01794 
01795   if ( adaptiveBlockSize_ ) {
01796     blockSize_ = numCurrRHS;
01797     currIdx.resize( numCurrRHS  );
01798     for (int i=0; i<numCurrRHS; ++i)
01799       currIdx[i] = startPtr+i;
01800   } 
01801   else {
01802     currIdx.resize( blockSize_ );
01803     for (int i=0; i<numCurrRHS; ++i)
01804       currIdx[i] = startPtr+i;
01805     for (int i=numCurrRHS; i<blockSize_; ++i)
01806       currIdx[i] = -1;
01807   }
01808 
01809   // Inform the linear problem of the linear systems to solve/generate.
01810   problem_->setLSIndex( currIdx );
01811 
01812   // Check the number of blocks and change them is necessary.
01813   int dim = MVT::GetVecLength( *(problem_->getRHS()) );
01814      
01815   //ADD ERROR CHECKING TO MAKE SURE SIZE OF BLOCK KRYLOV SUBSPACE NOT LARGER THAN dim
01816  
01817   // reset loss of orthogonality flag
01818   loaDetected_ = false;
01819      
01820   // Assume convergence is achieved, then let any failed convergence set this to false.
01821   bool isConverged = true;
01822 
01823   // Initialize storage for all state variables
01824   initializeStateStorage(); 
01825 
01826   // Parameter list MLP: WE WILL NEED TO ASSIGN SOME PARAMETER VALUES LATER
01827   Teuchos::ParameterList plist;
01828   
01829   while (numRHS2Solve > 0){ //This loops through each block of RHS's being solved
01830     // **************************************************************************************************************************
01831     // Begin initial solver preparations.  Either update U,C for new operator or generate an initial space using a cycle of GMRES
01832     // **************************************************************************************************************************
01833     int prime_iterations;
01834     if(keff > 0){//If there is already a subspace to recycle, then use it
01835       // Update U, C for the new operator
01836 
01837 //      TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,BlockGCRODRSolMgrRecyclingFailure,
01838 //                                 "Belos::BlockGCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
01839                                  printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
01840 
01841       // Compute image of U_ under the new operator
01842       index.resize(keff);
01843       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01844       RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01845       RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01846       problem_->apply( *Utmp, *Ctmp );
01847 
01848       RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01849 
01850       // Orthogonalize this block
01851       // Get a matrix to hold the orthonormalization coefficients.
01852       SDM Ftmp( Teuchos::View, *F_, keff, keff );
01853       int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,false));
01854       // Throw an error if we could not orthogonalize this block
01855       TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,BlockGCRODRSolMgrOrthoFailure,"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
01856 
01857       // U_ = U_*F^{-1}
01858       // First, compute LU factorization of R
01859       int info = 0;
01860       ipiv_.resize(Ftmp.numRows());
01861       lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
01862       TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01863   
01864       // Now, form inv(F)
01865       int lwork = Ftmp.numRows();
01866       work_.resize(lwork);
01867       lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01868       TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01869 
01870       // U_ = U1_; (via a swap)
01871       MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
01872       std::swap(U_, U1_);
01873 
01874       // Must reinitialize after swap
01875       index.resize(keff);
01876       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01877       Ctmp  = MVT::CloneViewNonConst( *C_, index );
01878       Utmp  = MVT::CloneView( *U_, index );
01879   
01880       // Compute C_'*R_
01881       SDM Ctr(keff,blockSize_);
01882       problem_->computeCurrPrecResVec( &*R_ );
01883       MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
01884 
01885       // Update solution ( x += U_*C_'*R_ )
01886       RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
01887       MVT::MvInit( *update, 0.0 );
01888       MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
01889       problem_->updateSolution( update, true );
01890   
01891       // Update residual norm ( r -= C_*C_'*R_ )
01892       MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
01893 
01894       // We recycled space from previous call
01895       prime_iterations = 0;
01896 
01897       // Since we have a previous recycle space, we do not need block_gmres_iter
01898       // and we must allocate V ourselves. See comments in initialize() routine for
01899       // further explanation.
01900       if (V_ == Teuchos::null) {
01901         // Check if there is a multivector to clone from.
01902         Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
01903         V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
01904       }
01905       else{
01906         // Generate V_ by cloning itself ONLY if more space is needed.
01907         if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
01908           Teuchos::RCP<const MV> tmp = V_;
01909           V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
01910         } 
01911       } 
01912     } //end if(keff > 0)
01913     else{ //if there was no subspace to recycle, then build one
01914       printer_->stream(Debug) << " No recycled subspace available for RHS index " << std::endl << std::endl;
01915       
01916       Teuchos::ParameterList primeList;
01917       //we set this as numBlocks-1 because that is the Krylov dimension we want
01918       //so that the size of the Hessenberg created matches that of what we were expecting.
01919       primeList.set("Num Blocks",numBlocks_-1);
01920       primeList.set("Block Size",blockSize_);
01921       primeList.set("Recycled Blocks",0);
01922       primeList.set("Keep Hessenberg",true);
01923       primeList.set("Initialize Hessenberg",true);
01924       
01925       int dim = MVT::GetVecLength( *(problem_->getRHS()) );
01926       if (blockSize_*numBlocks_ > dim) {//if user has selected a total subspace dimension larger than system dimension
01927         int tmpNumBlocks = 0;
01928         if (blockSize_ == 1)
01929           tmpNumBlocks = dim / blockSize_;  // Allow for a good breakdown.
01930         else{
01931           tmpNumBlocks = ( dim - blockSize_) / blockSize_;  // Allow for restarting.
01932           printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve():  Warning! Requested Krylov subspace dimension is larger than operator dimension!"
01933                                      << std::endl << "The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
01934           primeList.set("Num Blocks",tmpNumBlocks);
01935         }
01936       }
01937       else{
01938         primeList.set("Num Blocks",numBlocks_-1);
01939       }
01940       //Create Block GMRES iteration object to perform one cycle of GMRES
01941       Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
01942       block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
01943 
01944       // MLP: ADD LOGIC TO DEAL WITH USER ASKING TO GENERATE A LARGER SPACE THAN dim AS IN HEIDI'S BlockGmresSolMgr CODE (DIDN'T WE ALREADY DO THIS SOMEWHERE?)
01945       block_gmres_iter->setSize( blockSize_, numBlocks_-1 );    
01946 
01947       //Create the first block in the current BLOCK Krylov basis (using residual)
01948       Teuchos::RCP<MV> V_0;
01949       if (currIdx[blockSize_-1] == -1) {//if augmented with random RHS's
01950         V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
01951         problem_->computeCurrPrecResVec( &*V_0 );
01952       }
01953       else {
01954         V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
01955       }
01956        
01957       // Get a matrix to hold the orthonormalization coefficients.
01958       Teuchos::RCP<SDM > z_0 = Teuchos::rcp( new SDM( blockSize_, blockSize_ ) );
01959       
01960       // Orthonormalize the new V_0
01961       int rank = ortho_->normalize( *V_0, z_0 );
01962       // MLP: ADD EXCEPTION IF INITIAL BLOCK IS RANK DEFFICIENT
01963 
01964       // Set the new state and initialize the iteration.
01965       GmresIterationState<ScalarType,MV> newstate;
01966       newstate.V = V_0;
01967       newstate.z = z_0;
01968       newstate.curDim = 0;
01969       block_gmres_iter->initializeGmres(newstate);
01970 
01971       bool primeConverged = false;
01972 
01973       try { 
01974         printer_->stream(Debug) << " Preparing to Iterate!!!!" << std::endl << std::endl;
01975         block_gmres_iter->iterate();
01976         
01977         // *********************
01978         // check for convergence
01979         // *********************
01980         if ( convTest_->getStatus() == Passed ) {
01981           printer_->stream(Debug) << "We converged during the prime the pump stage" << std::endl << std::endl;
01982           primeConverged = !(expConvTest_->getLOADetected());
01983           if ( expConvTest_->getLOADetected() ) {
01984             // we don't have convergence
01985             loaDetected_ = true;
01986             printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
01987           }
01988         }
01989         // *******************************************
01990         // check for maximum iterations of block GMRES
01991         // *******************************************
01992         else if( maxIterTest_->getStatus() == Passed ) {
01993           // we don't have convergence
01994           primeConverged = false;
01995         }
01996         // ****************************************************************
01997         // We need to recycle and continue, print a message indicating this
01998         // ****************************************************************
01999         else{ //debug statements
02000           printer_->stream(Debug) << " We did not converge on priming cycle of Block GMRES.  Therefore we recycle and restart. " << std::endl << std::endl;
02001         }
02002       } //end try
02003       catch (const GmresIterationOrthoFailure &e) {
02004         // If the block size is not one, it's not considered a lucky breakdown.
02005         if (blockSize_ != 1) {
02006           printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
02007                                    << block_gmres_iter->getNumIters() << std::endl
02008                                    << e.what() << std::endl;
02009          if (convTest_->getStatus() != Passed)
02010            primeConverged = false;
02011         }
02012         else {
02013           // If the block size is one, try to recover the most recent least-squares solution
02014           block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
02015           // Check to see if the most recent least-squares solution yielded convergence.
02016           sTest_->checkStatus( &*block_gmres_iter );
02017           if (convTest_->getStatus() != Passed)
02018             isConverged = false;
02019         }
02020       } // end catch (const GmresIterationOrthoFailure &e) 
02021       catch (const std::exception &e) {
02022         printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
02023                                  << block_gmres_iter->getNumIters() << std::endl
02024                                  << e.what() << std::endl;
02025         throw;
02026       }
02027       
02028       // Record number of iterations in generating initial recycle spacec
02029       prime_iterations = block_gmres_iter->getNumIters();//instantiated here because it is not needed outside of else{} scope;  we'll see if this is true or not
02030        
02031       // Update the linear problem.
02032       RCP<MV> update = block_gmres_iter->getCurrentUpdate();
02033       problem_->updateSolution( update, true );
02034       
02035       // Update the block residual
02036       
02037       problem_->computeCurrPrecResVec( &*R_ );
02038       
02039       // Get the state.
02040       newstate = block_gmres_iter->getState();
02041       int p = newstate.curDim;
02042       H_->assign(*(newstate.H));//IS THIS A GOOD IDEA?  I DID IT SO MEMBERS WOULD HAVE ACCESS TO H.
02043       // Get new Krylov vectors, store in V_
02044       
02045       // Since the block_gmres_iter returns the state
02046       // with const RCP's we need to cast newstate.V as
02047       // a non const RCP.  This is okay because block_gmres_iter
02048       // is about to fall out of scope, so we are simply 
02049       // rescuing *newstate.V from being destroyed so we can 
02050       // use it for future block recycled GMRES cycles
02051       V_ = rcp_const_cast<MV>(newstate.V);
02052       newstate.V.release();
02053       //COMPUTE NEW RECYCLE SPACE SOMEHOW
02054       buildRecycleSpaceKryl(keff, block_gmres_iter);
02055       printer_->stream(Debug) << "Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
02056       
02057       // Return to outer loop if the priming solve 
02058       // converged, set the next linear system.
02059       if (primeConverged) {
02060         /* MLP:  POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
02061         // Inform the linear problem that we are 
02062         // finished with this block linear system.
02063         problem_->setCurrLS();
02064         
02065         // Update indices for the linear systems to be solved.
02066         // KMS: Fix the numRHS2Solve; Copy from Heidi's BlockGmres code
02067         numRHS2Solve -= 1;
02068         printer_->stream(Debug) << numRHS2Solve << " Right hand sides left to solve" << std::endl;
02069         if ( numRHS2Solve > 0 ) {
02070           currIdx[0]++;
02071           // Set the next indices
02072           problem_->setLSIndex( currIdx );
02073         }
02074         else {
02075           currIdx.resize( numRHS2Solve );
02076         }
02077         *****************************************************************************/ 
02078         
02079         // Inform the linear problem that we are finished with this block linear system.
02080         problem_->setCurrLS();
02081         
02082         // Update indices for the linear systems to be solved.
02083         startPtr += numCurrRHS;
02084         numRHS2Solve -= numCurrRHS;
02085         if ( numRHS2Solve > 0 ) {
02086           numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
02087           if ( adaptiveBlockSize_ ) {
02088             blockSize_ = numCurrRHS;
02089             currIdx.resize( numCurrRHS  );
02090             for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
02091           }
02092           else {
02093             currIdx.resize( blockSize_ );
02094             for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
02095             for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
02096           }
02097           // Set the next indices.
02098           problem_->setLSIndex( currIdx );
02099         }
02100         else {
02101           currIdx.resize( numRHS2Solve );
02102         }
02103         continue;//Begin solving the next block of RHS's
02104       }  //end if (primeConverged)
02105     } //end else [if(keff > 0)]
02106 
02107     // *****************************************
02108     // Initial subspace update/construction done
02109     // Start cycles of recycled GMRES
02110     // *****************************************
02111     Teuchos::ParameterList blockgcrodrList;
02112     blockgcrodrList.set("Num Blocks",numBlocks_);
02113     blockgcrodrList.set("Block Size",blockSize_);
02114     blockgcrodrList.set("Recycled Blocks",keff);
02115 
02116     Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
02117     block_gcrodr_iter = Teuchos::rcp( new BlockGCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,blockgcrodrList) );  
02118     BlockGCRODRIterState<ScalarType,MV> newstate;
02119 
02120     index.resize( blockSize_ ); 
02121     for(int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
02122     Teuchos::RCP<MV> V0 =  MVT::CloneViewNonConst( *V_,  index );
02123 
02124     // MLP: MVT::SetBlock(*R_,index,*V0);
02125     MVT::Assign(*R_,*V0);
02126 
02127     index.resize(keff);//resize to get appropriate recycle space vectors
02128     for(int i=0; i < keff; i++){ index[i] = i;};
02129     B_ = rcp(new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
02130     H_ = rcp(new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
02131 
02132     newstate.V = V_;
02133     newstate.B= B_;
02134     newstate.U = MVT::CloneViewNonConst(*U_, index);
02135     newstate.C = MVT::CloneViewNonConst(*C_, index);
02136     newstate.H = H_;
02137     newstate.curDim = blockSize_;
02138     block_gcrodr_iter -> initialize(newstate);
02139 
02140     int numRestarts = 0;
02141 
02142     while(1){//Each execution of this loop is a cycle of block GCRODR
02143       try{  
02144         block_gcrodr_iter -> iterate();
02145 
02146         // ***********************
02147         // Check Convergence First
02148         // ***********************
02149         if( convTest_->getStatus() == Passed ) {
02150           // we have convergence
02151           break;//from while(1)
02152         } //end if converged  
02153 
02154         // ***********************************
02155         // Check if maximum iterations reached
02156         // ***********************************
02157         else if(maxIterTest_->getStatus() == Passed ){
02158           // no convergence; hit maxit
02159           isConverged = false;
02160           break; // from while(1)
02161         } //end elseif reached maxit
02162 
02163         // **********************************************
02164         // Check if subspace full; do we need to restart?
02165         // **********************************************
02166         else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
02167 
02168           // Update recycled space even if we have reached max number of restarts
02169 
02170           // Update linear problem
02171           Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
02172           problem_->updateSolution(update, true);
02173           buildRecycleSpaceAugKryl(block_gcrodr_iter);
02174 
02175           printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
02176           // NOTE: If we have hit the maximum number of restarts, then we will quit
02177           if(numRestarts >= maxRestarts_) {
02178             isConverged = false;
02179             break; //from while(1)
02180           } //end if max restarts
02181 
02182           numRestarts++;
02183 
02184           printer_ -> stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
02185 
02186           // Create the restart vector (first block in the current Krylov basis)
02187           problem_->computeCurrPrecResVec( &*R_ );
02188           index.resize( blockSize_ );
02189           for (int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
02190           Teuchos::RCP<MV> V0 =  MVT::CloneViewNonConst( *V_,  index );
02191           MVT::SetBlock(*R_,index,*V0);
02192 
02193           // Set the new state and initialize the solver.
02194           BlockGCRODRIterState<ScalarType,MV> restartState;
02195           index.resize( numBlocks_*blockSize_ );
02196           for (int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
02197           restartState.V  = MVT::CloneViewNonConst( *V_,  index );
02198           index.resize( keff );
02199           for (int ii=0; ii<keff; ++ii) index[ii] = ii;
02200           restartState.U  = MVT::CloneViewNonConst( *U_,  index );
02201           restartState.C  = MVT::CloneViewNonConst( *C_,  index );
02202           B_ = rcp(new SDM(Teuchos::View, *G_, keff,                  (numBlocks_-1)*blockSize_,     0, keff));
02203           H_ = rcp(new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
02204           restartState.B = B_;
02205           restartState.H = H_;
02206           restartState.curDim = blockSize_;
02207           block_gcrodr_iter->initialize(restartState);
02208 
02209         } //end else if need to restart
02210       
02211         // ****************************************************************
02212         // We returned from iterate(), but none of our status tests passed.
02213         // Something is wrong, and it is probably our fault.
02214         // ****************************************************************
02215         else {
02216           TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
02217         } //end else (no status test passed)
02218 
02219       }//end try
02220       catch(const BlockGCRODRIterOrthoFailure &e){ //there was an exception
02221         // Try to recover the most recent least-squares solution
02222         block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
02223         // Check to see if the most recent least-squares solution yielded convergence.
02224         sTest_->checkStatus( &*block_gcrodr_iter );
02225         if (convTest_->getStatus() != Passed) isConverged = false;
02226         break;
02227       }  // end catch orthogonalization failure
02228       catch(const std::exception &e){
02229         printer_->stream(Errors) << "Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
02230                                  << block_gcrodr_iter->getNumIters() << std::endl
02231                                  << e.what() << std::endl;
02232         throw;
02233       } //end catch standard exception
02234     } //end while(1)
02235 
02236     // Compute the current solution.
02237     // Update the linear problem.
02238     Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
02239     problem_->updateSolution( update, true );
02240 
02241     /* MLP:  POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
02242     // Inform the linear problem that we are finished with this block linear system.
02243     problem_->setCurrLS();
02244 
02245     // Update indices for the linear systems to be solved.
02246     numRHS2Solve -= 1;
02247     if ( numRHS2Solve > 0 ) {
02248       currIdx[0]++;
02249       // Set the next indices.
02250       problem_->setLSIndex( currIdx );
02251     }
02252     else {
02253       currIdx.resize( numRHS2Solve );
02254     }
02255     ******************************************************************************/
02256 
02257     // Inform the linear problem that we are finished with this block linear system.
02258     problem_->setCurrLS();
02259 
02260     // Update indices for the linear systems to be solved.
02261     startPtr += numCurrRHS;
02262     numRHS2Solve -= numCurrRHS;
02263     if ( numRHS2Solve > 0 ) {
02264       numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
02265       if ( adaptiveBlockSize_ ) {
02266         blockSize_ = numCurrRHS;
02267         currIdx.resize( numCurrRHS  );
02268         for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
02269       }
02270       else {
02271         currIdx.resize( blockSize_ );
02272         for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
02273         for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
02274       }
02275       // Set the next indices.
02276       problem_->setLSIndex( currIdx );
02277     }
02278     else {
02279       currIdx.resize( numRHS2Solve );
02280     }
02281 
02282     // If we didn't build a recycle space this solve but ran at least k iterations, force build of new recycle space
02283     if (!builtRecycleSpace_) {
02284       buildRecycleSpaceAugKryl(block_gcrodr_iter);
02285       printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
02286     }
02287   }//end while (numRHS2Solve > 0)
02288      
02289   // print final summary
02290   sTest_->print( printer_->stream(FinalSummary) );
02291 
02292   // print timing information
02293   #ifdef BELOS_TEUCHOS_TIME_MONITOR
02294   // Calling summarize() can be expensive, so don't call unless the user wants to print out timing details. 
02295   // summarize() will do all the work even if it's passed a "black hole" output stream.
02296   if (verbosity_ & TimingDetails) Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
02297   #endif
02298   // get iteration information for this solve
02299   numIters_ = maxIterTest_->getNumIters();
02300 
02301   if (!isConverged) return Unconverged; // return from BlockGCRODRSolMgr::solve()
02302     return Converged; // return from BlockGCRODRSolMgr::solve()
02303 } //end solve()
02304 
02305 } //End Belos Namespace
02306 
02307 #endif /* BELOS_BLOCK_GCRODR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines