BelosGCRODRSolMgr.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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_GCRODR_SOLMGR_HPP
00030 #define BELOS_GCRODR_SOLMGR_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041 
00042 #include "BelosGCRODRIter.hpp"
00043 #include "BelosBlockFGmresIter.hpp"
00044 #include "BelosDGKSOrthoManager.hpp"
00045 #include "BelosICGSOrthoManager.hpp"
00046 #include "BelosIMGSOrthoManager.hpp"
00047 #include "BelosStatusTestMaxIters.hpp"
00048 #include "BelosStatusTestGenResNorm.hpp"
00049 #include "BelosStatusTestCombo.hpp"
00050 #include "BelosStatusTestOutputFactory.hpp"
00051 #include "BelosOutputManager.hpp"
00052 #include "Teuchos_BLAS.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055 
00072 namespace Belos {
00073   
00075 
00076   
00083   class GCRODRSolMgrLinearProblemFailure : public BelosError {
00084     public:
00085       GCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
00086   };
00087   
00094   class GCRODRSolMgrOrthoFailure : public BelosError {
00095     public:
00096       GCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00097   };
00098   
00105   class GCRODRSolMgrLAPACKFailure : public BelosError {
00106     public:
00107       GCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
00108   };
00109 
00116   class GCRODRSolMgrRecyclingFailure : public BelosError {
00117     public:
00118       GCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
00119   };
00120 
00122 
00123   template<class ScalarType, class MV, class OP>
00124   class GCRODRSolMgr : public SolverManager<ScalarType,MV,OP> {
00125     
00126   private:
00127     typedef MultiVecTraits<ScalarType,MV> MVT;
00128     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00129     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00130     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00131     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00132     
00133   public:
00134     
00136 
00137    
00143      GCRODRSolMgr();
00144  
00158     GCRODRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, const Teuchos::RCP<Teuchos::ParameterList> &pl );
00159     
00161     virtual ~GCRODRSolMgr() {};
00163     
00165 
00166     
00169     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00170       return *problem_;
00171     }
00172 
00175     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00176 
00179     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { 
00180       return params_; 
00181     }
00182  
00188     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00189       return tuple(timerSolve_);
00190     }
00191   
00193     int getNumIters() const {
00194       return numIters_;
00195     }
00196  
00199     bool isLOADetected() const { return false; }
00200  
00202     
00204 
00205    
00207     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { 
00208       problem_ = problem; 
00209     }
00210    
00212     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00213     
00215    
00217 
00218 
00222     void reset( const ResetType type ) { 
00223       if ((type & Belos::Problem) && !Teuchos::is_null(problem_))  
00224         problem_->setProblem();
00225     }
00227  
00229 
00230     
00247     ReturnType solve();
00248     
00250     
00253     
00255     std::string description() const;
00256     
00258     
00259   private:
00260 
00261     // Called by all constructors; Contains init instructions common to all constructors
00262     void init();
00263 
00264     // Initialize solver state storage
00265     void initializeStateStorage();
00266 
00267     //  Computes harmonic eigenpairs of projected matrix created during the priming solve.
00268     //  HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
00269     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00270     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00271     int getHarmonicVecs1(int m, 
00272        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00273        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00274 
00275     //  Computes harmonic eigenpairs of projected matrix created during one cycle.
00276     //  HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
00277     //  VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
00278     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00279     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00280     int getHarmonicVecs2(int keff, int m, 
00281        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00282        const Teuchos::RCP<const MV>& VV,
00283        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00284 
00285     // Sort list of n floating-point numbers and return permutation vector
00286     void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00287 
00288     // Method to convert string to enumerated type for residual.
00289     Belos::ScaleType convertStringToScaleType( std::string& scaleType ) {
00290       if (scaleType == "Norm of Initial Residual") {
00291   return Belos::NormOfInitRes;
00292       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00293   return Belos::NormOfPrecInitRes;
00294       } else if (scaleType == "Norm of RHS") {
00295   return Belos::NormOfRHS;
00296       } else if (scaleType == "None") {
00297   return Belos::None;
00298       } else 
00299   TEST_FOR_EXCEPTION( true ,std::logic_error,"Belos::GCRODRSolMgr(): Invalid residual scaling type.");
00300     }
00301 
00302     // Lapack interface
00303     Teuchos::LAPACK<int,ScalarType> lapack;
00304 
00305     // Linear problem.
00306     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00307     
00308     // Output manager.
00309     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00310     Teuchos::RCP<std::ostream> outputStream_;
00311 
00312     // Status test.
00313     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00314     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00315     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00316     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00317     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00318 
00319     // Orthogonalization manager.
00320     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00321     
00322     // Current parameter list.
00323     Teuchos::RCP<ParameterList> params_;
00324 
00325     // Default solver values.
00326     static const MagnitudeType convtol_default_;
00327     static const MagnitudeType orthoKappa_default_;
00328     static const int maxRestarts_default_;
00329     static const int maxIters_default_;
00330     static const int numBlocks_default_;
00331     static const int blockSize_default_;
00332     static const int recycledBlocks_default_;
00333     static const int verbosity_default_;
00334     static const int outputStyle_default_;
00335     static const int outputFreq_default_;
00336     static const std::string impResScale_default_; 
00337     static const std::string expResScale_default_; 
00338     static const std::string label_default_;
00339     static const std::string orthoType_default_;
00340     static const Teuchos::RCP<std::ostream> outputStream_default_;
00341 
00342     // Current solver values.
00343     MagnitudeType convtol_, orthoKappa_;
00344     int maxRestarts_, maxIters_, numIters_;
00345     int verbosity_, outputStyle_, outputFreq_;
00346     std::string orthoType_; 
00347     std::string impResScale_, expResScale_;
00348 
00350     // Solver State Storage
00352     // 
00353     // The number of blocks and recycle blocks (m and k, respectively)
00354     int numBlocks_, recycledBlocks_;
00355     // Current size of recycled subspace 
00356     int keff;
00357     //
00358     // Residual vector
00359     Teuchos::RCP<MV> r_;
00360     //
00361     // Search space
00362     Teuchos::RCP<MV> V_;
00363     //
00364     // Recycled subspace and its image
00365     Teuchos::RCP<MV> U_, C_;
00366     //
00367     // Updated recycle space and its image
00368     Teuchos::RCP<MV> U1_, C1_;
00369     //
00370     // Storage used in constructing harmonic Ritz values/vectors
00371     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
00372     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00373     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00374     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
00375     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
00376     std::vector<ScalarType> tau_;
00377     std::vector<ScalarType> work_;
00378     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00379     std::vector<int> ipiv_;
00381 
00382     // Timers.
00383     std::string label_;
00384     Teuchos::RCP<Teuchos::Time> timerSolve_;
00385 
00386     // Internal state variables.
00387     bool isSet_;
00388   };
00389 
00390 
00391 // Default solver values.
00392 template<class ScalarType, class MV, class OP>
00393 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00394 
00395 template<class ScalarType, class MV, class OP>
00396 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00397 
00398 template<class ScalarType, class MV, class OP>
00399 const int GCRODRSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 100;
00400 
00401 template<class ScalarType, class MV, class OP>
00402 const int GCRODRSolMgr<ScalarType,MV,OP>::maxIters_default_ = 5000;
00403 
00404 template<class ScalarType, class MV, class OP>
00405 const int GCRODRSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 50;
00406 
00407 template<class ScalarType, class MV, class OP>
00408 const int GCRODRSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00409 
00410 template<class ScalarType, class MV, class OP>
00411 const int GCRODRSolMgr<ScalarType,MV,OP>::recycledBlocks_default_ = 5;
00412 
00413 template<class ScalarType, class MV, class OP>
00414 const int GCRODRSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00415 
00416 template<class ScalarType, class MV, class OP>
00417 const int GCRODRSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00418 
00419 template<class ScalarType, class MV, class OP>
00420 const int GCRODRSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00421 
00422 template<class ScalarType, class MV, class OP>
00423 const std::string GCRODRSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00424 
00425 template<class ScalarType, class MV, class OP>
00426 const std::string GCRODRSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00427 
00428 template<class ScalarType, class MV, class OP>
00429 const std::string GCRODRSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00430 
00431 template<class ScalarType, class MV, class OP>
00432 const std::string GCRODRSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00433 
00434 template<class ScalarType, class MV, class OP>
00435 const Teuchos::RCP<std::ostream> GCRODRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00436 
00437 
00438 // Empty Constructor
00439 template<class ScalarType, class MV, class OP>
00440 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr() {
00441   init();
00442 }
00443 
00444 
00445 // Basic Constructor
00446 template<class ScalarType, class MV, class OP>
00447 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr( 
00448                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00449                const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
00450 
00451   // initialize local pointers to null and local vars to default values
00452   init();
00453 
00454   TEST_FOR_EXCEPTION(problem == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00455   problem_ = problem;
00456 
00457   // set the parameters using the list that was passed in.
00458   if (!is_null(pl))
00459     setParameters( pl );  
00460 }
00461 
00462 // Common instructions executed in all constructors
00463 template<class ScalarType, class MV, class OP>
00464 void GCRODRSolMgr<ScalarType,MV,OP>::init() {
00465   outputStream_ = outputStream_default_;
00466   convtol_ = convtol_default_;
00467   orthoKappa_ = orthoKappa_default_;
00468   maxRestarts_ = maxRestarts_default_;
00469   maxIters_ = maxIters_default_;
00470   numBlocks_ = numBlocks_default_;
00471   recycledBlocks_ = recycledBlocks_default_;
00472   verbosity_ = verbosity_default_;
00473   outputStyle_ = outputStyle_default_;
00474   outputFreq_ = outputFreq_default_;
00475   orthoType_ = orthoType_default_;
00476   impResScale_ = impResScale_default_;
00477   expResScale_ = expResScale_default_;
00478   label_ = label_default_;
00479   isSet_ = false;
00480   numBlocks_ = 0;
00481   recycledBlocks_ = 0;
00482   keff = 0;
00483   r_ = Teuchos::null;
00484   V_ = Teuchos::null;
00485   U_ = Teuchos::null; 
00486   C_ = Teuchos::null;
00487   U1_ = Teuchos::null;
00488   C1_ = Teuchos::null;
00489   PP_ = Teuchos::null;
00490   HP_ = Teuchos::null;
00491   H2_ = Teuchos::null;
00492   R_ = Teuchos::null;
00493   H_ = Teuchos::null;
00494   B_ = Teuchos::null;
00495 }
00496 
00497 template<class ScalarType, class MV, class OP>
00498 void GCRODRSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00499 {
00500   // Create the internal parameter list if ones doesn't already exist.
00501   if (params_ == Teuchos::null) {
00502     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00503   }
00504   else {
00505     params->validateParameters(*getValidParameters());
00506   }
00507 
00508   // Check for maximum number of restarts
00509   if (params->isParameter("Maximum Restarts")) {
00510     maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
00511 
00512     // Update parameter in our list.
00513     params_->set("Maximum Restarts", maxRestarts_);
00514   }
00515 
00516   // Check for maximum number of iterations
00517   if (params->isParameter("Maximum Iterations")) {
00518     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00519 
00520     // Update parameter in our list and in status test.
00521     params_->set("Maximum Iterations", maxIters_);
00522     if (maxIterTest_!=Teuchos::null)
00523       maxIterTest_->setMaxIters( maxIters_ );
00524   }
00525 
00526   // Check for the maximum number of blocks.
00527   if (params->isParameter("Num Blocks")) {
00528     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00529     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00530            "Belos::GCRODRSolMgr: \"Num Blocks\" must be strictly positive.");
00531 
00532     // Update parameter in our list.
00533     params_->set("Num Blocks", numBlocks_);
00534   }
00535 
00536   // Check for the maximum number of blocks.
00537   if (params->isParameter("Num Recycled Blocks")) {
00538     recycledBlocks_ = params->get("Num Recycled Blocks",recycledBlocks_default_);
00539     TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
00540            "Belos::GCRODRSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
00541 
00542     TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
00543            "Belos::GCRODRSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
00544 
00545     // Update parameter in our list.
00546     params_->set("Num Recycled Blocks", recycledBlocks_);
00547   }
00548 
00549   // Check to see if the timer label changed.
00550   if (params->isParameter("Timer Label")) {
00551     std::string tempLabel = params->get("Timer Label", label_default_);
00552 
00553     // Update parameter in our list and solver timer
00554     if (tempLabel != label_) {
00555       label_ = tempLabel;
00556       params_->set("Timer Label", label_);
00557       std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00558       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00559     }
00560   }
00561 
00562   // Check if the orthogonalization changed.
00563   if (params->isParameter("Orthogonalization")) {
00564     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00565     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00566       std::invalid_argument,
00567       "Belos::GCRODRSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00568     if (tempOrthoType != orthoType_) {
00569       orthoType_ = tempOrthoType;
00570       // Create orthogonalization manager
00571       if (orthoType_=="DGKS") {
00572   if (orthoKappa_ <= 0) {
00573     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00574   }
00575   else {
00576     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00577     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00578   }
00579       }
00580       else if (orthoType_=="ICGS") {
00581   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00582       } 
00583       else if (orthoType_=="IMGS") {
00584   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00585       } 
00586     }  
00587   }
00588 
00589   // Check which orthogonalization constant to use.
00590   if (params->isParameter("Orthogonalization Constant")) {
00591     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00592 
00593     // Update parameter in our list.
00594     params_->set("Orthogonalization Constant",orthoKappa_);
00595     if (orthoType_=="DGKS") {
00596       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00597   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00598       }
00599     } 
00600   }
00601 
00602   // Check for a change in verbosity level
00603   if (params->isParameter("Verbosity")) {
00604     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00605       verbosity_ = params->get("Verbosity", verbosity_default_);
00606     } else {
00607       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00608     }
00609 
00610     // Update parameter in our list.
00611     params_->set("Verbosity", verbosity_);
00612     if (printer_ != Teuchos::null)
00613       printer_->setVerbosity(verbosity_);
00614   }
00615 
00616   // Check for a change in output style
00617   if (params->isParameter("Output Style")) {
00618     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00619       outputStyle_ = params->get("Output Style", outputStyle_default_);
00620     } else {
00621       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00622     }
00623 
00624     // Update parameter in our list.
00625     params_->set("Output Style", outputStyle_);
00626     outputTest_ = Teuchos::null;
00627   }
00628 
00629   // output stream
00630   if (params->isParameter("Output Stream")) {
00631     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00632 
00633     // Update parameter in our list.
00634     params_->set("Output Stream", outputStream_);
00635     if (printer_ != Teuchos::null)
00636       printer_->setOStream( outputStream_ );
00637   }
00638 
00639   // frequency level
00640   if (verbosity_ & Belos::StatusTestDetails) {
00641     if (params->isParameter("Output Frequency")) {
00642       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00643     }
00644 
00645     // Update parameter in out list and output status test.
00646     params_->set("Output Frequency", outputFreq_);
00647     if (outputTest_ != Teuchos::null)
00648       outputTest_->setOutputFrequency( outputFreq_ );
00649   }
00650 
00651   // Create output manager if we need to.
00652   if (printer_ == Teuchos::null) {
00653     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00654   }  
00655   
00656   // Convergence
00657   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00658   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00659 
00660   // Check for convergence tolerance
00661   if (params->isParameter("Convergence Tolerance")) {
00662     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00663 
00664     // Update parameter in our list and residual tests.
00665     params_->set("Convergence Tolerance", convtol_);
00666     if (impConvTest_ != Teuchos::null)
00667       impConvTest_->setTolerance( convtol_ );
00668     if (expConvTest_ != Teuchos::null)
00669       expConvTest_->setTolerance( convtol_ );
00670   }
00671  
00672   // Check for a change in scaling, if so we need to build new residual tests.
00673   if (params->isParameter("Implicit Residual Scaling")) {
00674     std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
00675 
00676     // Only update the scaling if it's different.
00677     if (impResScale_ != tempImpResScale) {
00678       Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
00679       impResScale_ = tempImpResScale;
00680 
00681       // Update parameter in our list and residual tests
00682       params_->set("Implicit Residual Scaling", impResScale_);
00683       if (impConvTest_ != Teuchos::null) {
00684         try { 
00685           impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
00686         }
00687         catch (std::exception& e) { 
00688           // Delete the convergence test so it gets constructed again.
00689     impConvTest_ = Teuchos::null;
00690           convTest_ = Teuchos::null;
00691         }
00692       }
00693     }      
00694   }
00695   
00696   if (params->isParameter("Explicit Residual Scaling")) {
00697     std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
00698 
00699     // Only update the scaling if it's different.
00700     if (expResScale_ != tempExpResScale) {
00701       Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
00702       expResScale_ = tempExpResScale;
00703 
00704       // Update parameter in our list and residual tests
00705       params_->set("Explicit Residual Scaling", expResScale_);
00706       if (expConvTest_ != Teuchos::null) {
00707         try { 
00708           expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
00709         }
00710         catch (std::exception& e) {
00711           // Delete the convergence test so it gets constructed again.
00712     expConvTest_ = Teuchos::null;
00713           convTest_ = Teuchos::null;
00714         }
00715       }
00716     }      
00717   }
00718 
00719   // Create status tests if we need to.
00720 
00721   // Basic test checks maximum iterations and native residual.
00722   if (maxIterTest_ == Teuchos::null)
00723     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00724 
00725   // Implicit residual test, using the native residual to determine if convergence was achieved.
00726   if (impConvTest_ == Teuchos::null) {
00727     impConvTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_ ) );
00728     impConvTest_->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
00729   }
00730 
00731   // Explicit residual test once the native residual is below the tolerance
00732   if (expConvTest_ == Teuchos::null) {
00733     expConvTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_ ) );
00734     expConvTest_->defineResForm( StatusTestResNorm_t::Explicit, Belos::TwoNorm );
00735     expConvTest_->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
00736   }
00737 
00738   if (convTest_ == Teuchos::null) {
00739     convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
00740   }
00741 
00742   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00743   
00744   // Create the status test output class.
00745   // This class manages and formats the output from the status test.
00746   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00747   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00748 
00749   // Set the solver string for the output test
00750   std::string solverDesc = " GCRODR ";
00751   outputTest_->setSolverDesc( solverDesc );
00752 
00753   // Create orthogonalization manager if we need to.
00754   if (ortho_ == Teuchos::null) {
00755     if (orthoType_=="DGKS") {
00756       if (orthoKappa_ <= 0) {
00757   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00758       }
00759       else {
00760   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00761   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00762       }
00763     }
00764     else if (orthoType_=="ICGS") {
00765       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00766     } 
00767     else if (orthoType_=="IMGS") {
00768       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00769     } 
00770     else {
00771       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00772        "Belos::GCRODRSolMgr(): Invalid orthogonalization type.");
00773     }  
00774   }
00775 
00776   // Create the timer if we need to.
00777   if (timerSolve_ == Teuchos::null) {
00778     std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00779     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00780   }
00781 
00782   // Inform the solver manager that the current parameters were set.
00783   isSet_ = true;
00784 }
00785 
00786     
00787 template<class ScalarType, class MV, class OP>
00788 Teuchos::RCP<const Teuchos::ParameterList> GCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const {
00789 
00790   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00791   if (is_null(validPL)) {
00792     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00793     // Set all the valid parameters and their default values.
00794     pl->set("Convergence Tolerance", convtol_default_,
00795       "The relative residual tolerance that needs to be achieved by the\n"
00796       "iterative solver in order for the linear system to be declared converged.");
00797     pl->set("Maximum Restarts", maxRestarts_default_,
00798       "The maximum number of cycles allowed for each\n"
00799       "set of RHS solved.");
00800     pl->set("Maximum Iterations", maxIters_default_,
00801       "The maximum number of iterations allowed for each\n"
00802       "set of RHS solved.");
00803     pl->set("Block Size", blockSize_default_,
00804       "Block Size Parameter -- currently must be 1 for GCRODR");
00805     pl->set("Num Blocks", numBlocks_default_,
00806       "The maximum number of vectors allowed in the Krylov subspace\n"
00807       "for each set of RHS solved.");
00808     pl->set("Num Recycled Blocks", recycledBlocks_default_,
00809       "The maximum number of vectors in the recycled subspace." );
00810     pl->set("Verbosity", verbosity_default_,
00811       "What type(s) of solver information should be outputted\n"
00812       "to the output stream.");
00813     pl->set("Output Style", outputStyle_default_,
00814       "What style is used for the solver information outputted\n"
00815       "to the output stream.");
00816     pl->set("Output Frequency", outputFreq_default_,
00817       "How often convergence information should be outputted\n"
00818       "to the output stream.");  
00819     pl->set("Output Stream", outputStream_default_,
00820       "A reference-counted pointer to the output stream where all\n"
00821       "solver output is sent.");
00822     pl->set("Implicit Residual Scaling", impResScale_default_,
00823       "The type of scaling used in the implicit residual convergence test.");
00824     pl->set("Explicit Residual Scaling", expResScale_default_,
00825       "The type of scaling used in the explicit residual convergence test.");
00826     pl->set("Timer Label", label_default_,
00827       "The string to use as a prefix for the timer labels.");
00828     //  pl->set("Restart Timers", restartTimers_);
00829     pl->set("Orthogonalization", orthoType_default_,
00830       "The type of orthogonalization to use: DGKS, ICGS, IMGS");
00831     pl->set("Orthogonalization Constant",orthoKappa_default_,
00832       "The constant used by DGKS orthogonalization to determine\n"
00833       "whether another step of classical Gram-Schmidt is necessary.");
00834     validPL = pl;
00835   }
00836   return validPL;
00837 
00838 }
00839 
00840 // initializeStateStorage
00841 template<class ScalarType, class MV, class OP>
00842 void GCRODRSolMgr<ScalarType,MV,OP>::initializeStateStorage() {
00843 
00844     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00845 
00846     // Check if there is any multivector to clone from.
00847     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
00848     if (rhsMV == Teuchos::null) {
00849       // Nothing to do
00850       return;
00851     }
00852     else {
00853 
00854       // Initialize the state storage
00855       TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
00856                          "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
00857 
00858       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00859       if (U_ == Teuchos::null) {
00860         U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
00861       }
00862       else {
00863         // Generate U_ by cloning itself ONLY if more space is needed.
00864         if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
00865           Teuchos::RCP<const MV> tmp = U_;
00866           U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
00867         }
00868       }
00869 
00870       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00871       if (C_ == Teuchos::null) {
00872         C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
00873       }
00874       else {
00875         // Generate C_ by cloning itself ONLY if more space is needed.
00876         if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
00877           Teuchos::RCP<const MV> tmp = C_;
00878           C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
00879         }
00880       }
00881 
00882       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00883       if (V_ == Teuchos::null) {
00884         V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
00885       }
00886       else {
00887         // Generate V_ by cloning itself ONLY if more space is needed.
00888         if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
00889           Teuchos::RCP<const MV> tmp = V_;
00890           V_ = MVT::Clone( *tmp, numBlocks_+1 );
00891         }
00892       }
00893 
00894       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00895       if (U1_ == Teuchos::null) {
00896         U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
00897       }
00898       else {
00899         // Generate U1_ by cloning itself ONLY if more space is needed.
00900         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
00901           Teuchos::RCP<const MV> tmp = U1_;
00902           U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
00903         }
00904       }
00905 
00906       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00907       if (C1_ == Teuchos::null) {
00908         C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
00909       }
00910       else {
00911         // Generate C1_ by cloning itself ONLY if more space is needed.
00912         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
00913           Teuchos::RCP<const MV> tmp = C1_;
00914           C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
00915         }
00916       }
00917 
00918       // Generate r_ only if it doesn't exist
00919       if (r_ == Teuchos::null)
00920         r_ = MVT::Clone( *rhsMV, 1 );
00921 
00922       // Size of tau_ will change during computation, so just be sure it starts with appropriate size
00923       tau_.resize(recycledBlocks_+1);
00924 
00925       // Size of work_ will change during computation, so just be sure it starts with appropriate size
00926       work_.resize(recycledBlocks_+1);
00927 
00928       // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
00929       ipiv_.resize(recycledBlocks_+1);
00930 
00931       // Generate H2_ only if it doesn't exist, otherwise resize it.
00932       if (H2_ == Teuchos::null)
00933         H2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
00934       else {
00935         if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
00936           H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
00937       }
00938       H2_->putScalar(zero);
00939 
00940       // Generate R_ only if it doesn't exist, otherwise resize it.
00941       if (R_ == Teuchos::null)
00942         R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
00943       else {
00944         if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
00945           R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
00946       }
00947       R_->putScalar(zero);
00948 
00949       // Generate PP_ only if it doesn't exist, otherwise resize it.
00950       if (PP_ == Teuchos::null)
00951         PP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
00952       else {
00953         if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
00954           PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
00955       }
00956 
00957       // Generate HP_ only if it doesn't exist, otherwise resize it.
00958       if (HP_ == Teuchos::null)
00959         HP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
00960       else {
00961         if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
00962           HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
00963       }
00964 
00965     } // end else
00966 }
00967 
00968   
00969 // solve()
00970 template<class ScalarType, class MV, class OP>
00971 ReturnType GCRODRSolMgr<ScalarType,MV,OP>::solve() {
00972   using Teuchos::RCP;
00973 
00974   // Set the current parameters if they were not set before.
00975   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00976   // then didn't set any parameters using setParameters().
00977   if (!isSet_) { setParameters( params_ ); }
00978 
00979   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00980   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00981   std::vector<int> index(numBlocks_+1);
00982   
00983   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
00984 
00985   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00986 
00987   // Create indices for the linear systems to be solved.
00988   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00989   std::vector<int> currIdx(1);
00990   currIdx[0] = 0;
00991 
00992   // Inform the linear problem of the current linear system to solve.
00993   problem_->setLSIndex( currIdx );
00994 
00995   // Check the number of blocks and change them is necessary.
00996   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
00997   if (numBlocks_ > dim) {
00998     numBlocks_ = dim;
00999     printer_->stream(Warnings) << 
01000       "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
01001       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
01002     params_->set("Num Blocks", numBlocks_);
01003   } 
01004 
01005   // Assume convergence is achieved, then let any failed convergence set this to false.
01006   bool isConverged = true;  
01007 
01008   // Initialize storage for all state variables
01009   initializeStateStorage();
01010 
01012   // Parameter list
01013   Teuchos::ParameterList plist;
01014   
01015   plist.set("Num Blocks",numBlocks_);
01016   plist.set("Recycled Blocks",recycledBlocks_);
01017   
01019   // GCRODR solver
01020   
01021   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
01022   gcrodr_iter = Teuchos::rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01023   // Number of iterations required to generate initial recycle space (if needed)
01024   int prime_iterations = 0;
01025 
01026   // Enter solve() iterations
01027   {
01028     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01029     
01030     while ( numRHS2Solve > 0 ) {
01031       
01032       // Reset the status test.  
01033       outputTest_->reset();
01034       
01036       // Initialize recycled subspace for GCRODR
01037       
01038       // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
01039       if (keff > 0) {
01040   TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
01041          "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
01042 
01043   printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
01044   // Compute image of U_ under the new operator
01045   index.resize(keff);
01046   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01047   RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01048   RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01049   problem_->apply( *Utmp, *Ctmp );
01050 
01051   RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01052 
01053   // Orthogonalize this block
01054   // Get a matrix to hold the orthonormalization coefficients.
01055         Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff ); 
01056   int rank = ortho_->normalize(*Ctmp, Teuchos::rcp(&Rtmp,false));
01057   // Throw an error if we could not orthogonalize this block
01058   TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
01059   
01060   // U_ = U_*R^{-1} 
01061   // First, compute LU factorization of R
01062   int info = 0;
01063         ipiv_.resize(Rtmp.numRows());
01064   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01065   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01066   
01067   // Now, form inv(R)
01068   int lwork = Rtmp.numRows();
01069         work_.resize(lwork);
01070   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01071   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01072   
01073         // U_ = U1_; (via a swap)
01074   MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
01075         std::swap(U_, U1_);
01076 
01077         // Must reinitialize after swap
01078   index.resize(keff);
01079   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01080   Ctmp  = MVT::CloneViewNonConst( *C_, index );
01081   Utmp  = MVT::CloneView( *U_, index );
01082 
01083   // Compute C_'*r_
01084   Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
01085   problem_->computeCurrPrecResVec( &*r_ );
01086   MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
01087 
01088   // Update solution ( x += U_*C_'*r_ )
01089   MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *problem_->getCurrLHSVec() );
01090   
01091   // Update residual norm ( r -= C_*C_'*r_ )  
01092   MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
01093 
01094         // We recycled space from previous call
01095         prime_iterations = 0;
01096 
01097       }
01098       else {
01099   
01100   // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
01101   printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
01102   
01103   Teuchos::ParameterList primeList;
01104   
01105   // Tell the block solver that the block size is one.
01106         primeList.set("Num Blocks",numBlocks_);
01107         primeList.set("Recycled Blocks",0);
01108   
01109   //  Create GCRODR iterator object to perform one cycle of GMRES.
01110   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
01111   gcrodr_prime_iter = Teuchos::rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
01112   
01113   // Create the first block in the current Krylov basis (residual).
01114   problem_->computeCurrPrecResVec( &*r_ );
01115         index.resize( 1 ); index[0] = 0;
01116   RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01117   MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01118 
01119   // Set the new state and initialize the solver.
01120   GCRODRIterState<ScalarType,MV> newstate;
01121         index.resize( numBlocks_+1 );
01122         for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01123         newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01124   newstate.U = Teuchos::null;
01125   newstate.C = Teuchos::null;
01126         newstate.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
01127   newstate.B = Teuchos::null;
01128   newstate.curDim = 0;
01129   gcrodr_prime_iter->initialize(newstate);
01130   
01131   // Perform one cycle of GMRES
01132   bool primeConverged = false;
01133   try {
01134     gcrodr_prime_iter->iterate();
01135 
01136     // Check convergence first
01137     if ( convTest_->getStatus() == Passed ) {
01138       // we have convergence
01139       primeConverged = true;
01140     }
01141   }
01142   catch (const GCRODRIterOrthoFailure &e) {
01143     // Try to recover the most recent least-squares solution
01144     gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
01145     
01146     // Check to see if the most recent least-squares solution yielded convergence.
01147     sTest_->checkStatus( &*gcrodr_prime_iter );
01148     if (convTest_->getStatus() == Passed)
01149       primeConverged = true;
01150   }
01151   catch (const std::exception &e) {
01152     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01153            << gcrodr_prime_iter->getNumIters() << std::endl 
01154            << e.what() << std::endl;
01155     throw;
01156   }
01157         // Record number of iterations in generating initial recycle spacec
01158         prime_iterations = gcrodr_prime_iter->getNumIters();
01159            
01160   // Update the linear problem.
01161   RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
01162   problem_->updateSolution( update, true );
01163 
01164   // Get the state.
01165   newstate = gcrodr_prime_iter->getState();
01166   int p = newstate.curDim;
01167 
01168   // Compute harmonic Ritz vectors 
01169   // NOTE:  The storage for the harmonic Ritz vectors (PP) is made one column larger 
01170         //        just in case we split a complex conjugate pair.
01171   // NOTE:  Generate a recycled subspace only if we have enough vectors.  If we converged
01172   //        too early, move on to the next linear system and try to generate a subspace again.
01173   if (recycledBlocks_ < p+1) {
01174     int info = 0;
01175           RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01176           // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
01177     keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
01178           // Hereafter, only keff columns of PP are needed
01179           PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
01180           // Now get views into C, U, V
01181           index.resize(keff);
01182           for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01183           RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01184           RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01185           RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01186     index.resize(p);
01187           for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
01188     RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01189 
01190     // Form U (the subspace to recycle)
01191           // U = newstate.V(:,1:p) * PP;
01192           MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01193 
01194     // Form orthonormalized C and adjust U so that C = A*U
01195 
01196     // First, compute [Q, R] = qr(H*P);
01197 
01198           // Step #1: Form HP = H*P
01199     Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
01200     Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
01201     HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
01202 
01203     // Step #1.5: Perform workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01204           int lwork = -1;
01205           tau_.resize(keff);
01206           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01207     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01208 
01209           // Step #2: Compute QR factorization of HP
01210           lwork = (int)work_[0];
01211           work_.resize(lwork);
01212           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01213     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01214 
01215           // Step #3: Explicitly construct Q and R factors 
01216     // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01217           Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
01218           for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01219           lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01220     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01221 
01222           // Now we have [Q,R] = qr(H*P)
01223 
01224     // Now compute C = V(:,1:p+1) * Q 
01225     index.resize( p+1 );
01226           for (int ii=0; ii < (p+1); ++ii) { index[ii] = ii; }
01227     Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
01228           MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01229         
01230     // Finally, compute U = U*R^{-1}. 
01231           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as 
01232           // backsolve capabilities don't exist in the Belos::MultiVec class
01233 
01234     // Step #1: First, compute LU factorization of R
01235     ipiv_.resize(Rtmp.numRows());
01236     lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01237     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01238     
01239     // Step #2: Form inv(R)
01240     lwork = Rtmp.numRows();
01241     work_.resize(lwork);
01242     lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01243     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01244     
01245           // Step #3: Let U = U * R^{-1}
01246     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01247 
01248     printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01249 
01250   }  // if (recycledBlocks_ < p+1)
01251 
01252         // Return to outer loop if the priming solve converged, set the next linear system.
01253   if (primeConverged) {
01254     // Inform the linear problem that we are finished with this block linear system.
01255     problem_->setCurrLS();
01256 
01257           // Update indices for the linear systems to be solved.
01258           numRHS2Solve -= 1;
01259           if ( numRHS2Solve > 0 ) {
01260             currIdx[0]++;
01261  
01262             // Set the next indices.
01263             problem_->setLSIndex( currIdx );
01264           }
01265           else {
01266             currIdx.resize( numRHS2Solve );
01267           }
01268 
01269     continue;
01270   }
01271       } // if (keff > 0) ...
01272       
01273       // Prepare for the Gmres iterations with the recycled subspace.
01274 
01275       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01276       gcrodr_iter->setSize( keff, numBlocks_ );
01277       
01278       // Reset the number of iterations.
01279       gcrodr_iter->resetNumIters(prime_iterations);
01280 
01281       // Reset the number of calls that the status test output knows about.
01282       outputTest_->resetNumCalls();
01283 
01284       // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
01285       problem_->computeCurrPrecResVec( &*r_ );
01286       index.resize( 1 ); index[0] = 0;
01287       RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01288       MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01289 
01290       // Set the new state and initialize the solver.
01291       GCRODRIterState<ScalarType,MV> newstate;
01292       index.resize( numBlocks_+1 );
01293       for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01294       newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01295       index.resize( keff );
01296       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01297       newstate.C  = MVT::CloneViewNonConst( *C_,  index );
01298       newstate.U  = MVT::CloneViewNonConst( *U_,  index );
01299       newstate.B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_,         keff, numBlocks_,    0, keff ) );
01300       newstate.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01301       newstate.curDim = 0;
01302       gcrodr_iter->initialize(newstate);
01303 
01304       // variables needed for inner loop
01305       int numRestarts = 0;
01306       std::vector<MagnitudeType> d(keff);
01307       RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp;
01308       while(1) {
01309   
01310   // tell gcrodr_iter to iterate
01311   try {
01312     gcrodr_iter->iterate();
01313     
01315     //
01316     // check convergence first
01317     //
01319     if ( convTest_->getStatus() == Passed ) {
01320       // we have convergence
01321       break;  // break from while(1){gcrodr_iter->iterate()}
01322     }
01324     //
01325     // check for maximum iterations
01326     //
01328     else if ( maxIterTest_->getStatus() == Passed ) {
01329       // we don't have convergence
01330       isConverged = false;
01331       break;  // break from while(1){gcrodr_iter->iterate()}
01332     }
01334     //
01335     // check for restarting, i.e. the subspace is full
01336     //
01338     else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
01339 
01340             // Update the recycled subspace even if we have hit the maximum number of restarts.
01341  
01342       // Get the state.
01343       GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
01344             int p = oldState.curDim;     
01345 
01346             // Update the linear problem.
01347             RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01348             problem_->updateSolution( update, true );
01349 
01350             // Take the norm of the recycled vectors
01351             {
01352               index.resize(keff);
01353               for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01354               RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01355               d.resize(keff);
01356               MVT::MvNorm( *Utmp, d );
01357               for (int i=0; i<keff; ++i) {
01358                 d[i] = one / d[i];
01359               }
01360               MVT::MvScale( *Utmp, d );
01361             }
01362 
01363             // Get view into current "full" upper Hessnburg matrix
01364             H2tmp = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
01365 
01366             // Insert D into the leading keff x keff  block of H2
01367             for (int i=0; i<keff; ++i) {
01368               (*H2tmp)(i,i) = d[i];
01369             }
01370             
01371             // Compute the harmoic Ritz pairs for the generalized eigenproblem
01372             // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
01373       int keff_new;
01374             {
01375               Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
01376               keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, PPtmp );
01377             }
01378 
01379       printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff_new << std::endl << std::endl;
01380 
01381       // Code to form new U, C
01382       // U = [U V(:,1:p)] * P; (in two steps)
01383 
01384       // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01385             RCP<MV> U1tmp;
01386             {
01387               index.resize( keff );
01388               for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01389               RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
01390               index.resize( keff_new );
01391               for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01392               U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
01393               Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
01394               MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
01395             }
01396 
01397             // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01398             {
01399               index.resize(p);
01400               for (int ii=0; ii < p; ii++) { index[ii] = ii; }
01401               RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01402               Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
01403               MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
01404             }
01405 
01406             // Form HP = H*P
01407             Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
01408             {
01409               Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
01410               HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
01411             }
01412 
01413       // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01414       int info = 0, lwork = -1;
01415       tau_.resize(keff_new);
01416       lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01417       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01418       
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,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01423 
01424       // Explicitly construct Q and R factors 
01425       // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01426             Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
01427             for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
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,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01430 
01431       // Form orthonormalized C and adjust U accordingly so that C = A*U      
01432       // C = [C V] * Q;
01433 
01434             // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01435             {
01436               RCP<MV> C1tmp;
01437               {
01438                 index.resize(keff);
01439                 for (int i=0; i < keff; i++) { index[i] = i; }
01440                 RCP<const MV> Ctmp  = MVT::CloneView( *C_,  index );
01441                 index.resize(keff_new);
01442                 for (int i=0; i < keff_new; i++) { index[i] = i; }
01443                 C1tmp  = MVT::CloneViewNonConst( *C1_,  index );
01444                 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
01445                 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
01446               }
01447 
01448               // Now compute C += V(:,1:p+1) * Q 
01449               {
01450                 index.resize( p+1 );
01451                 for (int i=0; i < p+1; ++i) { index[i] = i; }
01452                 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01453                 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
01454                 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
01455               }
01456             }
01457 
01458             // C_ = C1_; (via a swap)
01459             std::swap(C_, C1_);
01460 
01461       // Finally, compute U_ = U_*R^{-1}  
01462       // First, compute LU factorization of R
01463             ipiv_.resize(Rtmp.numRows());
01464       lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01465       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01466       
01467       // Now, form inv(R)
01468       lwork = Rtmp.numRows();
01469       work_.resize(lwork);
01470       lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01471       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01472     
01473             {
01474               index.resize(keff_new);
01475               for (int i=0; i < keff_new; i++) { index[i] = i; }
01476               RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_,  index );
01477               MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01478             }
01479 
01480             // NOTE:  If we have hit the maximum number of restarts then quit
01481       if ( numRestarts >= maxRestarts_ ) {
01482         isConverged = false;
01483         break; // break from while(1){gcrodr_iter->iterate()}
01484       }
01485       numRestarts++;
01486       
01487       printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01488 
01489             // Create the restart vector (first block in the current Krylov basis)
01490             problem_->computeCurrPrecResVec( &*r_ );
01491             index.resize( 1 ); index[0] = 0;
01492             RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01493             MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01494       
01495       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01496             if (keff != keff_new) {
01497           keff = keff_new;
01498         gcrodr_iter->setSize( keff, numBlocks_ );
01499               Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01500               b1.putScalar(zero);
01501             }
01502 
01503       // Set the new state and initialize the solver.
01504       GCRODRIterState<ScalarType,MV> restartState;
01505             index.resize( numBlocks_+1 );
01506             for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01507             restartState.V  = MVT::CloneViewNonConst( *V_,  index );
01508             index.resize( keff );
01509             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01510             restartState.U  = MVT::CloneViewNonConst( *U_,  index );
01511             restartState.C  = MVT::CloneViewNonConst( *C_,  index );
01512             restartState.B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff,         numBlocks_,    0, keff ) );
01513             restartState.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01514             restartState.curDim = 0;
01515             gcrodr_iter->initialize(restartState);
01516 
01517     } // end of restarting
01518     
01520     //
01521     // we returned from iterate(), but none of our status tests Passed.
01522     // something is wrong, and it is probably our fault.
01523     //
01525     
01526     else {
01527       TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::GCRODRSolMgr::solve(): Invalid return from GCRODRIter::iterate().");
01528     }
01529   }
01530         catch (const GCRODRIterOrthoFailure &e) {
01531     // Try to recover the most recent least-squares solution
01532     gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
01533     
01534     // Check to see if the most recent least-squares solution yielded convergence.
01535     sTest_->checkStatus( &*gcrodr_iter );
01536     if (convTest_->getStatus() != Passed)
01537       isConverged = false;
01538     break;
01539         }
01540         catch (const std::exception &e) {
01541     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01542                              << gcrodr_iter->getNumIters() << std::endl 
01543            << e.what() << std::endl;
01544           throw;
01545   }
01546       }
01547      
01548       // Compute the current solution.
01549       // Update the linear problem.
01550       RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01551       problem_->updateSolution( update, true );
01552 
01553       // Inform the linear problem that we are finished with this block linear system.
01554       problem_->setCurrLS();
01555 
01556       // Update indices for the linear systems to be solved.
01557       numRHS2Solve -= 1;
01558       if ( numRHS2Solve > 0 ) {
01559   currIdx[0]++;
01560 
01561         // Set the next indices.
01562         problem_->setLSIndex( currIdx );
01563       }
01564       else {
01565         currIdx.resize( numRHS2Solve );
01566       }
01567  
01568     }// while ( numRHS2Solve > 0 )
01569     
01570   }
01571   
01572   // print final summary
01573   sTest_->print( printer_->stream(FinalSummary) );
01574   
01575   // print timing information
01576   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01577  
01578   // get iteration information for this solve
01579   numIters_ = maxIterTest_->getNumIters();
01580  
01581   if (!isConverged) {
01582     return Unconverged; // return from GCRODRSolMgr::solve() 
01583   }
01584   return Converged; // return from GCRODRSolMgr::solve() 
01585 }
01586 
01587 //  Compute the harmonic eigenpairs of the projected, dense system.
01588 template<class ScalarType, class MV, class OP>
01589 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs1(int m, 
01590                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01591                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01592   int i, j;
01593   bool xtraVec = false;
01594   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01595   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01596 
01597   // Real and imaginary eigenvalue components
01598   std::vector<MagnitudeType> wr(m), wi(m);
01599 
01600   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01601   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,false);
01602 
01603   // Magnitude of harmonic Ritz values
01604   std::vector<MagnitudeType> w(m);
01605 
01606   // Sorted order of harmonic Ritz values, also used for DGEEV
01607   std::vector<int> iperm(m);
01608 
01609   // Size of workspace and workspace for DGEEV
01610   int lwork = 4*m;
01611   std::vector<ScalarType> work(lwork);
01612 
01613   // Output info
01614   int info = 0;
01615 
01616   // Solve linear system:  H_m^{-H}*e_m 
01617   Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
01618   Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
01619   e_m[m-1] = one;
01620   lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
01621   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01622 
01623   // Compute H_m + d*H_m^{-H}*e_m*e_m^H
01624   ScalarType d = HH(m, m-1) * HH(m, m-1);
01625   Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( HH );
01626   for( i=0; i<m; ++i )
01627     harmHH(i, m-1) += d * e_m[i];
01628 
01629   // Revise to do query for optimal workspace first
01630   // Create simple storage for the left eigenvectors, which we don't care about.
01631   const int ldvl = m;
01632   ScalarType* vl = 0;
01633   lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
01634         vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01635   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01636 
01637   // Construct magnitude of each harmonic Ritz value
01638   for( i=0; i<m; ++i )
01639     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01640 
01641   // Construct magnitude of each harmonic Ritz value
01642   this->sort(w, m, iperm);
01643 
01644   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01645   if (wi[iperm[recycledBlocks_-1]] != zero) {
01646     int countImag = 0;
01647     for ( i=0; i<recycledBlocks_; ++i ) {
01648       if (wi[iperm[i]] != zero)
01649   countImag++;
01650     }
01651     // Check to see if this count is even or odd:
01652     if (countImag % 2)
01653       xtraVec = true;
01654   }
01655 
01656   // Select recycledBlocks_ smallest eigenvectors
01657   for( i=0; i<recycledBlocks_; ++i ) {
01658     for( j=0; j<m; j++ ) {
01659       PP(j,i) = vr(j,iperm[i]);
01660     }
01661   }
01662   
01663   if (xtraVec) { // we need to store one more vector
01664     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
01665       for( j=0; j<m; ++j ) {   // so get the "imag" component
01666   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
01667       }
01668     }
01669     else { //  I picked the "imag" component
01670       for( j=0; j<m; ++j ) {   // so get the "real" component
01671   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
01672       }
01673     }
01674   }
01675 
01676   // Return whether we needed to store an additional vector
01677   if (xtraVec) {
01678     return recycledBlocks_+1;
01679   }
01680   return recycledBlocks_;
01681 }
01682 
01683 //  Compute the harmonic eigenpairs of the projected, dense system.
01684 template<class ScalarType, class MV, class OP>
01685 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs2(int keff, int m, 
01686                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01687                  const Teuchos::RCP<const MV>& VV,
01688                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01689   int i, j;
01690   int m2 = HH.numCols();
01691   bool xtraVec = false;
01692   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01693   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01694   std::vector<int> index;
01695   
01696   // Real and imaginary eigenvalue components
01697   std::vector<ScalarType> wr(m2), wi(m2);
01698 
01699   // Magnitude of harmonic Ritz values
01700   std::vector<MagnitudeType> w(m2);
01701 
01702   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01703   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,false);
01704 
01705   // Sorted order of harmonic Ritz values
01706   std::vector<int> iperm(m2);
01707 
01708   // Form matrices for generalized eigenproblem
01709   
01710   // B = H2' * H2; Don't zero out matrix when constructing
01711   Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,false);
01712   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
01713   
01714   // A_tmp = | C'*U        0 |
01715   //         | V_{m+1}'*U  I |
01716   Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keff+m+1, keff+m );
01717 
01718   // A_tmp(1:keff,1:keff) = C' * U;
01719   index.resize(keff);
01720   for (int i=0; i<keff; ++i) { index[i] = i; }
01721   Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_, index );
01722   Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01723   Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keff, keff );
01724   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
01725 
01726   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
01727   Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keff, keff );
01728   index.resize(m+1);
01729   for (i=0; i < m+1; i++) { index[i] = i; }
01730   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
01731   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
01732 
01733   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
01734   for( i=keff; i<keff+m; i++ ) {
01735     A_tmp(i,i) = one;
01736   }
01737 
01738   // A = H2' * A_tmp;
01739   Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
01740   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
01741 
01742   // Compute k smallest harmonic Ritz pairs
01743   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
01744   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
01745   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
01746   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
01747   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
01748   char balanc='P', jobvl='N', jobvr='V', sense='N'; 
01749   int ld = A.numRows();
01750   int lwork = 6*ld;
01751   int ldvl = ld, ldvr = ld;
01752   int info = 0,ilo = 0,ihi = 0;
01753   ScalarType abnrm = zero, bbnrm = zero;
01754   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
01755   std::vector<ScalarType> beta(ld);
01756   std::vector<ScalarType> work(lwork);
01757   std::vector<MagnitudeType> lscale(ld), rscale(ld);  
01758   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
01759   std::vector<int> iwork(ld+6);
01760   int *bwork = 0; // If sense == 'N', bwork is never referenced
01761   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0], 
01762                &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0], 
01763                &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
01764   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
01765   
01766   // Construct magnitude of each harmonic Ritz value
01767   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
01768   for( i=0; i<ld; i++ ) {
01769     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
01770   }
01771 
01772   // Construct magnitude of each harmonic Ritz value
01773   this->sort(w,ld,iperm);
01774 
01775   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01776   if (wi[iperm[ld-recycledBlocks_]] != zero) {
01777     int countImag = 0;
01778     for ( i=ld-recycledBlocks_; i<ld; i++ ) {
01779       if (wi[iperm[i]] != zero)
01780   countImag++;
01781     }
01782     // Check to see if this count is even or odd:
01783     if (countImag % 2)
01784       xtraVec = true;
01785   }
01786   
01787   // Select recycledBlocks_ smallest eigenvectors
01788   for( i=0; i<recycledBlocks_; i++ ) {
01789     for( j=0; j<ld; j++ ) {
01790       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
01791     }
01792   }
01793   
01794   if (xtraVec) { // we need to store one more vector
01795     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
01796       for( j=0; j<ld; j++ ) {   // so get the "imag" component
01797   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
01798       }
01799     }
01800     else { // I picked the "imag" component
01801       for( j=0; j<ld; j++ ) {   // so get the "real" component
01802   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
01803       }
01804     }
01805   }
01806   
01807   // Return whether we needed to store an additional vector
01808   if (xtraVec) {
01809     return recycledBlocks_+1;
01810   }
01811   return recycledBlocks_;
01812 }
01813 
01814 
01815 // This method sorts list of n floating-point numbers and return permutation vector
01816 template<class ScalarType, class MV, class OP>
01817 void GCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
01818   int l, r, j, i, flag;
01819   int    RR2;
01820   double dRR, dK;
01821   
01822   // Initialize the permutation vector.
01823   for(j=0;j<n;j++)
01824     iperm[j] = j;
01825   
01826   if (n <= 1) return;
01827   
01828   l    = n / 2 + 1;
01829   r    = n - 1;
01830   l    = l - 1;
01831   dRR  = dlist[l - 1];
01832   dK   = dlist[l - 1];
01833   
01834   RR2 = iperm[l - 1];
01835   while (r != 0) {
01836     j = l;
01837     flag = 1;
01838     
01839     while (flag == 1) {
01840       i = j;
01841       j = j + j;
01842       
01843       if (j > r + 1)
01844   flag = 0;
01845       else {
01846   if (j < r + 1)
01847     if (dlist[j] > dlist[j - 1]) j = j + 1;
01848   
01849   if (dlist[j - 1] > dK) {
01850     dlist[i - 1] = dlist[j - 1];
01851     iperm[i - 1] = iperm[j - 1];
01852   }
01853   else {
01854     flag = 0;
01855   }
01856       }
01857     }
01858     dlist[i - 1] = dRR;
01859     iperm[i - 1] = RR2;
01860     
01861     if (l == 1) {
01862       dRR  = dlist [r];
01863       RR2 = iperm[r];
01864       dK = dlist[r];
01865       dlist[r] = dlist[0];
01866       iperm[r] = iperm[0];
01867       r = r - 1;
01868     }
01869     else {
01870       l   = l - 1;
01871       dRR  = dlist[l - 1];
01872       RR2  = iperm[l - 1];
01873       dK   = dlist[l - 1];
01874     }
01875   }
01876   dlist[0] = dRR;
01877   iperm[0] = RR2; 
01878 }
01879 
01880 
01881 //  This method requires the solver manager to return a string that describes itself.
01882 template<class ScalarType, class MV, class OP>
01883 std::string GCRODRSolMgr<ScalarType,MV,OP>::description() const {
01884   std::ostringstream oss;
01885   oss << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01886   oss << "{";
01887   oss << "Ortho Type='"<<orthoType_;
01888   oss << ", Num Blocks=" <<numBlocks_;
01889   oss << ", Num Recycle Blocks=" << recycledBlocks_;
01890   oss << ", Max Restarts=" << maxRestarts_;
01891   oss << "}";
01892   return oss.str();
01893 }
01894   
01895 } // end Belos namespace
01896 
01897 #endif /* BELOS_GCRODR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:14 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3