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 
00973   // Set the current parameters if they were not set before.
00974   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00975   // then didn't set any parameters using setParameters().
00976   if (!isSet_) { setParameters( params_ ); }
00977 
00978   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00979   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00980   std::vector<int> index(numBlocks_+1);
00981   
00982   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
00983 
00984   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00985 
00986   // Create indices for the linear systems to be solved.
00987   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00988   std::vector<int> currIdx(1);
00989   currIdx[0] = 0;
00990 
00991   // Inform the linear problem of the current linear system to solve.
00992   problem_->setLSIndex( currIdx );
00993 
00994   // Check the number of blocks and change them is necessary.
00995   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
00996   if (numBlocks_ > dim) {
00997     numBlocks_ = dim;
00998     printer_->stream(Warnings) << 
00999       "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
01000       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
01001     params_->set("Num Blocks", numBlocks_);
01002   } 
01003 
01004   // Assume convergence is achieved, then let any failed convergence set this to false.
01005   bool isConverged = true;  
01006 
01007   // Initialize storage for all state variables
01008   initializeStateStorage();
01009 
01011   // Parameter list
01012   Teuchos::ParameterList plist;
01013   
01014   plist.set("Num Blocks",numBlocks_);
01015   plist.set("Recycled Blocks",recycledBlocks_);
01016   
01018   // GCRODR solver
01019   
01020   Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
01021   gcrodr_iter = Teuchos::rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01022   // Number of iterations required to generate initial recycle space (if needed)
01023   int prime_iterations = 0;
01024 
01025   // Enter solve() iterations
01026   {
01027     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01028     
01029     while ( numRHS2Solve > 0 ) {
01030       
01031       // Reset the status test.  
01032       outputTest_->reset();
01033       
01035       // Initialize recycled subspace for GCRODR
01036       
01037       // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
01038       if (keff > 0) {
01039   TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
01040          "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
01041 
01042   printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
01043   // Compute image of U_ under the new operator
01044   index.resize(keff);
01045   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01046   Teuchos::RCP<MV> Ctmp  = MVT::CloneView( *C_, index );
01047   Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_, index );
01048   Teuchos::RCP<MV> U1tmp = MVT::CloneView( *U1_, index );
01049   problem_->apply( *Utmp, *Ctmp );
01050 
01051   // Orthogonalize this block
01052   // Get a matrix to hold the orthonormalization coefficients.
01053         Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff ); 
01054   int rank = ortho_->normalize(*Ctmp, Teuchos::rcp(&Rtmp,false));
01055   // Throw an error if we could not orthogonalize this block
01056   TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
01057   
01058   // U_ = U_*R^{-1} 
01059   // First, compute LU factorization of R
01060   int info = 0;
01061         ipiv_.resize(Rtmp.numRows());
01062   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01063   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01064   
01065   // Now, form inv(R)
01066   int lwork = Rtmp.numRows();
01067         work_.resize(lwork);
01068   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01069   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01070   
01071         // U_ = U1_; (via a swap)
01072   MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
01073         Teuchos::RCP<MV> Uswap;
01074         Uswap = U_;
01075   U_ = U1_;
01076         U1_ = Uswap;
01077         Uswap = Teuchos::null;
01078 
01079         // Must reinitialize after swap
01080   index.resize(keff);
01081   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01082   Ctmp  = MVT::CloneView( *C_, index );
01083   Utmp  = MVT::CloneView( *U_, index );
01084 
01085   // Compute C_'*r_
01086   Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
01087   problem_->computeCurrPrecResVec( &*r_ );
01088   MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
01089 
01090   // Update solution ( x += U_*C_'*r_ )
01091   MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *problem_->getCurrLHSVec() );
01092   
01093   // Update residual norm ( r -= C_*C_'*r_ )  
01094   MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
01095 
01096         // We recycled space from previous call
01097         prime_iterations = 0;
01098 
01099       }
01100       else {
01101   
01102   // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
01103   printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
01104   
01105   Teuchos::ParameterList primeList;
01106   
01107   // Tell the block solver that the block size is one.
01108         primeList.set("Num Blocks",numBlocks_);
01109         primeList.set("Recycled Blocks",0);
01110   
01111   //  Create GCRODR iterator object to perform one cycle of GMRES.
01112   Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
01113   gcrodr_prime_iter = Teuchos::rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
01114   
01115   // Create the first block in the current Krylov basis (residual).
01116   problem_->computeCurrPrecResVec( &*r_ );
01117         index.resize( 1 ); index[0] = 0;
01118   Teuchos::RCP<MV> v0 =  MVT::CloneView( *V_,  index );
01119   MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01120 
01121   // Set the new state and initialize the solver.
01122   GCRODRIterState<ScalarType,MV> newstate;
01123         index.resize( numBlocks_+1 );
01124         for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01125         newstate.V  = MVT::CloneView( *V_,  index );
01126   newstate.U = Teuchos::null;
01127   newstate.C = Teuchos::null;
01128         newstate.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
01129   newstate.B = Teuchos::null;
01130   newstate.curDim = 0;
01131   gcrodr_prime_iter->initialize(newstate);
01132   
01133   // Perform one cycle of GMRES
01134   bool primeConverged = false;
01135   try {
01136     gcrodr_prime_iter->iterate();
01137 
01138     // Check convergence first
01139     if ( convTest_->getStatus() == Passed ) {
01140       // we have convergence
01141       primeConverged = true;
01142     }
01143   }
01144   catch (const GCRODRIterOrthoFailure &e) {
01145     // Try to recover the most recent least-squares solution
01146     gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
01147     
01148     // Check to see if the most recent least-squares solution yielded convergence.
01149     sTest_->checkStatus( &*gcrodr_prime_iter );
01150     if (convTest_->getStatus() == Passed)
01151       primeConverged = true;
01152   }
01153   catch (const std::exception &e) {
01154     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01155            << gcrodr_prime_iter->getNumIters() << std::endl 
01156            << e.what() << std::endl;
01157     throw;
01158   }
01159         // Record number of iterations in generating initial recycle spacec
01160         prime_iterations = gcrodr_prime_iter->getNumIters();
01161            
01162   // Update the linear problem.
01163   Teuchos::RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
01164   problem_->updateSolution( update, true );
01165 
01166   // Get the state.
01167   newstate = gcrodr_prime_iter->getState();
01168   int p = newstate.curDim;
01169 
01170   // Compute harmonic Ritz vectors 
01171   // NOTE:  The storage for the harmonic Ritz vectors (PP) is made one column larger 
01172         //        just in case we split a complex conjugate pair.
01173   // NOTE:  Generate a recycled subspace only if we have enough vectors.  If we converged
01174   //        too early, move on to the next linear system and try to generate a subspace again.
01175   if (recycledBlocks_ < p+1) {
01176     int info = 0;
01177           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01178           // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
01179     keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
01180           // Hereafter, only keff columns of PP are needed
01181           PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
01182           // Now get views into C, U, V
01183           index.resize(keff);
01184           for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01185           Teuchos::RCP<MV> Ctmp  = MVT::CloneView( *C_, index );
01186           Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_, index );
01187           Teuchos::RCP<MV> U1tmp = MVT::CloneView( *U1_, index );
01188     index.resize(p);
01189           for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
01190     Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01191 
01192     // Form U (the subspace to recycle)
01193           // U = newstate.V(:,1:p) * PP;
01194           MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01195 
01196     // Form orthonormalized C and adjust U so that C = A*U
01197 
01198     // First, compute [Q, R] = qr(H*P);
01199 
01200           // Step #1: Form HP = H*P
01201     Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
01202     Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
01203     HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
01204 
01205     // Step #1.5: Perform workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01206           int lwork = -1;
01207           tau_.resize(keff);
01208           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01209     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01210 
01211           // Step #2: Compute QR factorization of HP
01212           lwork = (int)work_[0];
01213           work_.resize(lwork);
01214           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01215     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01216 
01217           // Step #3: Explicitly construct Q and R factors 
01218     // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01219           Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
01220           for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01221           lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01222     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01223 
01224           // Now we have [Q,R] = qr(H*P)
01225 
01226     // Now compute C = V(:,1:p+1) * Q 
01227     index.resize( p+1 );
01228           for (int ii=0; ii < (p+1); ++ii) { index[ii] = ii; }
01229     Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
01230           MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01231         
01232     // Finally, compute U = U*R^{-1}. 
01233           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as 
01234           // backsolve capabilities don't exist in the Belos::MultiVec class
01235 
01236     // Step #1: First, compute LU factorization of R
01237     ipiv_.resize(Rtmp.numRows());
01238     lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01239     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01240     
01241     // Step #2: Form inv(R)
01242     lwork = Rtmp.numRows();
01243     work_.resize(lwork);
01244     lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01245     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01246     
01247           // Step #3: Let U = U * R^{-1}
01248     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01249 
01250     printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01251 
01252   }  // if (recycledBlocks_ < p+1)
01253 
01254         // Return to outer loop if the priming solve converged, set the next linear system.
01255   if (primeConverged) {
01256     // Inform the linear problem that we are finished with this block linear system.
01257     problem_->setCurrLS();
01258 
01259           // Update indices for the linear systems to be solved.
01260           numRHS2Solve -= 1;
01261           if ( numRHS2Solve > 0 ) {
01262             currIdx[0]++;
01263  
01264             // Set the next indices.
01265             problem_->setLSIndex( currIdx );
01266           }
01267           else {
01268             currIdx.resize( numRHS2Solve );
01269           }
01270 
01271     continue;
01272   }
01273       } // if (keff > 0) ...
01274       
01275       // Prepare for the Gmres iterations with the recycled subspace.
01276 
01277       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01278       gcrodr_iter->setSize( keff, numBlocks_ );
01279       
01280       // Reset the number of iterations.
01281       gcrodr_iter->resetNumIters(prime_iterations);
01282 
01283       // Reset the number of calls that the status test output knows about.
01284       outputTest_->resetNumCalls();
01285 
01286       // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
01287       problem_->computeCurrPrecResVec( &*r_ );
01288       index.resize( 1 ); index[0] = 0;
01289       Teuchos::RCP<MV> v0 =  MVT::CloneView( *V_,  index );
01290       MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01291 
01292       // Set the new state and initialize the solver.
01293       GCRODRIterState<ScalarType,MV> newstate;
01294       index.resize( numBlocks_+1 );
01295       for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01296       newstate.V  = MVT::CloneView( *V_,  index );
01297       index.resize( keff );
01298       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01299       newstate.C  = MVT::CloneView( *C_,  index );
01300       newstate.U  = MVT::CloneView( *U_,  index );
01301       newstate.B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_,         keff, numBlocks_,    0, keff ) );
01302       newstate.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01303       newstate.curDim = 0;
01304       gcrodr_iter->initialize(newstate);
01305 
01306       // variables needed for inner loop
01307       int numRestarts = 0;
01308       std::vector<MagnitudeType> d(keff);
01309       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp;
01310       Teuchos::RCP<MV> Utmp;
01311       Teuchos::RCP<MV> Ctmp;
01312       Teuchos::RCP<MV> U1tmp;
01313       Teuchos::RCP<MV> C1tmp;
01314       Teuchos::RCP<MV> Vtmp;
01315       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp;
01316       while(1) {
01317   
01318   // tell gcrodr_iter to iterate
01319   try {
01320     gcrodr_iter->iterate();
01321     
01323     //
01324     // check convergence first
01325     //
01327     if ( convTest_->getStatus() == Passed ) {
01328       // we have convergence
01329       break;  // break from while(1){gcrodr_iter->iterate()}
01330     }
01332     //
01333     // check for maximum iterations
01334     //
01336     else if ( maxIterTest_->getStatus() == Passed ) {
01337       // we don't have convergence
01338       isConverged = false;
01339       break;  // break from while(1){gcrodr_iter->iterate()}
01340     }
01342     //
01343     // check for restarting, i.e. the subspace is full
01344     //
01346     else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
01347 
01348             // Update the recycled subspace even if we have hit the maximum number of restarts.
01349  
01350       // Get the state.
01351       GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
01352             int p = oldState.curDim;     
01353 
01354             // Update the linear problem.
01355             Teuchos::RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01356             problem_->updateSolution( update, true );
01357 
01358             // Take the norm of the recycled vectors
01359             index.resize(keff);
01360             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01361             Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_, index );
01362             d.resize(keff);
01363             MVT::MvNorm( *Utmp, d );
01364             for (int i=0; i<keff; ++i) {
01365         d[i] = one / d[i];
01366       }
01367       MVT::MvScale( *Utmp, d );
01368           
01369             // Get view into current "full" upper Hessnburg matrix
01370             H2tmp = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
01371 
01372             // Insert D into the leading keff x keff  block of H2
01373             for (int i=0; i<keff; ++i) {
01374               (*H2tmp)(i,i) = d[i];
01375             }
01376             
01377             // Compute the harmoic Ritz pairs for the generalized eigenproblem
01378             // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
01379             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 ) );
01380       int keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, *PPtmp );
01381 
01382       printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff_new << std::endl << std::endl;
01383 
01384       // Code to form new U, C
01385       // U = [U V(:,1:p)] * P; (in two steps)
01386 
01387       // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01388             index.resize( keff );
01389             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01390             Utmp  = MVT::CloneView( *U_,  index );
01391             index.resize( keff_new );
01392             for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01393             U1tmp  = MVT::CloneView( *U1_,  index );
01394             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, keff, keff_new ) );
01395       MVT::MvTimesMatAddMv( one, *Utmp, *PPtmp, zero, *U1tmp );
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       index.resize(p);
01399       for (int ii=0; ii < p; ii++) { index[ii] = ii; }
01400       Vtmp = MVT::CloneView( *V_, index );
01401             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff_new, keff ) );
01402       MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, one, *U1tmp );
01403      
01404       // Form HP = H*P
01405             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p+keff, keff_new ) );
01406             Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
01407       HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,*PPtmp,zero);
01408 
01409       // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01410       int info = 0, lwork = -1;
01411       tau_.resize(keff_new);
01412       lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01413       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01414       
01415       lwork = (int)work_[0];
01416       work_.resize(lwork);
01417       lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01418       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01419 
01420       // Explicitly construct Q and R factors 
01421       // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01422             Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
01423             for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
01424       lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01425       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01426 
01427       // Form orthonormalized C and adjust U accordingly so that C = A*U      
01428       // C = [C V] * Q;
01429 
01430       // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01431       index.resize(keff);
01432       for (int i=0; i < keff; i++) { index[i] = i; }
01433             Ctmp  = MVT::CloneView( *C_,  index );
01434       index.resize(keff_new);
01435       for (int i=0; i < keff_new; i++) { index[i] = i; }
01436             C1tmp  = MVT::CloneView( *C1_,  index );
01437             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *HP_, keff, keff_new ) );
01438       MVT::MvTimesMatAddMv( one, *Ctmp, *PPtmp, zero, *C1tmp );
01439 
01440       // Now compute C += V(:,1:p+1) * Q 
01441       index.resize( p+1 );
01442       for (int i=0; i < p+1; ++i) { index[i] = i; }
01443       Vtmp = MVT::CloneView( *V_, index );
01444             PPtmp = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *HP_, p+1, keff_new, keff, 0 ) );
01445       MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, one, *C1tmp );
01446 
01447             // C_ = C1_; (via a swap)
01448             Teuchos::RCP<MV> Cswap;
01449             Cswap = C_;
01450             C_ = C1_;
01451             C1_ = Cswap;
01452             Cswap = Teuchos::null;
01453 
01454       // Finally, compute U_ = U_*R^{-1}  
01455       // First, compute LU factorization of R
01456             ipiv_.resize(Rtmp.numRows());
01457       lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01458       TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01459       
01460       // Now, form inv(R)
01461       lwork = Rtmp.numRows();
01462       work_.resize(lwork);
01463       lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01464       TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01465     
01466       index.resize(keff_new);
01467       for (int i=0; i < keff_new; i++) { index[i] = i; }
01468             Utmp  = MVT::CloneView( *U_,  index );
01469       MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01470 
01471             // NOTE:  If we have hit the maximum number of restarts then quit
01472       if ( numRestarts >= maxRestarts_ ) {
01473         isConverged = false;
01474         break; // break from while(1){gcrodr_iter->iterate()}
01475       }
01476       numRestarts++;
01477       
01478       printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01479 
01480             // Create the restart vector (first block in the current Krylov basis)
01481             problem_->computeCurrPrecResVec( &*r_ );
01482             index.resize( 1 ); index[0] = 0;
01483             Teuchos::RCP<MV> v0 =  MVT::CloneView( *V_,  index );
01484             MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01485       
01486       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01487             if (keff != keff_new) {
01488           keff = keff_new;
01489         gcrodr_iter->setSize( keff, numBlocks_ );
01490               Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01491               b1.putScalar(zero);
01492             }
01493 
01494       // Set the new state and initialize the solver.
01495       GCRODRIterState<ScalarType,MV> restartState;
01496             index.resize( numBlocks_+1 );
01497             for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01498             restartState.V  = MVT::CloneView( *V_,  index );
01499             index.resize( keff );
01500             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01501             restartState.U  = MVT::CloneView( *U_,  index );
01502             restartState.C  = MVT::CloneView( *C_,  index );
01503             restartState.B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff,         numBlocks_,    0, keff ) );
01504             restartState.H = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01505             restartState.curDim = 0;
01506             gcrodr_iter->initialize(restartState);
01507 
01508     } // end of restarting
01509     
01511     //
01512     // we returned from iterate(), but none of our status tests Passed.
01513     // something is wrong, and it is probably our fault.
01514     //
01516     
01517     else {
01518       TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::GCRODRSolMgr::solve(): Invalid return from GCRODRIter::iterate().");
01519     }
01520   }
01521         catch (const GCRODRIterOrthoFailure &e) {
01522     // Try to recover the most recent least-squares solution
01523     gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
01524     
01525     // Check to see if the most recent least-squares solution yielded convergence.
01526     sTest_->checkStatus( &*gcrodr_iter );
01527     if (convTest_->getStatus() != Passed)
01528       isConverged = false;
01529     break;
01530         }
01531         catch (const std::exception &e) {
01532     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01533                              << gcrodr_iter->getNumIters() << std::endl 
01534            << e.what() << std::endl;
01535           throw;
01536   }
01537       }
01538      
01539       // Compute the current solution.
01540       // Update the linear problem.
01541       Teuchos::RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01542       problem_->updateSolution( update, true );
01543 
01544       // Inform the linear problem that we are finished with this block linear system.
01545       problem_->setCurrLS();
01546 
01547       // Update indices for the linear systems to be solved.
01548       numRHS2Solve -= 1;
01549       if ( numRHS2Solve > 0 ) {
01550   currIdx[0]++;
01551 
01552         // Set the next indices.
01553         problem_->setLSIndex( currIdx );
01554       }
01555       else {
01556         currIdx.resize( numRHS2Solve );
01557       }
01558  
01559     }// while ( numRHS2Solve > 0 )
01560     
01561   }
01562   
01563   // print final summary
01564   sTest_->print( printer_->stream(FinalSummary) );
01565   
01566   // print timing information
01567   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01568  
01569   // get iteration information for this solve
01570   numIters_ = maxIterTest_->getNumIters();
01571  
01572   if (!isConverged) {
01573     return Unconverged; // return from GCRODRSolMgr::solve() 
01574   }
01575   return Converged; // return from GCRODRSolMgr::solve() 
01576 }
01577 
01578 //  Compute the harmonic eigenpairs of the projected, dense system.
01579 template<class ScalarType, class MV, class OP>
01580 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs1(int m, 
01581                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01582                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01583   int i, j;
01584   bool xtraVec = false;
01585   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01586   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01587 
01588   // Real and imaginary eigenvalue components
01589   std::vector<MagnitudeType> wr(m), wi(m);
01590 
01591   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01592   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,false);
01593 
01594   // Magnitude of harmonic Ritz values
01595   std::vector<MagnitudeType> w(m);
01596 
01597   // Sorted order of harmonic Ritz values, also used for DGEEV
01598   std::vector<int> iperm(m);
01599 
01600   // Size of workspace and workspace for DGEEV
01601   int lwork = 4*m;
01602   std::vector<ScalarType> work(lwork);
01603 
01604   // Output info
01605   int info = 0;
01606 
01607   // Solve linear system:  H_m^{-H}*e_m 
01608   Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
01609   Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
01610   e_m[m-1] = one;
01611   lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
01612   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01613 
01614   // Compute H_m + d*H_m^{-H}*e_m*e_m^H
01615   ScalarType d = HH(m, m-1) * HH(m, m-1);
01616   Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( HH );
01617   for( i=0; i<m; ++i )
01618     harmHH(i, m-1) += d * e_m[i];
01619 
01620   // Revise to do query for optimal workspace first
01621   // Create simple storage for the left eigenvectors, which we don't care about.
01622   const int ldvl = m;
01623   ScalarType* vl = 0;
01624   lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
01625         vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01626   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01627 
01628   // Construct magnitude of each harmonic Ritz value
01629   for( i=0; i<m; ++i )
01630     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01631 
01632   // Construct magnitude of each harmonic Ritz value
01633   this->sort(w, m, iperm);
01634 
01635   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01636   if (wi[iperm[recycledBlocks_-1]] != zero) {
01637     int countImag = 0;
01638     for ( i=0; i<recycledBlocks_; ++i ) {
01639       if (wi[iperm[i]] != zero)
01640   countImag++;
01641     }
01642     // Check to see if this count is even or odd:
01643     if (countImag % 2)
01644       xtraVec = true;
01645   }
01646 
01647   // Select recycledBlocks_ smallest eigenvectors
01648   for( i=0; i<recycledBlocks_; ++i ) {
01649     for( j=0; j<m; j++ ) {
01650       PP(j,i) = vr(j,iperm[i]);
01651     }
01652   }
01653   
01654   if (xtraVec) { // we need to store one more vector
01655     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
01656       for( j=0; j<m; ++j ) {   // so get the "imag" component
01657   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
01658       }
01659     }
01660     else { //  I picked the "imag" component
01661       for( j=0; j<m; ++j ) {   // so get the "real" component
01662   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
01663       }
01664     }
01665   }
01666 
01667   // Return whether we needed to store an additional vector
01668   if (xtraVec) {
01669     return recycledBlocks_+1;
01670   }
01671   return recycledBlocks_;
01672 }
01673 
01674 //  Compute the harmonic eigenpairs of the projected, dense system.
01675 template<class ScalarType, class MV, class OP>
01676 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs2(int keff, int m, 
01677                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01678                  const Teuchos::RCP<const MV>& VV,
01679                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01680   int i, j;
01681   int m2 = HH.numCols();
01682   bool xtraVec = false;
01683   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01684   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01685   std::vector<int> index;
01686   
01687   // Real and imaginary eigenvalue components
01688   std::vector<ScalarType> wr(m2), wi(m2);
01689 
01690   // Magnitude of harmonic Ritz values
01691   std::vector<MagnitudeType> w(m2);
01692 
01693   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01694   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,false);
01695 
01696   // Sorted order of harmonic Ritz values
01697   std::vector<int> iperm(m2);
01698 
01699   // Form matrices for generalized eigenproblem
01700   
01701   // B = H2' * H2; Don't zero out matrix when constructing
01702   Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,false);
01703   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
01704   
01705   // A_tmp = | C'*U        0 |
01706   //         | V_{m+1}'*U  I |
01707   Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keff+m+1, keff+m );
01708 
01709   // A_tmp(1:keff,1:keff) = C' * U;
01710   index.resize(keff);
01711   for (int i=0; i<keff; ++i) { index[i] = i; }
01712   Teuchos::RCP<MV> Ctmp  = MVT::CloneView( *C_, index );
01713   Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_, index );
01714   Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keff, keff );
01715   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
01716 
01717   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
01718   Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keff, keff );
01719   index.resize(m+1);
01720   for (i=0; i < m+1; i++) { index[i] = i; }
01721   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
01722   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
01723 
01724   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
01725   for( i=keff; i<keff+m; i++ ) {
01726     A_tmp(i,i) = one;
01727   }
01728 
01729   // A = H2' * A_tmp;
01730   Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
01731   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
01732 
01733   // Compute k smallest harmonic Ritz pairs
01734   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
01735   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
01736   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
01737   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
01738   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
01739   char balanc='P', jobvl='N', jobvr='V', sense='N'; 
01740   int ld = A.numRows();
01741   int lwork = 6*ld;
01742   int ldvl = ld, ldvr = ld;
01743   int info = 0,ilo = 0,ihi = 0;
01744   ScalarType abnrm = zero, bbnrm = zero;
01745   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
01746   std::vector<ScalarType> beta(ld);
01747   std::vector<ScalarType> work(lwork);
01748   std::vector<MagnitudeType> lscale(ld), rscale(ld);  
01749   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
01750   std::vector<int> iwork(ld+6);
01751   int *bwork = 0; // If sense == 'N', bwork is never referenced
01752   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0], 
01753                &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0], 
01754                &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
01755   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
01756   
01757   // Construct magnitude of each harmonic Ritz value
01758   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
01759   for( i=0; i<ld; i++ ) {
01760     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
01761   }
01762 
01763   // Construct magnitude of each harmonic Ritz value
01764   this->sort(w,ld,iperm);
01765 
01766   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
01767   if (wi[iperm[ld-recycledBlocks_]] != zero) {
01768     int countImag = 0;
01769     for ( i=ld-recycledBlocks_; i<ld; i++ ) {
01770       if (wi[iperm[i]] != zero)
01771   countImag++;
01772     }
01773     // Check to see if this count is even or odd:
01774     if (countImag % 2)
01775       xtraVec = true;
01776   }
01777   
01778   // Select recycledBlocks_ smallest eigenvectors
01779   for( i=0; i<recycledBlocks_; i++ ) {
01780     for( j=0; j<ld; j++ ) {
01781       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
01782     }
01783   }
01784   
01785   if (xtraVec) { // we need to store one more vector
01786     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
01787       for( j=0; j<ld; j++ ) {   // so get the "imag" component
01788   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
01789       }
01790     }
01791     else { // I picked the "imag" component
01792       for( j=0; j<ld; j++ ) {   // so get the "real" component
01793   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
01794       }
01795     }
01796   }
01797   
01798   // Return whether we needed to store an additional vector
01799   if (xtraVec) {
01800     return recycledBlocks_+1;
01801   }
01802   return recycledBlocks_;
01803 }
01804 
01805 
01806 // This method sorts list of n floating-point numbers and return permutation vector
01807 template<class ScalarType, class MV, class OP>
01808 void GCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
01809   int l, r, j, i, flag;
01810   int    RR2;
01811   double dRR, dK;
01812   
01813   // Initialize the permutation vector.
01814   for(j=0;j<n;j++)
01815     iperm[j] = j;
01816   
01817   if (n <= 1) return;
01818   
01819   l    = n / 2 + 1;
01820   r    = n - 1;
01821   l    = l - 1;
01822   dRR  = dlist[l - 1];
01823   dK   = dlist[l - 1];
01824   
01825   RR2 = iperm[l - 1];
01826   while (r != 0) {
01827     j = l;
01828     flag = 1;
01829     
01830     while (flag == 1) {
01831       i = j;
01832       j = j + j;
01833       
01834       if (j > r + 1)
01835   flag = 0;
01836       else {
01837   if (j < r + 1)
01838     if (dlist[j] > dlist[j - 1]) j = j + 1;
01839   
01840   if (dlist[j - 1] > dK) {
01841     dlist[i - 1] = dlist[j - 1];
01842     iperm[i - 1] = iperm[j - 1];
01843   }
01844   else {
01845     flag = 0;
01846   }
01847       }
01848     }
01849     dlist[i - 1] = dRR;
01850     iperm[i - 1] = RR2;
01851     
01852     if (l == 1) {
01853       dRR  = dlist [r];
01854       RR2 = iperm[r];
01855       dK = dlist[r];
01856       dlist[r] = dlist[0];
01857       iperm[r] = iperm[0];
01858       r = r - 1;
01859     }
01860     else {
01861       l   = l - 1;
01862       dRR  = dlist[l - 1];
01863       RR2  = iperm[l - 1];
01864       dK   = dlist[l - 1];
01865     }
01866   }
01867   dlist[0] = dRR;
01868   iperm[0] = RR2; 
01869 }
01870 
01871 
01872 //  This method requires the solver manager to return a string that describes itself.
01873 template<class ScalarType, class MV, class OP>
01874 std::string GCRODRSolMgr<ScalarType,MV,OP>::description() const {
01875   std::ostringstream oss;
01876   oss << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01877   oss << "{";
01878   oss << "Ortho Type='"<<orthoType_;
01879   oss << ", Num Blocks=" <<numBlocks_;
01880   oss << ", Num Recycle Blocks=" << recycledBlocks_;
01881   oss << ", Max Restarts=" << maxRestarts_;
01882   oss << "}";
01883   return oss.str();
01884 }
01885   
01886 } // end Belos namespace
01887 
01888 #endif /* BELOS_GCRODR_SOLMGR_HPP */

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7