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 #include "Teuchos_as.hpp"
00065 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00066 #include "Teuchos_TimeMonitor.hpp"
00067 #endif // BELOS_TEUCHOS_TIME_MONITOR
00068 
00069 // MLP: Add links to examples here
00070 
00071 namespace Belos{
00072 
00074 
00075 
00080   class BlockGCRODRSolMgrLinearProblemFailure : public BelosError {
00081   public:
00082     BlockGCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
00083   };
00084 
00089   class BlockGCRODRSolMgrOrthoFailure : public BelosError {
00090   public:
00091     BlockGCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00092   };
00093 
00098   class BlockGCRODRSolMgrLAPACKFailure : public BelosError {
00099   public:
00100     BlockGCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
00101   };
00102 
00107   class BlockGCRODRSolMgrRecyclingFailure : public BelosError {
00108   public:
00109     BlockGCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
00110   };
00111 
00113 
00125 
00126 template<class ScalarType, class MV, class OP>
00127 class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP> {
00128 private:
00129 
00130   typedef MultiVecTraits<ScalarType,MV> MVT;
00131   typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00132   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00133   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00134   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00135   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00136   typedef Teuchos::ScalarTraits<MagnitudeType> SMT;
00137   typedef OrthoManagerFactory<ScalarType, MV, OP> ortho_factory_type;
00138   typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
00139   typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
00140 
00141 public:
00143 
00144 
00150   BlockGCRODRSolMgr();
00151     
00177   BlockGCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00178          const Teuchos::RCP<Teuchos::ParameterList> &pl);
00179 
00181   virtual ~BlockGCRODRSolMgr() {};
00183   
00186 
00188   std::string description() const;      
00189 
00191 
00192 
00194 
00195     
00197   const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00198     return *problem_;
00199   }
00200 
00202   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00203 
00205   Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00206 
00208   int getNumIters() const { return numIters_; }
00209 
00211   bool isLOADetected() const { return loaDetected_; }
00212 
00214 
00216 
00217 
00219   void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) {
00220     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00221       "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
00222 
00223     // Check state of problem object before proceeding
00224     if (! problem->isProblemSet()) {
00225       const bool success = problem->setProblem();
00226       TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
00227         "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed.  This likely means that the "
00228         "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.");
00229     }
00230 
00231     problem_ = problem;
00232   }
00233 
00235   void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00236 
00238 
00240 
00241     
00248   void reset (const ResetType type) {
00249     if ((type & Belos::Problem) && ! problem_.is_null ())
00250       problem_->setProblem();
00251     else if (type & Belos::RecycleSubspace)
00252       keff = 0;
00253   }
00254 
00256 
00258 
00259 
00275   // \returns ::ReturnType specifying:
00278   ReturnType solve();
00279 
00281 
00282 private: 
00283 
00284   // Called by all constructors; Contains init instructions common to all constructors
00285   void init();
00286 
00287   // Initialize solver state storage
00288   void initializeStateStorage();
00289 
00290   // Recycling Methods
00291   // Appending Function name by:
00292   //  "Kryl" indicates it is specialized for building a recycle space after an 
00293   //         initial run of Block GMRES which generates an initial unaugmented block Krylov subspace
00294   //  "AugKryl" indicates  it is specialized for building a recycle space from the augmented Krylov subspace
00295 
00296   // Functions which control the building of a recycle space
00297   void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter);
00298   void buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
00299 
00300   // Recycling with Harmonic Ritz Vectors
00301   // Computes harmonic eigenpairs of projected matrix created during the priming solve.
00302   // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00303 
00304   // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension (numBlocks+1)*blockSize x numBlocks.
00305   // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00306   int getHarmonicVecsKryl (int m, const SDM& HH, SDM& PP);
00307 
00308   // 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.
00309   // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) (numBlocks+1)*blockSize vectors.
00310   // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00311   int getHarmonicVecsAugKryl (int keff, 
00312             int m,
00313             const SDM& HH,
00314             const Teuchos::RCP<const MV>& VV,
00315             SDM& PP);
00316 
00317   // Sort list of n floating-point numbers and return permutation vector
00318   void sort (std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00319 
00320   // Lapack interface
00321   Teuchos::LAPACK<int,ScalarType> lapack;
00322 
00324   Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00325 
00326   //Output Manager
00327   Teuchos::RCP<OutputManager<ScalarType> > printer_;
00328   Teuchos::RCP<std::ostream> outputStream_;
00329 
00330   //Status Test
00331   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00332   Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00333   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00334   Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00335   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00336 
00338   ortho_factory_type orthoFactory_;
00339 
00345   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00346 
00348   Teuchos::RCP<Teuchos::ParameterList> params_;
00349 
00360   mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
00361 
00362   //Default Solver Values
00363   static const bool adaptiveBlockSize_default_;
00364   static const std::string recycleMethod_default_;
00365 
00366   //Current Solver Values
00367   MagnitudeType convTol_, orthoKappa_;
00368   int blockSize_, maxRestarts_, maxIters_, numIters_;
00369   int verbosity_, outputStyle_, outputFreq_;
00370   bool adaptiveBlockSize_;
00371   std::string orthoType_, recycleMethod_;
00372   std::string impResScale_, expResScale_;
00373   std::string label_;
00374 
00376   // Solver State Storage
00378   //
00379   // The number of blocks and recycle blocks (m and k, respectively)
00380   int numBlocks_, recycledBlocks_;
00381   // Current size of recycled subspace
00382   int keff;
00383   //
00384   // Residual Vector
00385   Teuchos::RCP<MV> R_;
00386   //
00387   // Search Space
00388   Teuchos::RCP<MV> V_;
00389   //
00390   // Recycle subspace and its image under action of the operator
00391   Teuchos::RCP<MV> U_, C_;
00392   //
00393   // Updated recycle Space and its image under action of the operator
00394   Teuchos::RCP<MV> U1_, C1_;
00395   //
00396   // Storage used in constructing recycle space
00397   Teuchos::RCP<SDM > G_;
00398   Teuchos::RCP<SDM > H_;
00399   Teuchos::RCP<SDM > B_;
00400   Teuchos::RCP<SDM > PP_;
00401   Teuchos::RCP<SDM > HP_;
00402   std::vector<ScalarType> tau_;
00403   std::vector<ScalarType> work_;
00404   Teuchos::RCP<SDM > F_;
00405   std::vector<int> ipiv_;
00406 
00408   Teuchos::RCP<Teuchos::Time> timerSolve_;
00409 
00411   bool isSet_;
00412 
00414   bool loaDetected_;
00415 
00417   bool builtRecycleSpace_; 
00418 };
00419 
00420   //
00421   // Set default solver values
00422   //
00423   template<class ScalarType, class MV, class OP>
00424   const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ = true;
00425 
00426   template<class ScalarType, class MV, class OP>
00427   const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ = "harmvecs";
00428 
00429   // 
00430   // Method definitions 
00431   //
00432 
00433   template<class ScalarType, class MV, class OP>
00434   BlockGCRODRSolMgr<ScalarType,MV,OP>::BlockGCRODRSolMgr() {
00435     init();
00436   }
00437 
00438   //Basic Constructor
00439   template<class ScalarType, class MV, class OP>
00440   BlockGCRODRSolMgr<ScalarType,MV,OP>::
00441   BlockGCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00442         const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
00443     // Initialize local pointers to null, and initialize local
00444     // variables to default values.
00445     init ();
00446 
00447     TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
00448       "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
00449       "the linear problem argument 'problem' to be nonnull.");
00450 
00451     problem_ = problem;
00452 
00453     // Set the parameters using the list that was passed in.  If null, we defer initialization until a non-null list is set (by the
00454     // client calling setParameters(), or by calling solve() -- in either case, a null parameter list indicates that default parameters should be used).
00455     if (! pl.is_null())
00456       setParameters (pl);
00457   }
00458 
00459   template<class ScalarType, class MV, class OP>
00460   void BlockGCRODRSolMgr<ScalarType,MV,OP>::init() {
00461     adaptiveBlockSize_ = adaptiveBlockSize_default_;
00462     recycleMethod_ = recycleMethod_default_;
00463     isSet_ = false;
00464     loaDetected_ = false;
00465     builtRecycleSpace_ = false;
00466     keff = 0;//Effective Size of Recycle Space
00467     //The following are all RCP smart pointers to indicated matrices/vectors.
00468     //Some MATLAB notation used in comments.
00469     R_ = Teuchos::null;//The Block Residual
00470     V_ = Teuchos::null;//Block Arnoldi Vectors
00471     U_ = Teuchos::null;//Recycle Space
00472     C_ = Teuchos::null;//Image of U Under Action of Operator
00473     U1_ = Teuchos::null;//Newly Computed Recycle Space
00474     C1_ = Teuchos::null;//Image of Newly Computed Recycle Space
00475     PP_ = Teuchos::null;//Coordinates of New Recyc. Vectors in Augmented Space
00476     HP_ = Teuchos::null;//H_*PP_ or G_*PP_
00477     G_ = Teuchos::null;//G_ such that A*[U V(:,1:numBlocks_*blockSize_)] = [C V_]*G_
00478     F_ = Teuchos::null;//Upper Triangular portion of QR factorization for HP_
00479     H_ = Teuchos::null;//H_ such that A*V(:,1:numBlocks_*blockSize_) = V_*H_ + C_*C_'*V_
00480     B_ = Teuchos::null;//B_ = C_'*V_
00481     
00482     //THIS BLOCK OF CODE IS COMMENTED OUT AND PLACED ELSEWHERE IN THE CODE
00483 /*
00484     //WE TEMPORARILY INITIALIZE STATUS TESTS HERE FOR TESTING PURPOSES, BUT THEY SHOULD BE 
00485     //INITIALIZED SOMEWHERE ELSE, LIKE THE setParameters() FUNCTION
00486 
00487     //INSTANTIATE AND INITIALIZE TEST OBJECTS AS NEEDED
00488     if (maxIterTest_.is_null()){
00489       maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
00490     }
00491     //maxIterTest_->setMaxIters (maxIters_);
00492 
00493     //INSTANTIATE THE PRINTER
00494     if (printer_.is_null()) {
00495       printer_ = Teuchos::rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00496     }
00497 
00498     if (ortho_.is_null()) // || changedOrthoType) %%%In other codes, this is also triggered if orthogonalization type changed {
00499       // Create orthogonalization manager.  This requires that the OutputManager (printer_) already be initialized.
00500       Teuchos::RCP<const Teuchos::ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
00501       ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, orthoParams);
00502     }
00503 
00504     // Convenience typedefs
00505     typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00506     typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00507 
00508     if (impConvTest_.is_null()) {
00509       impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
00510       impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), Belos::TwoNorm);
00511       impConvTest_->setShowMaxResNormOnly( true );
00512     }
00513 
00514     if (expConvTest_.is_null()) {
00515       expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
00516       expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
00517       expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), Belos::TwoNorm);
00518       expConvTest_->setShowMaxResNormOnly( true );
00519     }
00520 
00521     if (convTest_.is_null()) {
00522       convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,  impConvTest_, expConvTest_));
00523     }
00524 
00525     sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
00526 
00527     StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
00528     outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
00529     Passed+Failed+Undefined);
00530 
00531  */
00532 
00533   }
00534 
00535   //  This method requires the solver manager to return a string that describes itself.
00536   template<class ScalarType, class MV, class OP>
00537   std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::description() const {
00538     std::ostringstream oss;
00539     oss << "Belos::BlockGCRODRSolMgr<" << SCT::name() << ", ...>";
00540     oss << "{";
00541     oss << "Ortho Type='"<<orthoType_ ;
00542     oss << ", Num Blocks=" <<numBlocks_;
00543     oss << ", Block Size=" <<blockSize_;
00544     oss << ", Num Recycle Blocks=" << recycledBlocks_;
00545     oss << ", Max Restarts=" << maxRestarts_;
00546     oss << "}";
00547     return oss.str();
00548   }
00549    
00550    template<class ScalarType, class MV, class OP>
00551    Teuchos::RCP<const Teuchos::ParameterList> 
00552    BlockGCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const {
00553      using Teuchos::ParameterList;
00554      using Teuchos::parameterList;
00555      using Teuchos::RCP;
00556      using Teuchos::rcpFromRef;
00557 
00558      if (defaultParams_.is_null()) {
00559        RCP<ParameterList> pl = parameterList ();
00560 
00561        const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
00562        const int maxRestarts = 1000;
00563        const int maxIters = 5000;
00564        const int blockSize = 2;
00565        const int numBlocks = 100;
00566        const int numRecycledBlocks = 25;
00567        const int verbosity = Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::StatusTestDetails;
00568        const Belos::OutputType outputStyle = Belos::General;
00569        const int outputFreq = 1;
00570        RCP<std::ostream> outputStream = rcpFromRef (std::cout);
00571        const std::string impResScale ("Norm of Preconditioned Initial Residual");
00572        const std::string expResScale ("Norm of Initial Residual");
00573        const std::string timerLabel ("Belos");
00574        const std::string orthoType ("DGKS");
00575        RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
00576        //const MagnitudeType orthoKappa = SCT::magnitude (SCT::eps());
00577        //
00578        // mfh 31 Dec 2011: Negative values of orthoKappa are ignored.
00579        // Setting this to a negative value by default ensures that
00580        // this parameter is only _not_ ignored if the user
00581        // specifically sets a valid value.
00582        const MagnitudeType orthoKappa = -SCT::one();
00583 
00584        // Set all the valid parameters and their default values.
00585        pl->set ("Convergence Tolerance", convTol,
00586     "The tolerance that the solver needs to achieve in order for "
00587     "the linear system(s) to be declared converged.  The meaning "
00588     "of this tolerance depends on the convergence test details.");
00589        pl->set("Maximum Restarts", maxRestarts,
00590          "The maximum number of restart cycles (not counting the first) "
00591          "allowed for each set of right-hand sides solved.");
00592        pl->set("Maximum Iterations", maxIters,
00593          "The maximum number of iterations allowed for each set of "
00594          "right-hand sides solved.");
00595        pl->set("Block Size", blockSize,
00596          "How many linear systems to solve at once.");
00597        pl->set("Num Blocks", numBlocks,
00598          "The maximum number of blocks allowed in the Krylov subspace "
00599          "for each set of right-hand sides solved.  This includes the "
00600          "initial block.  Total Krylov basis storage, not counting the "
00601          "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
00602        pl->set("Num Recycled Blocks", numRecycledBlocks,
00603          "The maximum number of vectors in the recycled subspace." );
00604        pl->set("Verbosity", verbosity,
00605          "What type(s) of solver information should be written "
00606          "to the output stream.");
00607        pl->set("Output Style", outputStyle,
00608          "What style is used for the solver information to write "
00609          "to the output stream.");
00610        pl->set("Output Frequency", outputFreq,
00611          "How often convergence information should be written "
00612          "to the output stream.");
00613        pl->set("Output Stream", outputStream,
00614          "A reference-counted pointer to the output stream where all "
00615          "solver output is sent.");
00616        pl->set("Implicit Residual Scaling", impResScale,
00617          "The type of scaling used in the implicit residual convergence test.");
00618        pl->set("Explicit Residual Scaling", expResScale,
00619          "The type of scaling used in the explicit residual convergence test.");
00620        pl->set("Timer Label", timerLabel,
00621          "The string to use as a prefix for the timer labels.");
00622        //  pl->set("Restart Timers", restartTimers_);
00623        pl->set("Orthogonalization", orthoType,
00624          "The orthogonalization method to use.  Valid options: " +
00625          orthoFactory_.validNamesString());
00626        pl->set ("Orthogonalization Parameters", *orthoParams,
00627     "Sublist of parameters specific to the orthogonalization method to use.");
00628        pl->set("Orthogonalization Constant", orthoKappa,
00629          "When using DGKS orthogonalization: the \"depTol\" constant, used "
00630          "to determine whether another step of classical Gram-Schmidt is "
00631          "necessary.  Otherwise ignored.  Nonpositive values are ignored.");
00632        defaultParams_ = pl;
00633      }
00634      return defaultParams_;
00635    }
00636    
00637    template<class ScalarType, class MV, class OP>
00638    void
00639    BlockGCRODRSolMgr<ScalarType,MV,OP>::
00640    setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params) {
00641      using Teuchos::isParameterType;
00642      using Teuchos::getParameter;
00643      using Teuchos::null;
00644      using Teuchos::ParameterList;
00645      using Teuchos::parameterList;
00646      using Teuchos::RCP;
00647      using Teuchos::rcp;
00648      using Teuchos::rcp_dynamic_cast;
00649      using Teuchos::rcpFromRef;
00650      using Teuchos::Exceptions::InvalidParameter;
00651      using Teuchos::Exceptions::InvalidParameterName;
00652      using Teuchos::Exceptions::InvalidParameterType;
00653 
00654      // The default parameter list contains all parameters that BlockGCRODRSolMgr understands, and none that it doesn't
00655      RCP<const ParameterList> defaultParams = getValidParameters();
00656 
00657      if (params.is_null()) {
00658        // Users really should supply a non-null input ParameterList,
00659        // so that we can fill it in.  However, if they really did give
00660        // us a null list, be nice and use default parameters.  In this
00661        // case, we don't have to do any validation.
00662        params_ = parameterList (*defaultParams);
00663      } 
00664      else {
00665        // Validate the user's parameter list and fill in any missing parameters with default values.
00666        //
00667        // Currently, validation is quite strict.  The following line
00668        // will throw an exception for misspelled or extra parameters.
00669        // However, it will _not_ throw an exception if there are
00670        // missing parameters: these will be filled in with their
00671        // default values.  There is additional discussion of other
00672        // validation strategies in the comments of this function for
00673        // Belos::GCRODRSolMgr.
00674        params->validateParametersAndSetDefaults (*defaultParams);
00675        // No side effects on the solver until after validation passes.
00676        params_ = params;
00677      }
00678 
00679      // 
00680      // Validate and set parameters relating to the Krylov subspace
00681      // dimensions and the number of blocks.  
00682      //
00683      // We try to read and validate all these parameters' values
00684      // before setting them in the solver.  This is because the
00685      // validity of each may depend on the others' values.  Our goal
00686      // is to avoid incorrect side effects, so we don't set the values
00687      // in the solver until we know they are right.
00688      //
00689      {
00690        const int maxRestarts = params_->get<int> ("Maximum Restarts");
00691        TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
00692          "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
00693          "must be nonnegative, but you specified a negative value of " 
00694          << maxRestarts << ".");
00695 
00696        const int maxIters = params_->get<int> ("Maximum Iterations");
00697        TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
00698          "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
00699          "must be positive, but you specified a nonpositive value of " 
00700          << maxIters << ".");
00701 
00702        const int numBlocks = params_->get<int> ("Num Blocks");
00703        TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
00704          "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
00705          "positive, but you specified a nonpositive value of " << numBlocks 
00706          << ".");
00707 
00708        const int blockSize = params_->get<int> ("Block Size");
00709        TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
00710          "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
00711          "positive, but you specified a nonpositive value of " << blockSize 
00712          << ".");
00713 
00714        const int recycledBlocks = params_->get<int> ("Num Recycled Blocks");
00715        TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
00716          "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
00717          "be positive, but you specified a nonpositive value of " 
00718          << recycledBlocks << ".");
00719        TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks, 
00720          std::invalid_argument, "Belos::BlockGCRODRSolMgr: The \"Num Recycled "
00721          "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
00722          "but you specified \"Num Recycled Blocks\" = " << recycledBlocks 
00723          << " and \"Num Blocks\" = " << numBlocks << ".");
00724   
00725        // Now that we've validated the various dimension-related
00726        // parameters, we can set them in the solvers.
00727        maxRestarts_ = maxRestarts;
00728        maxIters_ = maxIters;
00729        numBlocks_ = numBlocks;
00730        blockSize_ = blockSize;
00731        recycledBlocks_ = recycledBlocks;
00732      }
00733 
00734      // Check to see if the timer label changed.  If it did, update it in
00735      // the parameter list, and create a new timer with that label (if
00736      // Belos was compiled with timers enabled).
00737      {
00738        std::string tempLabel = params_->get<std::string> ("Timer Label");
00739        const bool labelChanged = (tempLabel != label_);
00740 
00741 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00742        std::string solveLabel = label_ + ": BlockGCRODRSolMgr total solve time";
00743        if (timerSolve_.is_null()) {
00744    // We haven't created a timer yet.
00745    timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00746        } else if (labelChanged) {
00747    // We created a timer before with a different label.  In that
00748    // case, clear the old timer and create a new timer with the
00749    // new label.  Clearing old timers prevents them from piling
00750    // up, since they are stored statically, possibly without
00751    // checking for duplicates.
00752    Teuchos::TimeMonitor::clearCounter (solveLabel);
00753    timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
00754        }
00755 #endif // BELOS_TEUCHOS_TIME_MONITOR
00756      }
00757 
00758      // Check for a change in verbosity level.  Just in case, we allow
00759      // the type of this parameter to be either int or MsgType (an
00760      // enum).
00761      if (params_->isParameter ("Verbosity")) {
00762        if (isParameterType<int> (*params_, "Verbosity")) {
00763    verbosity_ = params_->get<int> ("Verbosity");
00764        } 
00765        else {
00766    verbosity_ = (int) getParameter<MsgType> (*params_, "Verbosity");
00767        }
00768      }
00769 
00770      // Check for a change in output style.  Just in case, we allow
00771      // the type of this parameter to be either int or OutputType (an
00772      // enum).
00773      if (params_->isParameter ("Output Style")) {
00774        if (isParameterType<int> (*params_, "Output Style")) {
00775    outputStyle_ = params_->get<int> ("Output Style");
00776        } 
00777        else {
00778    outputStyle_ = (int) getParameter<OutputType> (*params_, "Output Style");
00779        }
00780 
00781        // Currently, the only way we can update the output style is to
00782        // create a new output test.  This is because we must first
00783        // instantiate a new StatusTestOutputFactory with the new
00784        // style, and then use the factory to make a new output test.
00785        // Thus, we set outputTest_ to null for now, so we remember
00786        // later to create a new one.
00787        outputTest_ = null;
00788      }
00789   
00790      // Get the output stream for the output manager.
00791      //
00792      // It has been pointed out (mfh 28 Feb 2011 in GCRODRSolMgr code)
00793      // that including an RCP<std::ostream> parameter makes it
00794      // impossible to serialize the parameter list.  If you serialize
00795      // the list and read it back in from the serialized
00796      // representation, you might not even get a valid
00797      // RCP<std::ostream>, let alone the same output stream that you
00798      // serialized.
00799      //
00800      // However, existing Belos users depend on setting the output
00801      // stream in the parameter list.  We retain this behavior for
00802      // backwards compatibility.
00803      //
00804      // In case the output stream can't be read back in, we default to
00805      // stdout (std::cout), just to ensure reasonable behavior.
00806      if (params_->isParameter ("Output Stream")) {
00807        try {
00808    outputStream_ = getParameter<RCP<std::ostream> > (*params_, "Output Stream");
00809        } 
00810        catch (InvalidParameter&) {
00811    outputStream_ = rcpFromRef (std::cout);
00812        }
00813        // We assume that a null output stream indicates that the user
00814        // doesn't want to print anything, so we replace it with a
00815        // "black hole" stream that prints nothing sent to it.  (We
00816        // can't use outputStream_ = Teuchos::null, since the output
00817        // manager assumes that operator<< works on its output stream.)
00818        if (outputStream_.is_null()) {
00819    outputStream_ = rcp (new Teuchos::oblackholestream);
00820        }
00821      }
00822 
00823      outputFreq_ = params_->get<int> ("Output Frequency");
00824      if (verbosity_ & Belos::StatusTestDetails) {
00825        // Update parameter in our output status test.
00826        if (! outputTest_.is_null()) {
00827    outputTest_->setOutputFrequency (outputFreq_);
00828        }
00829      }
00830 
00831      // Create output manager if we need to, using the verbosity level
00832      // and output stream that we fetched above.  We do this here because
00833      // instantiating an OrthoManager using OrthoManagerFactory requires
00834      // a valid OutputManager.
00835      if (printer_.is_null()) {
00836        printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00837      } else {
00838        printer_->setVerbosity (verbosity_);
00839        printer_->setOStream (outputStream_);
00840      }
00841 
00842      // Get the orthogonalization manager name ("Orthogonalization").
00843      //
00844      // Getting default values for the orthogonalization manager
00845      // parameters ("Orthogonalization Parameters") requires knowing the
00846      // orthogonalization manager name.  Save it for later, and also
00847      // record whether it's different than before.
00848      bool changedOrthoType = false;
00849      if (params_->isParameter ("Orthogonalization")) {
00850        const std::string& tempOrthoType = 
00851    params_->get<std::string> ("Orthogonalization");
00852        // Ensure that the specified orthogonalization type is valid.
00853        if (! orthoFactory_.isValidName (tempOrthoType)) {
00854    std::ostringstream os;
00855    os << "Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
00856       << tempOrthoType << "\".  The following are valid options "
00857       << "for the \"Orthogonalization\" name parameter: ";
00858    orthoFactory_.printValidNames (os);
00859    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
00860        }
00861        if (tempOrthoType != orthoType_) {
00862    changedOrthoType = true;
00863    orthoType_ = tempOrthoType;
00864        }
00865      }
00866 
00867      // Get any parameters for the orthogonalization
00868      // ("Orthogonalization Parameters").  If not supplied, the
00869      // orthogonalization manager factory will supply default values.
00870      //
00871      // NOTE (mfh 12 Jan 2011, summer 2011, 31 Dec 2011) For backwards
00872      // compatibility with other Belos GMRES-type solvers, if params
00873      // has an "Orthogonalization Constant" parameter and the DGKS
00874      // orthogonalization manager is to be used, the value of this
00875      // parameter will override DGKS's "depTol" parameter.
00876      //
00877      // Users must supply the orthogonalization manager parameters as
00878      // a sublist (supplying it as an RCP<ParameterList> would make
00879      // the resulting parameter list not serializable).
00880      RCP<ParameterList> orthoParams = sublist (params, "Orthogonalization Parameters", true);
00881      TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error, 
00882         "Failed to get orthogonalization parameters.  "
00883         "Please report this bug to the Belos developers.");
00884 
00885      // Check if the desired orthogonalization method changed, or if
00886      // the orthogonalization manager has not yet been instantiated.
00887      // If either is the case, instantiate a new MatOrthoManager
00888      // subclass instance corresponding to the desired
00889      // orthogonalization method.  We've already fetched the
00890      // orthogonalization method name (orthoType_) and its parameters
00891      // (orthoParams) above.
00892      //
00893      // As suggested (by mfh 12 Jan 2011 in Belos::GCRODRSolMgr) In
00894      // order to ensure parameter changes get propagated to the
00895      // orthomanager we simply reinstantiate the OrthoManager every
00896      // time, whether or not the orthogonalization method name or
00897      // parameters have changed.  This is not efficient for some
00898      // orthogonalization managers that keep a lot of internal state.
00899      // A more general way to fix this inefficiency would be to
00900      //
00901      // 1. make each orthogonalization implement 
00902      //    Teuchos::ParameterListAcceptor, and
00903      //
00904      // 2. make each orthogonalization implement a way to set the
00905      //    OutputManager instance and timer label.
00906      //
00907      // Create orthogonalization manager.  This requires that the
00908      // OutputManager (printer_) already be initialized.
00909      ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
00910              label_, orthoParams);
00911 
00912      // The DGKS orthogonalization accepts a "Orthogonalization
00913      // Constant" parameter (also called kappa in the code, but not in
00914      // the parameter list).  If its value is provided in the given
00915      // parameter list and its value is positive, use it.  Ignore
00916      // negative values.
00917      //
00918      // The "Orthogonalization Constant" parameter is retained only
00919      // for backwards compatibility.
00920      bool gotValidOrthoKappa = false;
00921      if (params_->isParameter ("Orthogonalization Constant")) {
00922        const MagnitudeType orthoKappa =
00923    params_->get<MagnitudeType> ("Orthogonalization Constant");
00924        if (orthoKappa > 0) {
00925    orthoKappa_ = orthoKappa;
00926    gotValidOrthoKappa = true;
00927    // Only DGKS currently accepts this parameter.
00928    if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
00929      typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
00930      // This cast should always succeed; it's a bug otherwise.
00931      // (If the cast fails, then orthoType_ doesn't correspond
00932      // to the OrthoManager subclass instance that we think we
00933      // have, so we initialized the wrong subclass somehow.)
00934      rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
00935    }
00936        }
00937      }
00938 
00939      // Convergence
00940      typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00941      typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00942 
00943      // Check for convergence tolerance
00944      convTol_ = params_->get<MagnitudeType> ("Convergence Tolerance");
00945      if (! impConvTest_.is_null()) {
00946        impConvTest_->setTolerance (convTol_);
00947      }
00948      if (! expConvTest_.is_null()) {
00949        expConvTest_->setTolerance (convTol_);
00950      }
00951 
00952      // Check for a change in scaling, if so we need to build new residual tests.
00953      if (params_->isParameter ("Implicit Residual Scaling")) {
00954        std::string tempImpResScale =
00955    getParameter<std::string> (*params_, "Implicit Residual Scaling");
00956 
00957        // Only update the scaling if it's different.
00958        if (impResScale_ != tempImpResScale) {
00959    ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
00960    impResScale_ = tempImpResScale;
00961 
00962    if (! impConvTest_.is_null()) {
00963      try {
00964        impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
00965      }
00966      catch (StatusTestError&) {
00967        // Delete the convergence test so it gets constructed again.
00968        impConvTest_ = null;
00969        convTest_ = null;
00970      }
00971    }
00972        }
00973      }
00974 
00975      if (params_->isParameter("Explicit Residual Scaling")) {
00976        std::string tempExpResScale =
00977    getParameter<std::string> (*params_, "Explicit Residual Scaling");
00978 
00979        // Only update the scaling if it's different.
00980        if (expResScale_ != tempExpResScale) {
00981    ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
00982    expResScale_ = tempExpResScale;
00983 
00984    if (! expConvTest_.is_null()) {
00985      try {
00986        expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
00987      }
00988      catch (StatusTestError&) {
00989        // Delete the convergence test so it gets constructed again.
00990        expConvTest_ = null;
00991        convTest_ = null;
00992      }
00993    }
00994        }
00995      }
00996 
00997      //
00998      // Create iteration stopping criteria ("status tests") if we need
00999      // to, by combining three different stopping criteria.
01000      //
01001      // First, construct maximum-number-of-iterations stopping criterion.
01002      if (maxIterTest_.is_null()) {
01003        maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
01004      } else {
01005        maxIterTest_->setMaxIters (maxIters_);
01006      }
01007 
01008      // Implicit residual test, using the native residual to determine if
01009      // convergence was achieved.
01010      if (impConvTest_.is_null()) {
01011        impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01012        impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
01013               Belos::TwoNorm);
01014      }
01015 
01016      // Explicit residual test once the native residual is below the tolerance
01017      if (expConvTest_.is_null()) {
01018        expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01019        expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
01020        expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
01021               Belos::TwoNorm);
01022      }
01023      // Convergence test first tests the implicit residual, then the
01024      // explicit residual if the implicit residual test passes.
01025      if (convTest_.is_null()) {
01026        convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
01027                  impConvTest_,
01028                  expConvTest_));
01029      }
01030      // Construct the complete iteration stopping criterion:
01031      //
01032      // "Stop iterating if the maximum number of iterations has been
01033      // reached, or if the convergence test passes."
01034      sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, 
01035             maxIterTest_, convTest_));
01036      // Create the status test output class.
01037      // This class manages and formats the output from the status test.
01038      StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
01039      outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
01040               Passed+Failed+Undefined);
01041 
01042      // Set the solver string for the output test
01043      std::string solverDesc = "Block GCRODR ";
01044      outputTest_->setSolverDesc (solverDesc);
01045 
01046      // Inform the solver manager that the current parameters were set.
01047      isSet_ = true;
01048    }
01049    
01050   // initializeStateStorage.
01051   template<class ScalarType, class MV, class OP>
01052   void 
01053   BlockGCRODRSolMgr<ScalarType,MV,OP>::initializeStateStorage()
01054   {
01055 
01056     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01057 
01058     // Check if there is any multivector to clone from.
01059     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
01060  
01061     //The Dimension of the Krylov Subspace 
01062     int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
01063 
01064     //Number of columns in [U_ V_(:,1:KrylSpaDim)]
01065     int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;// + 1 is for possible extra recycle vector
01066  
01067     //Number of columns in [C_ V_]
01068     int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
01069 
01070     //TEMPORARY SKELETON DEFINITION OF THIS FUNCTION TO GET THINGS WORKING
01071     //NOT EVERYTHING IS INITIALIZE CORRECTLY YET.
01072 
01073     //INITIALIZE RECYCLE SPACE VARIABLES HERE
01074 
01075     //WE DO NOT ALLOCATE V HERE IN THIS SOLVERMANAGER. If no recycle space exists, then block_gmres_iter
01076     //will allocated V for us.  If a recycle space already exists, then we will allocate V after updating the
01077     //recycle space for the current problem.
01078     // If the block Krylov subspace has not been initialized before, generate it using the RHS from lp_.
01079     /*if (V_ == Teuchos::null) {
01080       V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
01081       }
01082       else{
01083       // Generate V_ by cloning itself ONLY if more space is needed.
01084       if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
01085       Teuchos::RCP<const MV> tmp = V_;
01086       V_ = MVT::Clone( *tmp, numBlocks_+1 );
01087       }
01088       }*/
01089 
01090     //INTITIALIZE SPACE FOR THE NEWLY COMPUTED RECYCLE SPACE VARIABLES HERE
01091 
01092     if (U_ == Teuchos::null) {
01093       U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01094     }
01095     else {
01096       // Generate U_ by cloning itself ONLY if more space is needed.
01097       if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
01098   Teuchos::RCP<const MV> tmp = U_;
01099   U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01100       }
01101     }
01102 
01103     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01104     if (C_ == Teuchos::null) {
01105       C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01106     }
01107     else {
01108       // Generate C_ by cloning itself ONLY if more space is needed.
01109       if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
01110   Teuchos::RCP<const MV> tmp = C_;
01111   C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01112       }
01113     }
01114 
01115     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01116     if (U1_ == Teuchos::null) {
01117       U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01118     }
01119     else {
01120       // Generate U1_ by cloning itself ONLY if more space is needed.
01121       if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01122   Teuchos::RCP<const MV> tmp = U1_;
01123   U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01124       }
01125     }
01126 
01127     // If the subspace has not been initialized before, generate it using the RHS from lp_.
01128     if (C1_ == Teuchos::null) {
01129       C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01130     }
01131     else {
01132       // Generate C1_ by cloning itself ONLY if more space is needed.
01133       if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01134   Teuchos::RCP<const MV> tmp = C1_;
01135   C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01136       }
01137     }
01138 
01139     // Generate R_ only if it doesn't exist
01140     if (R_ == Teuchos::null){
01141       R_ = MVT::Clone( *rhsMV, blockSize_ );
01142     }
01143 
01144     //INITIALIZE SOME WORK VARIABLES
01145         
01146     // Generate G_ only if it doesn't exist, otherwise resize it.
01147     if (G_ == Teuchos::null){
01148       G_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
01149     }
01150     else{
01151       if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )       
01152   {
01153     G_->reshape( augSpaImgDim, augSpaDim );
01154   }
01155       G_->putScalar(zero);
01156     }
01157 
01158     // Generate H_ only if it doesn't exist by pointing it to a view of G_.
01159     if (H_ == Teuchos::null){
01160       H_ = Teuchos::rcp (new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
01161     }
01162 
01163     // Generate F_ only if it doesn't exist, otherwise resize it.
01164     if (F_ == Teuchos::null){
01165       F_ = Teuchos::rcp( new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
01166     }
01167     else {
01168       if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
01169   F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
01170       }
01171     }
01172     F_->putScalar(zero);
01173 
01174     // Generate PP_ only if it doesn't exist, otherwise resize it.
01175     if (PP_ == Teuchos::null){
01176       PP_ = Teuchos::rcp( new SDM( augSpaImgDim, recycledBlocks_+1 ) );
01177     }
01178     else{
01179       if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
01180   PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
01181       }
01182     }
01183 
01184     // Generate HP_ only if it doesn't exist, otherwise resize it.
01185     if (HP_ == Teuchos::null)
01186       HP_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
01187     else{
01188       if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
01189   HP_->reshape( augSpaImgDim, augSpaDim );
01190       }
01191     }
01192 
01193     // Size of tau_ will change during computation, so just be sure it starts with appropriate size
01194     tau_.resize(recycledBlocks_+1);
01195 
01196     // Size of work_ will change during computation, so just be sure it starts with appropriate size
01197     work_.resize(recycledBlocks_+1);
01198 
01199     // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
01200     ipiv_.resize(recycledBlocks_+1);
01201 
01202   }
01203 
01204 template<class ScalarType, class MV, class OP>
01205 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
01206 
01207   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01208   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01209 
01210 
01211   int p = block_gmres_iter->getState().curDim; //Dimension of the Krylov space generated
01212   std::vector<int> index(keff);//we use this to index certain columns of U, C, and V to 
01213   //get views into pieces of these matrices.
01214 
01215   //GET CORRECT PIECE OF MATRIX H APPROPRIATE TO SIZE OF KRYLOV SUBSPACE
01216   SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
01217   if(recycledBlocks_ >= p + blockSize_){//keep whole block Krylov subspace
01218     //IF THIS HAS HAPPENED, THIS MEANS WE CONVERGED DURING THIS CYCLE
01219     //THEREFORE, WE DO NOT CONSTRUCT C = A*U;
01220     index.resize(p);
01221     for (int ii=0; ii < p; ++ii) index[ii] = ii;
01222     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01223     MVT::SetBlock(*V_, index, *Utmp);
01224     keff = p;
01225   }
01226   else{ //use a subspace selection method to get recycle space
01227     int info = 0; 
01228     Teuchos::RCP<SDM > PPtmp = rcp (new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01229     if(recycleMethod_ == "harmvecs"){
01230       keff = getHarmonicVecsKryl(p, HH, *PPtmp);
01231       printer_->stream(Debug) << "keff = " << keff << std::endl;
01232     }
01233 // Hereafter, only keff columns of PP are needed
01234 PPtmp = rcp (new SDM ( Teuchos::View, *PP_, p, keff ) );
01235 // Now get views into C, U, V
01236 index.resize(keff);
01237 for (int ii=0; ii<keff; ++ii) index[ii] = ii;
01238 Teuchos::RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01239 Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01240 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01241 index.resize(p);
01242 for (int ii=0; ii < p; ++ii) index[ii] = ii;
01243 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01244 
01245 // Form U (the subspace to recycle)
01246 // U = newstate.V(:,1:p) * PP;
01247 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01248 
01249 // Step #1: Form HP = H*P
01250 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
01251 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
01252 // Step #1.5: Perform workspace size query for QR factorization of HP
01253 int lwork = -1;
01254 tau_.resize(keff);
01255 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01256 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01257 
01258 // Step #2: Compute QR factorization of HP
01259 lwork = (int)work_[0];
01260 work_.resize(lwork);
01261 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01262 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,  "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01263 
01264 // Step #3: Explicitly construct Q and R factors
01265 // The upper triangular part of HP is copied into R and HP becomes Q.
01266 SDM Rtmp( Teuchos::View, *F_, keff, keff );
01267 for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01268 lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01269 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01270     // Now we have [Q,R] = qr(H*P)
01271 
01272             // Now compute C = V(:,1:p+blockSize_) * Q
01273                 index.resize( p+blockSize_ );
01274             for (int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
01275             Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+blockSize_ vectors now; needed p above)
01276             MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01277 
01278     // Finally, compute U = U*R^{-1}.
01279           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
01280                 // backsolve capabilities don't exist in the Belos::MultiVec class
01281 
01282 // Step #1: First, compute LU factorization of R
01283 ipiv_.resize(Rtmp.numRows());
01284 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01285 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01286 // Step #2: Form inv(R)
01287 lwork = Rtmp.numRows();
01288 work_.resize(lwork);
01289 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01290 //TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01291 // Step #3: Let U = U * R^{-1}
01292 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01293 
01294 } //end else from if(recycledBlocks_ >= p + 1)
01295 return;
01296 } // end buildRecycleSpaceKryl defnition
01297 
01298 template<class ScalarType, class MV, class OP>
01299 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
01300   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01301   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01302 
01303   std::vector<MagnitudeType> d(keff);
01304   std::vector<int> index(numBlocks_+1);
01305 
01306   // Get the state
01307   BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
01308   int p = oldState.curDim;
01309 
01310   // insufficient new information to update recycle space
01311   if (p<1) return;
01312 
01313   if(recycledBlocks_ >= keff + p){ //we add new Krylov vectors to existing recycle space
01314      // If this has happened, we have converged during this cycle. Therefore, do not construct C = A*C;
01315      index.resize(p);
01316      for (int ii=0; ii < p; ++ii) { index[ii] = keff+ii; } //get a view after current reycle vectors
01317      Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01318      for (int ii=0; ii < p; ++ii) { index[ii] =ii; }
01319      MVT::SetBlock(*V_, index, *Utmp);
01320      keff += p;
01321   }
01322 
01323   // Take the norm of the recycled vectors
01324   {
01325     index.resize(keff);
01326     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01327     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01328     d.resize(keff);
01329     MVT::MvNorm( *Utmp, d );
01330     for (int i=0; i<keff; ++i)
01331       d[i] = one / d[i];
01332     MVT::MvScale( *Utmp, d );
01333   }
01334 
01335   // Get view into current "full" upper Hessnburg matrix
01336   // 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
01337   Teuchos::RCP<SDM> Gtmp = Teuchos::rcp( new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
01338 
01339   // Insert D into the leading keff x keff  block of G 
01340   for (int i=0; i<keff; ++i)
01341     (*Gtmp)(i,i) = d[i];
01342 
01343   // Compute the harmoic Ritz pairs for the generalized eigenproblem
01344   // getHarmonicVecsKryl assumes PP has recycledBlocks_+1 columns available
01345   // See previous block of comments for why we subtract p-blockSize_
01346   int keff_new;
01347   {
01348     SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
01349     keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
01350   }
01351 
01352   // Code to form new U, C
01353   // U = [U V(:,1:p)] * P; (in two steps)
01354   
01355   // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01356   Teuchos::RCP<MV> U1tmp;
01357   {
01358     index.resize( keff );
01359     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01360     Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
01361     index.resize( keff_new );
01362     for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01363     U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
01364     SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
01365     MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
01366   }
01367 
01368   // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01369   {
01370     index.resize(p-blockSize_);
01371     for (int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
01372     Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01373     SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
01374     MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
01375   }
01376 
01377   // Form GP = G*P
01378   SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
01379   {
01380     SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
01381     HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
01382   }
01383 
01384   // Workspace size query for QR factorization HP= QF (the worksize will be placed in work_[0])
01385   int info = 0, lwork = -1;
01386   tau_.resize(keff_new);
01387   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01388   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01389 
01390   lwork = (int)work_[0];
01391   work_.resize(lwork);
01392   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01393   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01394 
01395   // Explicitly construct Q and F factors
01396   // NOTE:  The upper triangular part of HP is copied into F and HP becomes Q.
01397   SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
01398   for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
01399   lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01400   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01401 
01402   // Form orthonormalized C and adjust U accordingly so that C = A*U
01403   // C = [C V] * Q;
01404 
01405   // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01406   {
01407     Teuchos::RCP<MV> C1tmp;
01408     {
01409       index.resize(keff);
01410       for (int i=0; i < keff; i++) { index[i] = i; }
01411       Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_,  index );
01412       index.resize(keff_new);
01413       for (int i=0; i < keff_new; i++) { index[i] = i; }
01414       C1tmp  = MVT::CloneViewNonConst( *C1_,  index );
01415       SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
01416       MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
01417     }
01418     // Now compute C += V(:,1:p+1) * Q
01419     {
01420       index.resize( p );
01421       for (int i=0; i < p; ++i) { index[i] = i; }
01422       Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01423       SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
01424       MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
01425     }
01426   }
01427 
01428   // C_ = C1_; (via a swap)
01429   std::swap(C_, C1_);
01430 
01431   // Finally, compute U_ = U_*R^{-1}
01432   // First, compute LU factorization of R
01433   ipiv_.resize(Ftmp.numRows());
01434   lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
01435   TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01436 
01437   // Now, form inv(R)
01438   lwork = Ftmp.numRows();
01439   work_.resize(lwork);
01440   lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01441   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01442 
01443   {
01444     index.resize(keff_new);
01445     for (int i=0; i < keff_new; i++) { index[i] = i; }
01446     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_,  index );
01447     MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
01448   }
01449 
01450   // Set the current number of recycled blocks and subspace 
01451   // dimension with the Block GCRO-DR iteration.
01452   if (keff != keff_new) {
01453     keff = keff_new;
01454     block_gcrodr_iter->setSize( keff, numBlocks_ );
01455     // Important to zero this out before next cyle
01456     SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01457     b1.putScalar(zero);
01458   }
01459 
01460 } //end buildRecycleSpaceAugKryl definition
01461 
01462 template<class ScalarType, class MV, class OP>
01463 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(int keff, int m, const SDM& GG, const Teuchos::RCP<const MV>& VV, SDM& PP){
01464   int i, j;
01465   int m2 = GG.numCols();
01466   bool xtraVec = false;
01467   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01468   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01469   std::vector<int> index;
01470 
01471   // Real and imaginary eigenvalue components
01472   std::vector<ScalarType> wr(m2), wi(m2);
01473 
01474   // Magnitude of harmonic Ritz values
01475   std::vector<MagnitudeType> w(m2);
01476 
01477   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01478   SDM vr(m2,m2,false);
01479 
01480   // Sorted order of harmonic Ritz values
01481   std::vector<int> iperm(m2);
01482 
01483   // Set flag indicating recycle space has been generated this solve
01484   builtRecycleSpace_ = true;
01485 
01486   // Form matrices for generalized eigenproblem
01487 
01488   // B = G' * G; Don't zero out matrix when constructing
01489   SDM B(m2,m2,false);
01490   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
01491 
01492   // A_tmp = | C'*U        0 |
01493   //         | V_{m+1}'*U  I |
01494   SDM A_tmp( keff+m+blockSize_, keff+m );
01495 
01496 
01497   // A_tmp(1:keff,1:keff) = C' * U;
01498   index.resize(keff);
01499   for (int i=0; i<keff; ++i) { index[i] = i; }
01500   Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_, index );
01501   Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01502   SDM A11( Teuchos::View, A_tmp, keff, keff );
01503   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
01504 
01505   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
01506   SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
01507   index.resize(m+blockSize_);
01508   for (i=0; i < m+blockSize_; i++) { index[i] = i; }
01509   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
01510   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
01511 
01512   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
01513   for( i=keff; i<keff+m; i++)
01514     A_tmp(i,i) = one;
01515 
01516   // A = G' * A_tmp;
01517   SDM A( m2, A_tmp.numCols() );
01518   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
01519 
01520   // Compute k smallest harmonic Ritz pairs
01521   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
01522   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
01523   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
01524   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
01525   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
01526   char balanc='P', jobvl='N', jobvr='V', sense='N';
01527   int ld = A.numRows();
01528   int lwork = 6*ld;
01529   int ldvl = ld, ldvr = ld;
01530   int info = 0,ilo = 0,ihi = 0;
01531   ScalarType abnrm = zero, bbnrm = zero;
01532   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
01533   std::vector<ScalarType> beta(ld);
01534   std::vector<ScalarType> work(lwork);
01535   std::vector<MagnitudeType> lscale(ld), rscale(ld);
01536   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
01537   std::vector<int> iwork(ld+6);
01538   int *bwork = 0; // If sense == 'N', bwork is never referenced
01539   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
01540   &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
01541   &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
01542   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
01543 
01544   // Construct magnitude of each harmonic Ritz value
01545   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
01546   for( i=0; i<ld; i++ ) // Construct magnitude of each harmonic Ritz value
01547     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
01548 
01549   this->sort(w,ld,iperm);
01550 
01551   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01552   if (wi[iperm[ld-recycledBlocks_]] != zero) {
01553     int countImag = 0;
01554     for ( i=ld-recycledBlocks_; i<ld; i++ )
01555       if (wi[iperm[i]] != zero) countImag++;
01556     // Check to see if this count is even or odd:
01557     if (countImag % 2) xtraVec = true;
01558   }
01559 
01560   // Select recycledBlocks_ smallest eigenvectors
01561 
01562   for( i=0; i<recycledBlocks_; i++ )
01563     for( j=0; j<ld; j++ )
01564       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
01565 
01566   if (xtraVec) { // we need to store one more vector
01567     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
01568       for( j=0; j<ld; j++ )                  // so get the "imag" component
01569         PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
01570     }
01571     else {                  // I picked the "imag" component
01572       for( j=0; j<ld; j++ ) // so get the "real" component
01573         PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
01574     }
01575   }
01576 
01577   // Return whether we needed to store an additional vector
01578   if (xtraVec)
01579     return recycledBlocks_+1;
01580   else
01581     return recycledBlocks_;
01582   
01583 } //end getHarmonicVecsAugKryl definition
01584 
01585 template<class ScalarType, class MV, class OP>
01586 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(int m, const SDM& HH, SDM& PP){
01587   bool xtraVec = false;
01588   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01589   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01590 
01591   // Real and imaginary eigenvalue components
01592   std::vector<MagnitudeType> wr(m), wi(m);
01593 
01594   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01595   SDM vr(m,m,false);
01596 
01597   // Magnitude of harmonic Ritz values
01598   std::vector<MagnitudeType> w(m);
01599 
01600   // Sorted order of harmonic Ritz values, also used for DGEEV
01601   std::vector<int> iperm(m);
01602 
01603   // Size of workspace and workspace for DGEEV
01604   int lwork = 4*m;
01605   std::vector<ScalarType> work(lwork);
01606 
01607   // Output info
01608   int info = 0;
01609 
01610   // Set flag indicating recycle space has been generated this solve
01611   builtRecycleSpace_ = true;
01612 
01613   // Solve linear system:  H_m^{-H}*E_m where E_m is the last blockSize_ columns of the identity matrix
01614   SDM HHt( HH, Teuchos::TRANS );
01615   Teuchos::RCP<SDM> harmRitzMatrix = rcp( new SDM( m, blockSize_));
01616 
01617   //Initialize harmRitzMatrix as E_m
01618   for(int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
01619 
01620   //compute harmRitzMatrix <- H_m^{-H}*E_m
01621   lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
01622 
01623   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01624   // Compute H_m + H_m^{-H}*E_m*H_lbl^{H}*H_lbl  
01625   // H_lbl is bottom-right block of H_, which is a blockSize_ x blockSize_ matrix
01626 
01627   Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
01628   Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl, Teuchos::TRANS );
01629   
01630   { //So that HTemp will fall out of scope
01631 
01632     // HH_lbl_t <- H_lbl_t*H_lbl
01633     Teuchos::RCP<SDM> Htemp = Teuchos::null;
01634     Htemp = Teuchos::rcp(new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
01635     Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
01636     H_lbl_t.assign(*Htemp);
01637     // harmRitzMatrix <- harmRitzMatrix*HH_lbl_t
01638     Htemp = Teuchos::rcp(new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
01639     Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
01640     harmRitzMatrix -> assign(*Htemp);
01641 
01642     // We need to add harmRitzMatrix to the last blockSize_ columns of HH and store the result harmRitzMatrix
01643     int harmColIndex, HHColIndex;
01644     Htemp = Teuchos::rcp(new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
01645     for(int i = 0; i<blockSize_; i++){
01646       harmColIndex = harmRitzMatrix -> numCols() - i -1;
01647       HHColIndex = m-i-1;
01648       for(int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
01649     } 
01650     harmRitzMatrix = Htemp;
01651   }
01652 
01653   // Revise to do query for optimal workspace first
01654   // Create simple storage for the left eigenvectors, which we don't care about.
01655 
01656   const int ldvl = m;
01657   ScalarType* vl = 0;
01658   lapack.GEEV('N', 'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
01659               vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01660   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01661 
01662   // Construct magnitude of each harmonic Ritz value
01663   for( int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01664 
01665   this->sort(w, m, iperm);
01666 
01667   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01668   if (wi[iperm[recycledBlocks_-1]] != zero) {
01669     int countImag = 0;
01670     for (int i=0; i<recycledBlocks_; ++i )
01671       if (wi[iperm[i]] != zero) countImag++;
01672     // Check to see if this count is even or odd:
01673     if (countImag % 2) xtraVec = true;
01674   }
01675 
01676   // Select recycledBlocks_ smallest eigenvectors
01677   for( int i=0; i<recycledBlocks_; ++i ) 
01678     for(int j=0; j<m; j++ )
01679       PP(j,i) = vr(j,iperm[i]);
01680 
01681   if (xtraVec) { // we need to store one more vector
01682     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
01683       for(int j=0; j<m; ++j )               // so get the "imag" component
01684         PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
01685     }
01686     else{                         // I picked the "imag" component
01687       for(int j=0; j<m; ++j )     // so get the "real" component
01688         PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
01689     }
01690   }
01691 
01692   // Return whether we needed to store an additional vector
01693   if (xtraVec) {
01694     printer_->stream(Debug) << "Recycled " << recycledBlocks_+1 << " vectors" << std::endl;
01695     return recycledBlocks_+1;
01696   }
01697   printer_->stream(Debug) << "Recycled " << recycledBlocks_ << " vectors" << std::endl;
01698   return recycledBlocks_;
01699 
01700 } //end getHarmonicVecsKryl
01701 
01702 
01703 // This method sorts list of n floating-point numbers and return permutation vector
01704 template<class ScalarType, class MV, class OP>
01705 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
01706   int l, r, j, i, flag;
01707   int    RR2;
01708   ScalarType dRR, dK;
01709 
01710   // Initialize the permutation vector.
01711   for(j=0;j<n;j++)
01712     iperm[j] = j;
01713 
01714   if (n <= 1) return;
01715 
01716   l    = n / 2 + 1;
01717   r    = n - 1;
01718   l    = l - 1;
01719   dRR  = dlist[l - 1];
01720   dK   = dlist[l - 1];
01721 
01722   RR2 = iperm[l - 1];
01723   while (r != 0) {
01724     j = l;
01725     flag = 1;
01726     while (flag == 1) {
01727       i = j;
01728       j = j + j;
01729       if (j > r + 1)
01730         flag = 0;
01731       else {
01732         if (j < r + 1)
01733           if (dlist[j] > dlist[j - 1]) j = j + 1;
01734         if (dlist[j - 1] > dK) {
01735           dlist[i - 1] = dlist[j - 1];
01736           iperm[i - 1] = iperm[j - 1];
01737         }
01738         else {
01739           flag = 0;
01740         }
01741       }
01742     }
01743     dlist[i - 1] = dRR;
01744     iperm[i - 1] = RR2;
01745     if (l == 1) {
01746       dRR  = dlist [r];
01747       RR2 = iperm[r];
01748       dK = dlist[r];
01749       dlist[r] = dlist[0];
01750       iperm[r] = iperm[0];
01751       r = r - 1;
01752     }
01753     else {
01754       l   = l - 1;
01755       dRR  = dlist[l - 1];
01756       RR2  = iperm[l - 1];
01757       dK   = dlist[l - 1];
01758     }
01759   } 
01760   dlist[0] = dRR;
01761   iperm[0] = RR2;
01762 } //end sort() definition
01763 
01764 template<class ScalarType, class MV, class OP>
01765 ReturnType BlockGCRODRSolMgr<ScalarType,MV,OP>::solve() {
01766   using Teuchos::RCP;
01767   using Teuchos::rcp;
01768   using Teuchos::rcp_const_cast;
01769 
01770   // MLP: NEED TO ADD CHECK IF PARAMETERS ARE SET LATER
01771      
01772   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01773   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01774   std::vector<int> index(numBlocks_+1);
01775 
01776   // MLP: EXCEPTION TESTING SHOULD GO HERE
01777 
01778   //THE FOLLOWING BLOCK OF CODE IS TO INFORM THE PROBLEM HOW MANY RIGHT HAND
01779   //SIDES WILL BE SOLVED.  IF THE CHOSEN BLOCK SIZE IS THE SAME AS THE NUMBER
01780   //OF RIGHT HAND SIDES IN THE PROBLEM THEN WE PASS THE PROBLEM 
01781   //INDICES FOR ALL RIGHT HAND SIDES.  IF BLOCK SIZE IS GREATER THAN
01782   //THE NUMBER OF RIGHT HAND SIDES AND THE USER HAS ADAPTIVE BLOCK SIZING  
01783   //TURNED OFF, THEN THE PROBLEM WILL GENERATE RANDOM RIGHT HAND SIDES SO WE HAVE AS MANY
01784   //RIGHT HAND SIDES AS BLOCK SIZE INDICATES.  THIS MAY NEED TO BE CHANGED
01785   //LATER TO ACCOMADATE SMARTER CHOOSING OF FICTITIOUS RIGHT HAND SIDES.
01786   //THIS CODE IS DIRECTLY LIFTED FROM BelosBlockGmresSolMgr.hpp.
01787      
01788   // Create indices for the linear systems to be solved.
01789   int startPtr = 0;
01790   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01791   int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
01792 
01793   //currIdx holds indices to all right-hand sides to be solved
01794   //or -1's to indicate that random right-hand sides should be generated
01795   std::vector<int> currIdx;
01796 
01797   if ( adaptiveBlockSize_ ) {
01798     blockSize_ = numCurrRHS;
01799     currIdx.resize( numCurrRHS  );
01800     for (int i=0; i<numCurrRHS; ++i)
01801       currIdx[i] = startPtr+i;
01802   } 
01803   else {
01804     currIdx.resize( blockSize_ );
01805     for (int i=0; i<numCurrRHS; ++i)
01806       currIdx[i] = startPtr+i;
01807     for (int i=numCurrRHS; i<blockSize_; ++i)
01808       currIdx[i] = -1;
01809   }
01810 
01811   // Inform the linear problem of the linear systems to solve/generate.
01812   problem_->setLSIndex( currIdx );
01813 
01814   //ADD ERROR CHECKING TO MAKE SURE SIZE OF BLOCK KRYLOV SUBSPACE NOT LARGER THAN dim
01815   //ptrdiff_t dim = MVText::GetGlobalLength( *(problem_->getRHS()) );
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       ptrdiff_t dim = MVText::GetGlobalLength( *(problem_->getRHS()) );
01926       if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {//if user has selected a total subspace dimension larger than system dimension
01927         ptrdiff_t 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",Teuchos::as<int>(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