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