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