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