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