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 #include "BelosOrthoManagerFactory.hpp"
00052 
00053 #include "BelosLinearProblem.hpp"
00054 #include "BelosSolverManager.hpp"
00055 
00056 #include "BelosGCRODRIter.hpp"
00057 #include "BelosBlockFGmresIter.hpp"
00058 #include "BelosStatusTestMaxIters.hpp"
00059 #include "BelosStatusTestGenResNorm.hpp"
00060 #include "BelosStatusTestCombo.hpp"
00061 #include "BelosStatusTestOutputFactory.hpp"
00062 #include "BelosOutputManager.hpp"
00063 #include "Teuchos_BLAS.hpp"
00064 #include "Teuchos_LAPACK.hpp"
00065 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00066 #include "Teuchos_TimeMonitor.hpp"
00067 #endif // BELOS_TEUCHOS_TIME_MONITOR
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     typedef OrthoManagerFactory<ScalarType, MV, OP> ortho_factory_type;
00146 
00147     
00148   public:
00150 
00151    
00157      GCRODRSolMgr();
00158  
00212     GCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00213       const Teuchos::RCP<Teuchos::ParameterList> &pl);
00214     
00216     virtual ~GCRODRSolMgr() {};
00218     
00220 
00221     
00224     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00225       return *problem_;
00226     }
00227 
00230     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00231 
00234     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { 
00235       return params_; 
00236     }
00237  
00243     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00244       return Teuchos::tuple(timerSolve_);
00245     }
00246   
00248     int getNumIters() const {
00249       return numIters_;
00250     }
00251  
00254     bool isLOADetected() const { return false; }
00255  
00257     
00259 
00260    
00262     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { 
00263       problem_ = problem; 
00264     }
00265    
00267     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00268     
00270    
00272 
00273 
00277     void reset( const ResetType type ) { 
00278       if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
00279       else if (type & Belos::RecycleSubspace) keff = 0;
00280     }
00282  
00284 
00285     
00302     ReturnType solve();
00303     
00305     
00308     
00310     std::string description() const;
00311     
00313     
00314   private:
00315 
00316     // Called by all constructors; Contains init instructions common to all constructors
00317     void init();
00318 
00319     // Initialize solver state storage
00320     void initializeStateStorage();
00321 
00322     // Compute updated recycle space given existing recycle space and newly generated Krylov space
00323     void buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
00324 
00325     //  Computes harmonic eigenpairs of projected matrix created during the priming solve.
00326     //  HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
00327     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00328     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00329     int getHarmonicVecs1(int m, 
00330        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00331        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00332 
00333     //  Computes harmonic eigenpairs of projected matrix created during one cycle.
00334     //  HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
00335     //  VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
00336     //  PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
00337     //  The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
00338     int getHarmonicVecs2(int keff, int m, 
00339        const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
00340        const Teuchos::RCP<const MV>& VV,
00341        Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
00342 
00343     // Sort list of n floating-point numbers and return permutation vector
00344     void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00345 
00346     // Method to convert string to enumerated type for residual.
00347     Belos::ScaleType convertStringToScaleType( std::string& scaleType ) {
00348       if (scaleType == "Norm of Initial Residual") {
00349   return Belos::NormOfInitRes;
00350       } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00351   return Belos::NormOfPrecInitRes;
00352       } else if (scaleType == "Norm of RHS") {
00353   return Belos::NormOfRHS;
00354       } else if (scaleType == "None") {
00355   return Belos::None;
00356       } else 
00357   TEST_FOR_EXCEPTION( true ,std::logic_error,"Belos::GCRODRSolMgr(): Invalid residual scaling type.");
00358     }
00359 
00360     // Lapack interface
00361     Teuchos::LAPACK<int,ScalarType> lapack;
00362 
00363     // Linear problem.
00364     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00365     
00366     // Output manager.
00367     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00368     Teuchos::RCP<std::ostream> outputStream_;
00369 
00370     // Status test.
00371     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00372     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00373     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
00374     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
00375     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00376 
00379     ortho_factory_type orthoFactory_;
00380 
00384     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00385     
00386     // Current parameter list.
00387     Teuchos::RCP<Teuchos::ParameterList> params_;
00388 
00389     // Default solver values.
00390     static const MagnitudeType convTol_default_;
00391     static const MagnitudeType orthoKappa_default_;
00392     static const int maxRestarts_default_;
00393     static const int maxIters_default_;
00394     static const int numBlocks_default_;
00395     static const int blockSize_default_;
00396     static const int recycledBlocks_default_;
00397     static const int verbosity_default_;
00398     static const int outputStyle_default_;
00399     static const int outputFreq_default_;
00400     static const std::string impResScale_default_; 
00401     static const std::string expResScale_default_; 
00402     static const std::string label_default_;
00403     static const std::string orthoType_default_;
00404     static const Teuchos::RCP<std::ostream> outputStream_default_;
00405 
00406     // Current solver values.
00407     MagnitudeType convTol_, orthoKappa_;
00408     int maxRestarts_, maxIters_, numIters_;
00409     int verbosity_, outputStyle_, outputFreq_;
00410     std::string orthoType_; 
00411     std::string impResScale_, expResScale_;
00412 
00414     // Solver State Storage
00416     // 
00417     // The number of blocks and recycle blocks (m and k, respectively)
00418     int numBlocks_, recycledBlocks_;
00419     // Current size of recycled subspace 
00420     int keff;
00421     //
00422     // Residual vector
00423     Teuchos::RCP<MV> r_;
00424     //
00425     // Search space
00426     Teuchos::RCP<MV> V_;
00427     //
00428     // Recycled subspace and its image
00429     Teuchos::RCP<MV> U_, C_;
00430     //
00431     // Updated recycle space and its image
00432     Teuchos::RCP<MV> U1_, C1_;
00433     //
00434     // Storage used in constructing harmonic Ritz values/vectors
00435     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
00436     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00437     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00438     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
00439     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
00440     std::vector<ScalarType> tau_;
00441     std::vector<ScalarType> work_;
00442     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
00443     std::vector<int> ipiv_;
00445 
00446     // Timers.
00447     std::string label_;
00448     Teuchos::RCP<Teuchos::Time> timerSolve_;
00449 
00450     // Internal state variables.
00451     bool isSet_;
00452 
00453     // Have we generated or regenerated a recycle space yet this solve?
00454     bool builtRecycleSpace_;
00455   };
00456 
00457 
00458 // Default solver values.
00459 template<class ScalarType, class MV, class OP>
00460 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::convTol_default_ = 1e-8;
00461 
00462 template<class ScalarType, class MV, class OP>
00463 const typename GCRODRSolMgr<ScalarType,MV,OP>::MagnitudeType GCRODRSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = 0.0;
00464 
00465 template<class ScalarType, class MV, class OP>
00466 const int GCRODRSolMgr<ScalarType,MV,OP>::maxRestarts_default_ = 100;
00467 
00468 template<class ScalarType, class MV, class OP>
00469 const int GCRODRSolMgr<ScalarType,MV,OP>::maxIters_default_ = 5000;
00470 
00471 template<class ScalarType, class MV, class OP>
00472 const int GCRODRSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 50;
00473 
00474 template<class ScalarType, class MV, class OP>
00475 const int GCRODRSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
00476 
00477 template<class ScalarType, class MV, class OP>
00478 const int GCRODRSolMgr<ScalarType,MV,OP>::recycledBlocks_default_ = 5;
00479 
00480 template<class ScalarType, class MV, class OP>
00481 const int GCRODRSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00482 
00483 template<class ScalarType, class MV, class OP>
00484 const int GCRODRSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00485 
00486 template<class ScalarType, class MV, class OP>
00487 const int GCRODRSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00488 
00489 template<class ScalarType, class MV, class OP>
00490 const std::string GCRODRSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
00491 
00492 template<class ScalarType, class MV, class OP>
00493 const std::string GCRODRSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
00494 
00495 template<class ScalarType, class MV, class OP>
00496 const std::string GCRODRSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00497 
00498 template<class ScalarType, class MV, class OP>
00499 const std::string GCRODRSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00500 
00501 template<class ScalarType, class MV, class OP>
00502 const Teuchos::RCP<std::ostream> GCRODRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00503 
00504 
00505 // Empty Constructor
00506 template<class ScalarType, class MV, class OP>
00507 GCRODRSolMgr<ScalarType,MV,OP>::GCRODRSolMgr() {
00508   init();
00509 }
00510 
00511 
00512 // Basic Constructor
00513 template<class ScalarType, class MV, class OP>
00514 GCRODRSolMgr<ScalarType,MV,OP>::
00515 GCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00516        const Teuchos::RCP<Teuchos::ParameterList> &pl ) 
00517 {
00518   // Initialize local pointers to null, and initialize local variables
00519   // to default values.
00520   init();
00521 
00522   TEST_FOR_EXCEPTION(problem == Teuchos::null, std::invalid_argument, 
00523          "Belos::GCRODRSolMgr constructor: The solver manager's "
00524          "constructor needs the linear problem argument 'problem' "
00525          "to be non-null.");
00526   problem_ = problem;
00527 
00528   // Set the parameters using the list that was passed in.  If null,
00529   // we defer initialization until a non-null list is set (by the
00530   // client calling setParameters(), or by calling solve() -- in
00531   // either case, a null parameter list indicates that default
00532   // parameters should be used).
00533   if (! is_null (pl))
00534     setParameters (pl);
00535 }
00536 
00537 // Common instructions executed in all constructors
00538 template<class ScalarType, class MV, class OP>
00539 void GCRODRSolMgr<ScalarType,MV,OP>::init() {
00540   outputStream_ = outputStream_default_;
00541   convTol_ = convTol_default_;
00542   orthoKappa_ = orthoKappa_default_;
00543   maxRestarts_ = maxRestarts_default_;
00544   maxIters_ = maxIters_default_;
00545   numBlocks_ = numBlocks_default_;
00546   recycledBlocks_ = recycledBlocks_default_;
00547   verbosity_ = verbosity_default_;
00548   outputStyle_ = outputStyle_default_;
00549   outputFreq_ = outputFreq_default_;
00550   orthoType_ = orthoType_default_;
00551   impResScale_ = impResScale_default_;
00552   expResScale_ = expResScale_default_;
00553   label_ = label_default_;
00554   isSet_ = false;
00555   builtRecycleSpace_ = false;
00556   keff = 0;
00557   r_ = Teuchos::null;
00558   V_ = Teuchos::null;
00559   U_ = Teuchos::null; 
00560   C_ = Teuchos::null;
00561   U1_ = Teuchos::null;
00562   C1_ = Teuchos::null;
00563   PP_ = Teuchos::null;
00564   HP_ = Teuchos::null;
00565   H2_ = Teuchos::null;
00566   R_ = Teuchos::null;
00567   H_ = Teuchos::null;
00568   B_ = Teuchos::null;
00569 }
00570 
00571 template<class ScalarType, class MV, class OP>
00572 void 
00573 GCRODRSolMgr<ScalarType,MV,OP>::
00574 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
00575 {
00576   using Teuchos::isParameterType;
00577   using Teuchos::getParameter;
00578   using Teuchos::null;
00579   using Teuchos::ParameterList;
00580   using Teuchos::parameterList;
00581   using Teuchos::RCP;
00582   using Teuchos::rcp;
00583   using Teuchos::rcp_dynamic_cast;
00584   using Teuchos::rcpFromRef;
00585   using Teuchos::Exceptions::InvalidParameter;
00586   using Teuchos::Exceptions::InvalidParameterName;
00587   using Teuchos::Exceptions::InvalidParameterType;
00588 
00589   // The default parameter list contains all parameters that
00590   // GCRODRSolMgr understands, and none that it doesn't understand.
00591   RCP<const ParameterList> defaultParams = getValidParameters();
00592 
00593   // Create the internal parameter list if one doesn't already exist.
00594   //
00595   // (mfh 28 Feb 2011, 10 Mar 2011) At the time this code was written,
00596   // ParameterList did not have validators or validateParameters().
00597   // This is why the code below carefully validates the parameters one
00598   // by one and fills in defaults.  This code could be made a lot
00599   // shorter by using validators.  To do so, we would have to define
00600   // appropriate validators for all the parameters.  (This would more
00601   // or less just move all that validation code out of this routine
00602   // into to getValidParameters().)
00603   //
00604   // For an analogous reason, GCRODRSolMgr defines default parameter
00605   // values as class data, as well as in the default ParameterList.
00606   // This redundancy could be removed by defining the default
00607   // parameter values only in the default ParameterList (which
00608   // documents each parameter as well -- handy!).
00609   if (params_.is_null()) {
00610     params_ = parameterList (*defaultParams);
00611   } else {
00612     // A common case for setParameters() is for it to be called at the
00613     // beginning of the solve() routine.  This follows the Belos
00614     // pattern of delaying initialization until the last possible
00615     // moment (when the user asks Belos to perform the solve).  In
00616     // this common case, we save ourselves a deep copy of the input
00617     // parameter list.
00618     if (params_ != params)
00619       {
00620   // Make a deep copy of the input parameter list.  This allows
00621   // the caller to modify or change params later, without
00622   // affecting the behavior of this solver.  This solver will
00623   // then only change its internal parameters if setParameters()
00624   // is called again.
00625   params_ = parameterList (*params);
00626       }
00627 
00628     // Fill in any missing parameters and their default values.  Also,
00629     // throw an exception if the parameter list has any misspelled or
00630     // "extra" parameters.  If you don't like this behavior, you'll
00631     // want to replace the line of code below with your desired
00632     // validation scheme.  Note that Teuchos currently only implements
00633     // two options: 
00634     //
00635     // 1. validateParameters() requires that params_ has all the
00636     //    parameters that the default list has, and none that it
00637     //    doesn't have.
00638     //
00639     // 2. validateParametersAndSetDefaults() fills in missing
00640     //    parameters in params_ using the default list, but requires
00641     //    that any parameter provided in params_ is also in the
00642     //    default list.
00643     //
00644     // Here is an easy way to ignore any "extra" or misspelled
00645     // parameters: Make a deep copy of the default list, fill in any
00646     // "missing" parameters from the _input_ list, and then validate
00647     // the input list using the deep copy of the default list.  We
00648     // show this in code:
00649     //
00650     // RCP<ParameterList> defaultCopy = parameterList (*getValidParameters ());
00651     // defaultCopy->validateParametersAndSetDefaults (params);
00652     // params->validateParametersAndSetDefaults (defaultCopy);
00653     //
00654     // This method is not entirely robust, because the input list may
00655     // have incorrect validators set for existing parameters in the
00656     // default list.  This would then cause "validation" of the
00657     // default list to throw an exception.  As a result, we've chosen
00658     // for now to be intolerant of misspellings and "extra" parameters
00659     // in the input list.
00660     params_->validateParametersAndSetDefaults (*defaultParams);
00661   }
00662 
00663   // Check for maximum number of restarts.
00664   if (params->isParameter ("Maximum Restarts")) {
00665     maxRestarts_ = params->get("Maximum Restarts", maxRestarts_default_);
00666 
00667     // Update parameter in our list.
00668     params_->set ("Maximum Restarts", maxRestarts_);
00669   }
00670 
00671   // Check for maximum number of iterations
00672   if (params->isParameter ("Maximum Iterations")) {
00673     maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
00674 
00675     // Update parameter in our list and in status test.
00676     params_->set ("Maximum Iterations", maxIters_);
00677     if (! maxIterTest_.is_null())
00678       maxIterTest_->setMaxIters (maxIters_);
00679   }
00680 
00681   // Check for the maximum number of blocks.
00682   if (params->isParameter ("Num Blocks")) {
00683     numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
00684     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00685            "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
00686            "be strictly positive, but you specified a value of "
00687            << numBlocks_ << ".");
00688     // Update parameter in our list.
00689     params_->set ("Num Blocks", numBlocks_);
00690   }
00691 
00692   // Check for the maximum number of blocks.
00693   if (params->isParameter ("Num Recycled Blocks")) {
00694     recycledBlocks_ = params->get ("Num Recycled Blocks", 
00695            recycledBlocks_default_);
00696     TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
00697            "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
00698            "parameter must be strictly positive, but you specified "
00699            "a value of " << recycledBlocks_ << ".");
00700     TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
00701            "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
00702            "parameter must be less than the \"Num Blocks\" "
00703            "parameter, but you specified \"Num Recycled Blocks\" "
00704            "= " << recycledBlocks_ << " and \"Num Blocks\" = " 
00705            << numBlocks_ << ".");
00706     // Update parameter in our list.
00707     params_->set("Num Recycled Blocks", recycledBlocks_);
00708   }
00709 
00710   // Check to see if the timer label changed.  If it did, update it in
00711   // the parameter list, and create a new timer with that label (if
00712   // Belos was compiled with timers enabled).
00713   if (params->isParameter ("Timer Label")) {
00714     std::string tempLabel = params->get ("Timer Label", label_default_);
00715 
00716     // Update parameter in our list and solver timer
00717     if (tempLabel != label_) {
00718       label_ = tempLabel;
00719       params_->set ("Timer Label", label_);
00720       std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
00721 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00722       timerSolve_ = Teuchos::TimeMonitor::getNewTimer (solveLabel);
00723 #endif
00724     }
00725   }
00726 
00727   // Check for a change in verbosity level
00728   if (params->isParameter ("Verbosity")) {
00729     if (isParameterType<int> (*params, "Verbosity")) {
00730       verbosity_ = params->get ("Verbosity", verbosity_default_);
00731     } else {
00732       verbosity_ = (int) getParameter<Belos::MsgType> (*params, "Verbosity");
00733     }
00734     // Update parameter in our list.
00735     params_->set ("Verbosity", verbosity_);
00736     // If the output manager (printer_) is null, then we will
00737     // instantiate it later with the correct verbosity.
00738     if (! printer_.is_null())
00739       printer_->setVerbosity (verbosity_);
00740   }
00741 
00742   // Check for a change in output style
00743   if (params->isParameter ("Output Style")) {
00744     if (isParameterType<int> (*params, "Output Style")) {
00745       outputStyle_ = params->get ("Output Style", outputStyle_default_);
00746     } else {
00747       outputStyle_ = (int) getParameter<OutputType> (*params, "Output Style");
00748     }
00749 
00750     // Update parameter in our list.
00751     params_->set ("Output Style", outputStyle_);
00752     // We will (re)instantiate the output status test afresh below.
00753     outputTest_ = null;
00754   }
00755 
00756   // Get the output stream for the output manager.
00757   //
00758   // While storing the output stream in the parameter list (either as
00759   // an RCP or as a nonconst reference) is convenient and safe for
00760   // programming, it makes it impossible to serialize the parameter
00761   // list, read it back in from the serialized representation, and get
00762   // the same output stream as before.  This is because output streams
00763   // may be arbitrary constructed objects.
00764   //
00765   // In case the user tried reading in the parameter list from a
00766   // serialized representation and the output stream can't be read
00767   // back in, we set the output stream to point to std::cout.  This
00768   // ensures reasonable behavior.
00769   if (params->isParameter ("Output Stream")) {
00770     try {
00771       outputStream_ = getParameter<RCP<std::ostream> > (*params, "Output Stream");
00772     } catch (InvalidParameter&) {
00773       outputStream_ = rcpFromRef (std::cout);
00774     }
00775     // We assume that a null output stream indicates that the user
00776     // doesn't want to print anything, so we replace it with a "black
00777     // hole" stream that prints nothing sent to it.  (We can't use a
00778     // null output stream, since the output manager always sends
00779     // things it wants to print to the output stream.)
00780     if (outputStream_.is_null())
00781       outputStream_ = rcp (new Teuchos::oblackholestream);
00782 
00783     // Update parameter in our list.
00784     params_->set ("Output Stream", outputStream_);
00785     // If the output manager (printer_) is null, then we will
00786     // instantiate it later with the correct output stream.
00787     if (! printer_.is_null())
00788       printer_->setOStream (outputStream_);
00789   }
00790 
00791   // frequency level
00792   if (verbosity_ & Belos::StatusTestDetails) {
00793     if (params->isParameter ("Output Frequency")) {
00794       outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
00795     }
00796 
00797     // Update parameter in out list and output status test.
00798     params_->set("Output Frequency", outputFreq_);
00799     if (! outputTest_.is_null())
00800       outputTest_->setOutputFrequency (outputFreq_);
00801   }
00802 
00803   // Create output manager if we need to, using the verbosity level
00804   // and output stream that we fetched above.  We do this here because
00805   // instantiating an OrthoManager using OrthoManagerFactory requires
00806   // a valid OutputManager.
00807   if (printer_.is_null()) {
00808     printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
00809   }  
00810 
00811   // Get the orthogonalization manager name ("Orthogonalization").
00812   //
00813   // Getting default values for the orthogonalization manager
00814   // parameters ("Orthogonalization Parameters") requires knowing the
00815   // orthogonalization manager name.  Save it for later, and also
00816   // record whether it's different than before.
00817   bool changedOrthoType = false;
00818   if (params->isParameter ("Orthogonalization")) 
00819     {
00820       const std::string& tempOrthoType = 
00821   params->get ("Orthogonalization", orthoType_default_);
00822       // Ensure that the specified orthogonalization type is valid.
00823       if (! orthoFactory_.isValidName (tempOrthoType))
00824   {
00825     std::ostringstream os;
00826     os << "Belos::GCRODRSolMgr: Invalid orthogonalization name \"" 
00827        << tempOrthoType << "\".  The following are valid options "
00828        << "for the \"Orthogonalization\" name parameter: ";
00829     orthoFactory_.printValidNames (os);
00830     throw std::invalid_argument (os.str());
00831   }
00832       if (tempOrthoType != orthoType_)
00833   {
00834     changedOrthoType = true;
00835     orthoType_ = tempOrthoType;
00836     // Update parameter in our list.
00837     params_->set ("Orthogonalization", orthoType_);
00838   }
00839     }
00840 
00841   // Get any parameters for the orthogonalization ("Orthogonalization
00842   // Parameters").  If not supplied, the orthogonalization manager
00843   // factory will supply default values.
00844   //
00845   // NOTE (mfh 12 Jan 2011) For the sake of backwards compatibility,
00846   // if params has an "Orthogonalization Constant" parameter and the
00847   // DGKS orthogonalization manager is to be used, the value of this
00848   // parameter will override DGKS's "depTol" parameter.
00849   //
00850   // Users may supply the orthogonalization manager parameters either
00851   // as a sublist, or as an RCP.  We test for both.  Note that setting
00852   // the sublist as an RCP<const ParameterList> rather than a
00853   // ParameterList means that you can't correctly serialize the solver
00854   // manager's parameter list.
00855   RCP<const ParameterList> orthoParams;
00856   {
00857     bool gotOrthoParams = false;
00858     try { // Could it be an RCP?
00859       orthoParams = 
00860   params->get<RCP<const ParameterList> >("Orthogonalization Parameters");
00861       gotOrthoParams = true;
00862     } catch (InvalidParameter&) {
00863       // We didn't get orthoParams; gotOrthoParams stays false.
00864     }
00865     if (! gotOrthoParams) {
00866       try { // Could it be a sublist?
00867   const ParameterList& _orthoParams = 
00868     params_->sublist("Orthogonalization Parameters");
00869   // A deep copy is the only safe way to ensure that
00870   // orthoParams doesn't "go away," since params doesn't
00871   // belong to the solution manager and may fall out of
00872   // scope.
00873   orthoParams = rcp (new ParameterList (_orthoParams));
00874   gotOrthoParams = true;
00875       } catch (InvalidParameter&) {
00876   // We didn't get orthoParams; gotOrthoParams stays false.
00877       }
00878     }
00879     // We didn't get the parameter list from params, so get a default
00880     // parameter list from the OrthoManagerFactory.
00881     if (! gotOrthoParams)
00882       orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
00883     // Update parameter in our list.
00884     params_->set ("Orthogonalization Parameters", orthoParams);
00885   }
00886 
00887   // Check if the desired orthogonalization method changed, or if the
00888   // orthogonalization manager has not yet been instantiated.  If
00889   // either is the case, instantiate a new MatOrthoManager subclass
00890   // instance corresponding to the desired orthogonalization method.
00891   // We've already fetched the orthogonalization method name
00892   // (orthoType_) and its parameters (orthoParams) above.
00893   //
00894   // NOTE (mfh 12 Jan 2011) We only instantiate a new MatOrthoManager
00895   // subclass if the orthogonalization method name is different than
00896   // before.  Thus, for some orthogonalization managers, changes to
00897   // their parameters may not get propagated, if the manager type
00898   // itself didn't change.  The one exception is the "depTol"
00899   // (a.k.a. orthoKappa or "Orthogonalization Constant") parameter of
00900   // DGKS; changes to that _do_ get propagated down to the DGKS
00901   // instance.  
00902   //
00903   // The most general way to fix this issue would be to supply each
00904   // orthogonalization manager class with a setParameters() method
00905   // that takes a parameter list input, and changes the parameters as
00906   // appropriate.  A less efficient but correct way would be simply to
00907   // reinstantiate the OrthoManager every time, whether or not the
00908   // orthogonalization method name or parameters have changed.
00909   if (ortho_.is_null() || changedOrthoType)
00910     {
00911       // Create orthogonalization manager.  This requires that the
00912       // OutputManager (printer_) already be initialized.
00913       ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_, 
00914               label_, orthoParams);
00915     }
00916 
00917   // The DGKS orthogonalization accepts a "Orthogonalization Constant"
00918   // parameter (also called kappa in the code, but not in the
00919   // parameter list).  If its value is provided in the given parameter
00920   // list, and its value is positive, use it.  Ignore negative values.
00921   //
00922   // NOTE (mfh 12 Jan 2011) This overrides the "depTol" parameter that
00923   // may have been specified in "Orthogonalization Parameters".  We
00924   // retain this behavior for backwards compatibility.
00925   bool gotValidOrthoKappa = false;
00926   if (params->isParameter ("Orthogonalization Constant"))
00927     {
00928       const MagnitudeType orthoKappa = 
00929   params->get ("Orthogonalization Constant", orthoKappa_default_);
00930       if (orthoKappa > 0)
00931   {
00932     orthoKappa_ = orthoKappa;
00933     gotValidOrthoKappa = true;
00934     // Update parameter in our list.
00935     params_->set("Orthogonalization Constant", orthoKappa_);
00936     // Only DGKS currently accepts this parameter.
00937     if (orthoType_ == "DGKS" && ! ortho_.is_null())
00938       {
00939         typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
00940         // This cast should always succeed; it's a bug
00941         // otherwise.  (If the cast fails, then orthoType_
00942         // doesn't correspond to the OrthoManager subclass
00943         // instance that we think we have, so we initialized the
00944         // wrong subclass somehow.)
00945         rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
00946       }
00947   }
00948     }
00949   
00950   // Convergence
00951   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00952   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00953 
00954   // Check for convergence tolerance
00955   if (params->isParameter("Convergence Tolerance")) {
00956     convTol_ = params->get ("Convergence Tolerance", convTol_default_);
00957 
00958     // Update parameter in our list and residual tests.
00959     params_->set ("Convergence Tolerance", convTol_);
00960     if (! impConvTest_.is_null())
00961       impConvTest_->setTolerance (convTol_);
00962     if (! expConvTest_.is_null())
00963       expConvTest_->setTolerance (convTol_);
00964   }
00965  
00966   // Check for a change in scaling, if so we need to build new residual tests.
00967   if (params->isParameter ("Implicit Residual Scaling")) {
00968     std::string tempImpResScale = 
00969       getParameter<std::string> (*params, "Implicit Residual Scaling");
00970 
00971     // Only update the scaling if it's different.
00972     if (impResScale_ != tempImpResScale) {
00973       ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
00974       impResScale_ = tempImpResScale;
00975 
00976       // Update parameter in our list and residual tests
00977       params_->set("Implicit Residual Scaling", impResScale_);
00978       // NOTE (mfh 28 Feb 2011) StatusTestImpResNorm only lets you
00979       // call defineScaleForm() once.  The code below attempts to call
00980       // defineScaleForm(); if the scale form has already been
00981       // defined, it constructs a new StatusTestImpResNorm instance.
00982       // StatusTestImpResNorm really should not expose the
00983       // defineScaleForm() method, since it's serving an
00984       // initialization purpose; all initialization should happen in
00985       // the constructor whenever possible.  In that case, the code
00986       // below could be simplified into a single (re)instantiation.
00987       if (! impConvTest_.is_null()) {
00988         try { 
00989           impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
00990         }
00991         catch (StatusTestError&) {
00992           // Delete the convergence test so it gets constructed again.
00993     impConvTest_ = null;
00994           convTest_ = null;
00995         }
00996       }
00997     }      
00998   }
00999   
01000   if (params->isParameter("Explicit Residual Scaling")) {
01001     std::string tempExpResScale = 
01002       getParameter<std::string> (*params, "Explicit Residual Scaling");
01003 
01004     // Only update the scaling if it's different.
01005     if (expResScale_ != tempExpResScale) {
01006       ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
01007       expResScale_ = tempExpResScale;
01008 
01009       // Update parameter in our list and residual tests
01010       params_->set("Explicit Residual Scaling", expResScale_);
01011       // NOTE (mfh 28 Feb 2011) See note above on the (re)construction
01012       // of StatusTestImpResNorm.
01013       if (! expConvTest_.is_null()) {
01014         try { 
01015           expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
01016         }
01017         catch (StatusTestError&) {
01018           // Delete the convergence test so it gets constructed again.
01019     expConvTest_ = null;
01020           convTest_ = null;
01021         }
01022       }
01023     }      
01024   }
01025   //
01026   // Create iteration stopping criteria ("status tests") if we need
01027   // to, by combining three different stopping criteria.
01028   //
01029   // First, construct maximum-number-of-iterations stopping criterion.
01030   if (maxIterTest_.is_null())
01031     maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
01032 
01033   // Implicit residual test, using the native residual to determine if
01034   // convergence was achieved.
01035   if (impConvTest_.is_null()) {
01036     impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01037     impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), 
01038            Belos::TwoNorm);
01039   }
01040 
01041   // Explicit residual test once the native residual is below the tolerance
01042   if (expConvTest_.is_null()) {
01043     expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
01044     expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
01045     expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), 
01046            Belos::TwoNorm);
01047   }
01048   // Convergence test first tests the implicit residual, then the
01049   // explicit residual if the implicit residual test passes.
01050   if (convTest_.is_null()) {
01051     convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ, 
01052               impConvTest_, 
01053               expConvTest_));
01054   }
01055   // Construct the complete iteration stopping criterion:
01056   //
01057   // "Stop iterating if the maximum number of iterations has been
01058   // reached, or if the convergence test passes."
01059   sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, 
01060                maxIterTest_, 
01061                convTest_));
01062   // Create the status test output class.
01063   // This class manages and formats the output from the status test.
01064   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
01065   outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_, 
01066            Passed+Failed+Undefined);
01067 
01068   // Set the solver string for the output test
01069   std::string solverDesc = " GCRODR ";
01070   outputTest_->setSolverDesc( solverDesc );
01071 
01072   // Create the timer if we need to.
01073   if (timerSolve_.is_null()) {
01074     std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
01075 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01076     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
01077 #endif
01078   }
01079 
01080   // Inform the solver manager that the current parameters were set.
01081   isSet_ = true;
01082 }
01083 
01084     
01085 template<class ScalarType, class MV, class OP>
01086 Teuchos::RCP<const Teuchos::ParameterList> GCRODRSolMgr<ScalarType,MV,OP>::getValidParameters() const {
01087 
01088   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
01089   if (is_null(validPL)) {
01090     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
01091     // Set all the valid parameters and their default values.
01092     pl->set("Convergence Tolerance", convTol_default_,
01093       "The relative residual tolerance that needs to be achieved by the\n"
01094       "iterative solver in order for the linear system to be declared converged.");
01095     pl->set("Maximum Restarts", maxRestarts_default_,
01096       "The maximum number of cycles allowed for each\n"
01097       "set of RHS solved.");
01098     pl->set("Maximum Iterations", maxIters_default_,
01099       "The maximum number of iterations allowed for each\n"
01100       "set of RHS solved.");
01101     // mfh 25 Oct 2010: "Block Size" must be 1 because GCRODR is
01102     // currently not a block method: i.e., it does not work on
01103     // multiple right-hand sides at once.
01104     pl->set("Block Size", blockSize_default_,
01105       "Block Size Parameter -- currently must be 1 for GCRODR");
01106     pl->set("Num Blocks", numBlocks_default_,
01107       "The maximum number of vectors allowed in the Krylov subspace\n"
01108       "for each set of RHS solved.");
01109     pl->set("Num Recycled Blocks", recycledBlocks_default_,
01110       "The maximum number of vectors in the recycled subspace." );
01111     pl->set("Verbosity", verbosity_default_,
01112       "What type(s) of solver information should be outputted\n"
01113       "to the output stream.");
01114     pl->set("Output Style", outputStyle_default_,
01115       "What style is used for the solver information outputted\n"
01116       "to the output stream.");
01117     pl->set("Output Frequency", outputFreq_default_,
01118       "How often convergence information should be outputted\n"
01119       "to the output stream.");  
01120     pl->set("Output Stream", outputStream_default_,
01121       "A reference-counted pointer to the output stream where all\n"
01122       "solver output is sent.");
01123     pl->set("Implicit Residual Scaling", impResScale_default_,
01124       "The type of scaling used in the implicit residual convergence test.");
01125     pl->set("Explicit Residual Scaling", expResScale_default_,
01126       "The type of scaling used in the explicit residual convergence test.");
01127     pl->set("Timer Label", label_default_,
01128       "The string to use as a prefix for the timer labels.");
01129     //  pl->set("Restart Timers", restartTimers_);
01130     pl->set("Orthogonalization", orthoType_default_,
01131       "The type of orthogonalization to use.  Valid options: " + 
01132       orthoFactory_.validNamesString());
01133     {
01134       // We have to help out the C++ compiler's type inference a bit here.
01135       typedef Teuchos::RCP<const Teuchos::ParameterList> const_plist_ptr;
01136 #if 0
01137       const_plist_ptr orthoParams = 
01138   orthoFactory_.getDefaultParameters (orthoType_default_);
01139 #else
01140       const_plist_ptr orthoParams;
01141 #endif // 0
01142       pl->set< const_plist_ptr > ("Orthogonalization Parameters", orthoParams,
01143           "Parameters specific to the type of "
01144           "orthogonalization used.");
01145     }
01146     pl->set("Orthogonalization Constant", orthoKappa_default_,
01147       "When using DGKS orthogonalization: the \"depTol\" constant, used "
01148       "to determine whether another step of classical Gram-Schmidt is "
01149       "necessary.  Otherwise ignored.");
01150     validPL = pl;
01151   }
01152   return validPL;
01153 
01154 }
01155 
01156 // initializeStateStorage
01157 template<class ScalarType, class MV, class OP>
01158 void GCRODRSolMgr<ScalarType,MV,OP>::initializeStateStorage() {
01159 
01160     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01161 
01162     // Check if there is any multivector to clone from.
01163     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
01164     if (rhsMV == Teuchos::null) {
01165       // Nothing to do
01166       return;
01167     }
01168     else {
01169 
01170       // Initialize the state storage
01171       TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
01172                          "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
01173 
01174       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01175       if (U_ == Teuchos::null) {
01176         U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01177       }
01178       else {
01179         // Generate U_ by cloning itself ONLY if more space is needed.
01180         if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
01181           Teuchos::RCP<const MV> tmp = U_;
01182           U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01183         }
01184       }
01185 
01186       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01187       if (C_ == Teuchos::null) {
01188         C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01189       }
01190       else {
01191         // Generate C_ by cloning itself ONLY if more space is needed.
01192         if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
01193           Teuchos::RCP<const MV> tmp = C_;
01194           C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01195         }
01196       }
01197 
01198       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01199       if (V_ == Teuchos::null) {
01200         V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
01201       }
01202       else {
01203         // Generate V_ by cloning itself ONLY if more space is needed.
01204         if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
01205           Teuchos::RCP<const MV> tmp = V_;
01206           V_ = MVT::Clone( *tmp, numBlocks_+1 );
01207         }
01208       }
01209 
01210       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01211       if (U1_ == Teuchos::null) {
01212         U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01213       }
01214       else {
01215         // Generate U1_ by cloning itself ONLY if more space is needed.
01216         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01217           Teuchos::RCP<const MV> tmp = U1_;
01218           U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01219         }
01220       }
01221 
01222       // If the subspace has not been initialized before, generate it using the RHS from lp_.
01223       if (C1_ == Teuchos::null) {
01224         C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
01225       }
01226       else {
01227         // Generate C1_ by cloning itself ONLY if more space is needed.
01228         if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
01229           Teuchos::RCP<const MV> tmp = C1_;
01230           C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
01231         }
01232       }
01233 
01234       // Generate r_ only if it doesn't exist
01235       if (r_ == Teuchos::null)
01236         r_ = MVT::Clone( *rhsMV, 1 );
01237 
01238       // Size of tau_ will change during computation, so just be sure it starts with appropriate size
01239       tau_.resize(recycledBlocks_+1);
01240 
01241       // Size of work_ will change during computation, so just be sure it starts with appropriate size
01242       work_.resize(recycledBlocks_+1);
01243 
01244       // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
01245       ipiv_.resize(recycledBlocks_+1);
01246 
01247       // Generate H2_ only if it doesn't exist, otherwise resize it.
01248       if (H2_ == Teuchos::null)
01249         H2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
01250       else {
01251         if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
01252           H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
01253       }
01254       H2_->putScalar(zero);
01255 
01256       // Generate R_ only if it doesn't exist, otherwise resize it.
01257       if (R_ == Teuchos::null)
01258         R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
01259       else {
01260         if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
01261           R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
01262       }
01263       R_->putScalar(zero);
01264 
01265       // Generate PP_ only if it doesn't exist, otherwise resize it.
01266       if (PP_ == Teuchos::null)
01267         PP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
01268       else {
01269         if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
01270           PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
01271       }
01272 
01273       // Generate HP_ only if it doesn't exist, otherwise resize it.
01274       if (HP_ == Teuchos::null)
01275         HP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
01276       else {
01277         if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
01278           HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
01279       }
01280 
01281     } // end else
01282 }
01283 
01284   
01285 // solve()
01286 template<class ScalarType, class MV, class OP>
01287 ReturnType GCRODRSolMgr<ScalarType,MV,OP>::solve() {
01288   using Teuchos::RCP;
01289   using Teuchos::rcp;
01290 
01291   // Set the current parameters if they were not set before.
01292   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
01293   // then didn't set any parameters using setParameters().
01294   if (!isSet_) { setParameters( params_ ); }
01295 
01296   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01297   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01298   std::vector<int> index(numBlocks_+1);
01299   
01300   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
01301 
01302   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
01303 
01304   // Create indices for the linear systems to be solved.
01305   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
01306   std::vector<int> currIdx(1);
01307   currIdx[0] = 0;
01308 
01309   // Inform the linear problem of the current linear system to solve.
01310   problem_->setLSIndex( currIdx );
01311 
01312   // Check the number of blocks and change them is necessary.
01313   int dim = MVT::GetVecLength( *(problem_->getRHS()) );  
01314   if (numBlocks_ > dim) {
01315     numBlocks_ = dim;
01316     printer_->stream(Warnings) << 
01317       "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
01318       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
01319     params_->set("Num Blocks", numBlocks_);
01320   } 
01321 
01322   // Assume convergence is achieved, then let any failed convergence set this to false.
01323   bool isConverged = true;  
01324 
01325   // Initialize storage for all state variables
01326   initializeStateStorage();
01327 
01329   // Parameter list
01330   Teuchos::ParameterList plist;
01331   
01332   plist.set("Num Blocks",numBlocks_);
01333   plist.set("Recycled Blocks",recycledBlocks_);
01334   
01336   // GCRODR solver
01337   
01338   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
01339   gcrodr_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
01340   // Number of iterations required to generate initial recycle space (if needed)
01341   int prime_iterations = 0;
01342 
01343   // Enter solve() iterations
01344   {
01345 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01346     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01347 #endif
01348     
01349     while ( numRHS2Solve > 0 ) {
01350 
01351       // Set flag indicating recycle space has not been generated this solve
01352       builtRecycleSpace_ = false;
01353       
01354       // Reset the status test.  
01355       outputTest_->reset();
01356       
01358       // Initialize recycled subspace for GCRODR
01359       
01360       // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
01361       if (keff > 0) {
01362   TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
01363          "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
01364 
01365   printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
01366   // Compute image of U_ under the new operator
01367   index.resize(keff);
01368   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01369   RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
01370   RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01371   problem_->apply( *Utmp, *Ctmp );
01372 
01373   RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01374 
01375   // Orthogonalize this block
01376   // Get a matrix to hold the orthonormalization coefficients.
01377         Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff ); 
01378   int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,false));
01379   // Throw an error if we could not orthogonalize this block
01380   TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
01381   
01382   // U_ = U_*R^{-1} 
01383   // First, compute LU factorization of R
01384   int info = 0;
01385         ipiv_.resize(Rtmp.numRows());
01386   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01387   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01388   
01389   // Now, form inv(R)
01390   int lwork = Rtmp.numRows();
01391         work_.resize(lwork);
01392   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01393   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01394   
01395         // U_ = U1_; (via a swap)
01396   MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
01397         std::swap(U_, U1_);
01398 
01399         // Must reinitialize after swap
01400   index.resize(keff);
01401   for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01402   Ctmp  = MVT::CloneViewNonConst( *C_, index );
01403   Utmp  = MVT::CloneView( *U_, index );
01404 
01405   // Compute C_'*r_
01406   Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
01407   problem_->computeCurrPrecResVec( &*r_ );
01408   MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
01409 
01410   // Update solution ( x += U_*C_'*r_ )
01411         RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
01412         MVT::MvInit( *update, 0.0 );
01413         MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
01414         problem_->updateSolution( update, true );
01415   
01416   // Update residual norm ( r -= C_*C_'*r_ )  
01417   MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
01418 
01419         // We recycled space from previous call
01420         prime_iterations = 0;
01421 
01422       }
01423       else {
01424   
01425   // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
01426   printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
01427   
01428   Teuchos::ParameterList primeList;
01429   
01430   // Tell the block solver that the block size is one.
01431         primeList.set("Num Blocks",numBlocks_);
01432         primeList.set("Recycled Blocks",0);
01433   
01434   //  Create GCRODR iterator object to perform one cycle of GMRES.
01435   RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
01436   gcrodr_prime_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
01437   
01438   // Create the first block in the current Krylov basis (residual).
01439   problem_->computeCurrPrecResVec( &*r_ );
01440         index.resize( 1 ); index[0] = 0;
01441   RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01442   MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01443 
01444   // Set the new state and initialize the solver.
01445   GCRODRIterState<ScalarType,MV> newstate;
01446         index.resize( numBlocks_+1 );
01447         for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01448         newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01449   newstate.U = Teuchos::null;
01450   newstate.C = Teuchos::null;
01451         newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
01452   newstate.B = Teuchos::null;
01453   newstate.curDim = 0;
01454   gcrodr_prime_iter->initialize(newstate);
01455   
01456   // Perform one cycle of GMRES
01457   bool primeConverged = false;
01458   try {
01459     gcrodr_prime_iter->iterate();
01460 
01461     // Check convergence first
01462     if ( convTest_->getStatus() == Passed ) {
01463       // we have convergence
01464       primeConverged = true;
01465     }
01466   }
01467   catch (const GCRODRIterOrthoFailure &e) {
01468     // Try to recover the most recent least-squares solution
01469     gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
01470     
01471     // Check to see if the most recent least-squares solution yielded convergence.
01472     sTest_->checkStatus( &*gcrodr_prime_iter );
01473     if (convTest_->getStatus() == Passed)
01474       primeConverged = true;
01475   }
01476   catch (const std::exception &e) {
01477     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01478            << gcrodr_prime_iter->getNumIters() << std::endl 
01479            << e.what() << std::endl;
01480     throw;
01481   }
01482         // Record number of iterations in generating initial recycle spacec
01483         prime_iterations = gcrodr_prime_iter->getNumIters();
01484            
01485   // Update the linear problem.
01486   RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
01487   problem_->updateSolution( update, true );
01488 
01489   // Get the state.
01490   newstate = gcrodr_prime_iter->getState();
01491   int p = newstate.curDim;
01492 
01493   // Compute harmonic Ritz vectors 
01494   // NOTE:  The storage for the harmonic Ritz vectors (PP) is made one column larger 
01495         //        just in case we split a complex conjugate pair.
01496   // NOTE:  Generate a recycled subspace only if we have enough vectors.  If we converged
01497   //        too early, move on to the next linear system and try to generate a subspace again.
01498   if (recycledBlocks_ < p+1) {
01499     int info = 0;
01500           RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
01501           // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
01502     keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
01503           // Hereafter, only keff columns of PP are needed
01504           PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
01505           // Now get views into C, U, V
01506           index.resize(keff);
01507           for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01508           RCP<MV> Ctmp  = MVT::CloneViewNonConst( *C_, index );
01509           RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01510           RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
01511     index.resize(p);
01512           for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
01513     RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01514 
01515     // Form U (the subspace to recycle)
01516           // U = newstate.V(:,1:p) * PP;
01517           MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
01518 
01519     // Form orthonormalized C and adjust U so that C = A*U
01520 
01521     // First, compute [Q, R] = qr(H*P);
01522 
01523           // Step #1: Form HP = H*P
01524     Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
01525     Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
01526     HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
01527 
01528     // Step #1.5: Perform workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01529           int lwork = -1;
01530           tau_.resize(keff);
01531           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01532     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01533 
01534           // Step #2: Compute QR factorization of HP
01535           lwork = (int)work_[0];
01536           work_.resize(lwork);
01537           lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01538     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01539 
01540           // Step #3: Explicitly construct Q and R factors 
01541     // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01542           Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
01543           for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
01544           lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01545     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01546 
01547           // Now we have [Q,R] = qr(H*P)
01548 
01549     // Now compute C = V(:,1:p+1) * Q 
01550     index.resize( p+1 );
01551           for (int ii=0; ii < (p+1); ++ii) { index[ii] = ii; }
01552     Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
01553           MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
01554         
01555     // Finally, compute U = U*R^{-1}. 
01556           // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as 
01557           // backsolve capabilities don't exist in the Belos::MultiVec class
01558 
01559     // Step #1: First, compute LU factorization of R
01560     ipiv_.resize(Rtmp.numRows());
01561     lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01562     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01563     
01564     // Step #2: Form inv(R)
01565     lwork = Rtmp.numRows();
01566     work_.resize(lwork);
01567     lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01568     TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,  "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
01569     
01570           // Step #3: Let U = U * R^{-1}
01571     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01572 
01573     printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01574 
01575   }  // if (recycledBlocks_ < p+1)
01576 
01577         // Return to outer loop if the priming solve converged, set the next linear system.
01578   if (primeConverged) {
01579     // Inform the linear problem that we are finished with this block linear system.
01580     problem_->setCurrLS();
01581 
01582           // Update indices for the linear systems to be solved.
01583           numRHS2Solve -= 1;
01584           if ( numRHS2Solve > 0 ) {
01585             currIdx[0]++;
01586  
01587             // Set the next indices.
01588             problem_->setLSIndex( currIdx );
01589           }
01590           else {
01591             currIdx.resize( numRHS2Solve );
01592           }
01593 
01594     continue;
01595   }
01596       } // if (keff > 0) ...
01597       
01598       // Prepare for the Gmres iterations with the recycled subspace.
01599 
01600       // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01601       gcrodr_iter->setSize( keff, numBlocks_ );
01602       
01603       // Reset the number of iterations.
01604       gcrodr_iter->resetNumIters(prime_iterations);
01605 
01606       // Reset the number of calls that the status test output knows about.
01607       outputTest_->resetNumCalls();
01608 
01609       // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
01610       problem_->computeCurrPrecResVec( &*r_ );
01611       index.resize( 1 ); index[0] = 0;
01612       RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01613       MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01614 
01615       // Set the new state and initialize the solver.
01616       GCRODRIterState<ScalarType,MV> newstate;
01617       index.resize( numBlocks_+1 );
01618       for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01619       newstate.V  = MVT::CloneViewNonConst( *V_,  index );
01620       index.resize( keff );
01621       for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01622       newstate.C  = MVT::CloneViewNonConst( *C_,  index );
01623       newstate.U  = MVT::CloneViewNonConst( *U_,  index );
01624       newstate.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_,         keff, numBlocks_,    0, keff ) );
01625       newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01626       newstate.curDim = 0;
01627       gcrodr_iter->initialize(newstate);
01628 
01629       // variables needed for inner loop
01630       int numRestarts = 0;
01631       while(1) {
01632   
01633   // tell gcrodr_iter to iterate
01634   try {
01635     gcrodr_iter->iterate();
01636     
01638     //
01639     // check convergence first
01640     //
01642     if ( convTest_->getStatus() == Passed ) {
01643       // we have convergence
01644       break;  // break from while(1){gcrodr_iter->iterate()}
01645     }
01647     //
01648     // check for maximum iterations
01649     //
01651     else if ( maxIterTest_->getStatus() == Passed ) {
01652       // we don't have convergence
01653       isConverged = false;
01654       break;  // break from while(1){gcrodr_iter->iterate()}
01655     }
01657     //
01658     // check for restarting, i.e. the subspace is full
01659     //
01661     else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
01662 
01663             // Update the recycled subspace even if we have hit the maximum number of restarts.
01664  
01665             // Update the linear problem.
01666             RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01667             problem_->updateSolution( update, true );
01668 
01669             buildRecycleSpace2(gcrodr_iter);
01670 
01671             printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01672 
01673             // NOTE:  If we have hit the maximum number of restarts then quit
01674             if ( numRestarts >= maxRestarts_ ) {
01675               isConverged = false;
01676               break; // break from while(1){gcrodr_iter->iterate()}
01677             }
01678             numRestarts++;
01679 
01680             printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
01681 
01682             // Create the restart vector (first block in the current Krylov basis)
01683             problem_->computeCurrPrecResVec( &*r_ );
01684             index.resize( 1 ); index[0] = 0;
01685             RCP<MV> v0 =  MVT::CloneViewNonConst( *V_,  index );
01686             MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
01687 
01688             // Set the new state and initialize the solver.
01689             GCRODRIterState<ScalarType,MV> restartState;
01690             index.resize( numBlocks_+1 );
01691             for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01692             restartState.V  = MVT::CloneViewNonConst( *V_,  index );
01693             index.resize( keff );
01694             for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01695             restartState.U  = MVT::CloneViewNonConst( *U_,  index );
01696             restartState.C  = MVT::CloneViewNonConst( *C_,  index );
01697             restartState.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
01698             restartState.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
01699             restartState.curDim = 0;
01700             gcrodr_iter->initialize(restartState);
01701 
01702 
01703     } // end of restarting
01704     
01706     //
01707     // we returned from iterate(), but none of our status tests Passed.
01708     // something is wrong, and it is probably our fault.
01709     //
01711     
01712     else {
01713       TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::GCRODRSolMgr::solve(): Invalid return from GCRODRIter::iterate().");
01714     }
01715   }
01716         catch (const GCRODRIterOrthoFailure &e) {
01717     // Try to recover the most recent least-squares solution
01718     gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
01719     
01720     // Check to see if the most recent least-squares solution yielded convergence.
01721     sTest_->checkStatus( &*gcrodr_iter );
01722     if (convTest_->getStatus() != Passed)
01723       isConverged = false;
01724     break;
01725         }
01726         catch (const std::exception &e) {
01727     printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration " 
01728                              << gcrodr_iter->getNumIters() << std::endl 
01729            << e.what() << std::endl;
01730           throw;
01731   }
01732       }
01733      
01734       // Compute the current solution.
01735       // Update the linear problem.
01736       RCP<MV> update = gcrodr_iter->getCurrentUpdate();
01737       problem_->updateSolution( update, true );
01738 
01739       // Inform the linear problem that we are finished with this block linear system.
01740       problem_->setCurrLS();
01741 
01742       // Update indices for the linear systems to be solved.
01743       numRHS2Solve -= 1;
01744       if ( numRHS2Solve > 0 ) {
01745   currIdx[0]++;
01746 
01747         // Set the next indices.
01748         problem_->setLSIndex( currIdx );
01749       }
01750       else {
01751         currIdx.resize( numRHS2Solve );
01752       }
01753 
01754       // If we didn't build a recycle space this solve but ran at least k iterations,
01755       // force build of new recycle space
01756 
01757       if (!builtRecycleSpace_) {
01758         buildRecycleSpace2(gcrodr_iter);
01759         printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
01760       }
01761 
01762     }// while ( numRHS2Solve > 0 )
01763     
01764   }
01765   
01766   // print final summary
01767   sTest_->print( printer_->stream(FinalSummary) );
01768   
01769   // print timing information
01770 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01771   // Calling summarize() can be expensive, so don't call unless the
01772   // user wants to print out timing details.  summarize() will do all
01773   // the work even if it's passed a "black hole" output stream.
01774   if (verbosity_ & TimingDetails)
01775     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01776 #endif
01777  
01778   // get iteration information for this solve
01779   numIters_ = maxIterTest_->getNumIters();
01780  
01781   if (!isConverged) {
01782     return Unconverged; // return from GCRODRSolMgr::solve() 
01783   }
01784   return Converged; // return from GCRODRSolMgr::solve() 
01785 }
01786 
01787 //  Given existing recycle space and Krylov space, build new recycle space
01788 template<class ScalarType, class MV, class OP>
01789 void GCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter) {
01790 
01791   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01792   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01793 
01794   std::vector<MagnitudeType> d(keff);
01795   std::vector<int> index(numBlocks_+1);
01796 
01797   // Get the state
01798   GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
01799   int p = oldState.curDim;
01800 
01801   // insufficient new information to update recycle space
01802   if (p<1) return;
01803 
01804   // Take the norm of the recycled vectors
01805   {
01806     index.resize(keff);
01807     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01808     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_, index );
01809     d.resize(keff);
01810     MVT::MvNorm( *Utmp, d );
01811     for (int i=0; i<keff; ++i) {
01812       d[i] = one / d[i];
01813     }
01814     MVT::MvScale( *Utmp, d );
01815   }
01816 
01817   // Get view into current "full" upper Hessnburg matrix
01818   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
01819 
01820   // Insert D into the leading keff x keff  block of H2
01821   for (int i=0; i<keff; ++i) {
01822     (*H2tmp)(i,i) = d[i];
01823   }
01824 
01825   // Compute the harmoic Ritz pairs for the generalized eigenproblem
01826   // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
01827   int keff_new;
01828   {
01829     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
01830     keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, PPtmp );
01831   }
01832 
01833   // Code to form new U, C
01834   // U = [U V(:,1:p)] * P; (in two steps)
01835 
01836   // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
01837   Teuchos::RCP<MV> U1tmp;
01838   {
01839     index.resize( keff );
01840     for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
01841     Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
01842     index.resize( keff_new );
01843     for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
01844     U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
01845     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
01846     MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
01847   }
01848 
01849   // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
01850   {
01851     index.resize(p);
01852     for (int ii=0; ii < p; ii++) { index[ii] = ii; }
01853     Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01854     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
01855     MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
01856   }
01857 
01858   // Form HP = H*P
01859   Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
01860   {
01861     Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
01862     HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
01863   }
01864 
01865   // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
01866   int info = 0, lwork = -1;
01867   tau_.resize(keff_new);
01868   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01869   TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
01870 
01871   lwork = (int)work_[0];
01872   work_.resize(lwork);
01873   lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01874   TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
01875 
01876   // Explicitly construct Q and R factors
01877   // NOTE:  The upper triangular part of HP is copied into R and HP becomes Q.
01878   Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
01879   for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
01880   lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
01881   TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _ORGQR failed to construct the Q factor.");
01882 
01883   // Form orthonormalized C and adjust U accordingly so that C = A*U
01884   // C = [C V] * Q;
01885 
01886   // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
01887   {
01888     Teuchos::RCP<MV> C1tmp;
01889     {
01890       index.resize(keff);
01891       for (int i=0; i < keff; i++) { index[i] = i; }
01892       Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_,  index );
01893       index.resize(keff_new);
01894       for (int i=0; i < keff_new; i++) { index[i] = i; }
01895       C1tmp  = MVT::CloneViewNonConst( *C1_,  index );
01896       Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
01897       MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
01898     }
01899     // Now compute C += V(:,1:p+1) * Q
01900     {
01901       index.resize( p+1 );
01902       for (int i=0; i < p+1; ++i) { index[i] = i; }
01903       Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
01904       Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
01905       MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
01906     }
01907   }
01908 
01909   // C_ = C1_; (via a swap)
01910   std::swap(C_, C1_);
01911 
01912   // Finally, compute U_ = U_*R^{-1}
01913   // First, compute LU factorization of R
01914   ipiv_.resize(Rtmp.numRows());
01915   lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
01916   TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
01917 
01918   // Now, form inv(R)
01919   lwork = Rtmp.numRows();
01920   work_.resize(lwork);
01921   lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
01922   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
01923 
01924   {
01925     index.resize(keff_new);
01926     for (int i=0; i < keff_new; i++) { index[i] = i; }
01927     Teuchos::RCP<MV> Utmp  = MVT::CloneViewNonConst( *U_,  index );
01928     MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
01929   }
01930 
01931   // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
01932   if (keff != keff_new) {
01933     keff = keff_new;
01934     gcrodr_iter->setSize( keff, numBlocks_ );
01935     // Important to zero this out before next cyle
01936     Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
01937     b1.putScalar(zero);
01938   }
01939 
01940 }
01941 
01942 
01943 //  Compute the harmonic eigenpairs of the projected, dense system.
01944 template<class ScalarType, class MV, class OP>
01945 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs1(int m, 
01946                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
01947                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
01948   int i, j;
01949   bool xtraVec = false;
01950   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01951   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01952 
01953   // Real and imaginary eigenvalue components
01954   std::vector<MagnitudeType> wr(m), wi(m);
01955 
01956   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
01957   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,false);
01958 
01959   // Magnitude of harmonic Ritz values
01960   std::vector<MagnitudeType> w(m);
01961 
01962   // Sorted order of harmonic Ritz values, also used for DGEEV
01963   std::vector<int> iperm(m);
01964 
01965   // Size of workspace and workspace for DGEEV
01966   int lwork = 4*m;
01967   std::vector<ScalarType> work(lwork);
01968 
01969   // Output info
01970   int info = 0;
01971 
01972   // Set flag indicating recycle space has been generated this solve
01973   builtRecycleSpace_ = true;
01974 
01975   // Solve linear system:  H_m^{-H}*e_m 
01976   Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
01977   Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
01978   e_m[m-1] = one;
01979   lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
01980   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01981 
01982   // Compute H_m + d*H_m^{-H}*e_m*e_m^H
01983   ScalarType d = HH(m, m-1) * HH(m, m-1);
01984   Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( HH );
01985   for( i=0; i<m; ++i )
01986     harmHH(i, m-1) += d * e_m[i];
01987 
01988   // Revise to do query for optimal workspace first
01989   // Create simple storage for the left eigenvectors, which we don't care about.
01990   const int ldvl = m;
01991   ScalarType* vl = 0;
01992   lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
01993         vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &info);
01994   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
01995 
01996   // Construct magnitude of each harmonic Ritz value
01997   for( i=0; i<m; ++i )
01998     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
01999 
02000   // Construct magnitude of each harmonic Ritz value
02001   this->sort(w, m, iperm);
02002 
02003   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
02004   if (wi[iperm[recycledBlocks_-1]] != zero) {
02005     int countImag = 0;
02006     for ( i=0; i<recycledBlocks_; ++i ) {
02007       if (wi[iperm[i]] != zero)
02008   countImag++;
02009     }
02010     // Check to see if this count is even or odd:
02011     if (countImag % 2)
02012       xtraVec = true;
02013   }
02014 
02015   // Select recycledBlocks_ smallest eigenvectors
02016   for( i=0; i<recycledBlocks_; ++i ) {
02017     for( j=0; j<m; j++ ) {
02018       PP(j,i) = vr(j,iperm[i]);
02019     }
02020   }
02021   
02022   if (xtraVec) { // we need to store one more vector
02023     if (wi[iperm[recycledBlocks_-1]] > 0) { // I picked the "real" component
02024       for( j=0; j<m; ++j ) {   // so get the "imag" component
02025   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
02026       }
02027     }
02028     else { //  I picked the "imag" component
02029       for( j=0; j<m; ++j ) {   // so get the "real" component
02030   PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
02031       }
02032     }
02033   }
02034 
02035   // Return whether we needed to store an additional vector
02036   if (xtraVec) {
02037     return recycledBlocks_+1;
02038   }
02039   return recycledBlocks_;
02040 }
02041 
02042 //  Compute the harmonic eigenpairs of the projected, dense system.
02043 template<class ScalarType, class MV, class OP>
02044 int GCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecs2(int keff, int m, 
02045                  const Teuchos::SerialDenseMatrix<int,ScalarType>& HH, 
02046                  const Teuchos::RCP<const MV>& VV,
02047                  Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
02048   int i, j;
02049   int m2 = HH.numCols();
02050   bool xtraVec = false;
02051   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
02052   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
02053   std::vector<int> index;
02054   
02055   // Real and imaginary eigenvalue components
02056   std::vector<ScalarType> wr(m2), wi(m2);
02057 
02058   // Magnitude of harmonic Ritz values
02059   std::vector<MagnitudeType> w(m2);
02060 
02061   // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
02062   Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,false);
02063 
02064   // Sorted order of harmonic Ritz values
02065   std::vector<int> iperm(m2);
02066 
02067   // Set flag indicating recycle space has been generated this solve
02068   builtRecycleSpace_ = true;
02069 
02070   // Form matrices for generalized eigenproblem
02071   
02072   // B = H2' * H2; Don't zero out matrix when constructing
02073   Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,false);
02074   B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
02075   
02076   // A_tmp = | C'*U        0 |
02077   //         | V_{m+1}'*U  I |
02078   Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keff+m+1, keff+m );
02079 
02080   // A_tmp(1:keff,1:keff) = C' * U;
02081   index.resize(keff);
02082   for (int i=0; i<keff; ++i) { index[i] = i; }
02083   Teuchos::RCP<const MV> Ctmp  = MVT::CloneView( *C_, index );
02084   Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_, index );
02085   Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keff, keff );
02086   MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
02087 
02088   // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
02089   Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keff, keff );
02090   index.resize(m+1);
02091   for (i=0; i < m+1; i++) { index[i] = i; }
02092   Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
02093   MVT::MvTransMv( one, *Vp, *Utmp, A21 );
02094 
02095   // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
02096   for( i=keff; i<keff+m; i++ ) {
02097     A_tmp(i,i) = one;
02098   }
02099 
02100   // A = H2' * A_tmp;
02101   Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
02102   A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
02103 
02104   // Compute k smallest harmonic Ritz pairs
02105   // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
02106   //                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
02107   //                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
02108   //                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
02109   // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
02110   char balanc='P', jobvl='N', jobvr='V', sense='N'; 
02111   int ld = A.numRows();
02112   int lwork = 6*ld;
02113   int ldvl = ld, ldvr = ld;
02114   int info = 0,ilo = 0,ihi = 0;
02115   ScalarType abnrm = zero, bbnrm = zero;
02116   ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
02117   std::vector<ScalarType> beta(ld);
02118   std::vector<ScalarType> work(lwork);
02119   std::vector<MagnitudeType> lscale(ld), rscale(ld);  
02120   std::vector<MagnitudeType> rconde(ld), rcondv(ld);
02121   std::vector<int> iwork(ld+6);
02122   int *bwork = 0; // If sense == 'N', bwork is never referenced
02123   lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0], 
02124                &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0], 
02125                &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
02126   TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
02127   
02128   // Construct magnitude of each harmonic Ritz value
02129   // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
02130   for( i=0; i<ld; i++ ) {
02131     w[i] = Teuchos::ScalarTraits<ScalarType>::squareroot( (wr[i]/beta[i])*(wr[i]/beta[i]) + (wi[i]/beta[i])*(wi[i]/beta[i]) );
02132   }
02133 
02134   // Construct magnitude of each harmonic Ritz value
02135   this->sort(w,ld,iperm);
02136 
02137   // Determine exact size for PP (i.e., determine if we need to store an additional vector)
02138   if (wi[iperm[ld-recycledBlocks_]] != zero) {
02139     int countImag = 0;
02140     for ( i=ld-recycledBlocks_; i<ld; i++ ) {
02141       if (wi[iperm[i]] != zero)
02142   countImag++;
02143     }
02144     // Check to see if this count is even or odd:
02145     if (countImag % 2)
02146       xtraVec = true;
02147   }
02148   
02149   // Select recycledBlocks_ smallest eigenvectors
02150   for( i=0; i<recycledBlocks_; i++ ) {
02151     for( j=0; j<ld; j++ ) {
02152       PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
02153     }
02154   }
02155   
02156   if (xtraVec) { // we need to store one more vector
02157     if (wi[iperm[ld-recycledBlocks_]] > 0) { // I picked the "real" component
02158       for( j=0; j<ld; j++ ) {   // so get the "imag" component
02159   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
02160       }
02161     }
02162     else { // I picked the "imag" component
02163       for( j=0; j<ld; j++ ) {   // so get the "real" component
02164   PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
02165       }
02166     }
02167   }
02168   
02169   // Return whether we needed to store an additional vector
02170   if (xtraVec) {
02171     return recycledBlocks_+1;
02172   }
02173   return recycledBlocks_;
02174 }
02175 
02176 
02177 // This method sorts list of n floating-point numbers and return permutation vector
02178 template<class ScalarType, class MV, class OP>
02179 void GCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) {
02180   int l, r, j, i, flag;
02181   int    RR2;
02182   double dRR, dK;
02183   
02184   // Initialize the permutation vector.
02185   for(j=0;j<n;j++)
02186     iperm[j] = j;
02187   
02188   if (n <= 1) return;
02189   
02190   l    = n / 2 + 1;
02191   r    = n - 1;
02192   l    = l - 1;
02193   dRR  = dlist[l - 1];
02194   dK   = dlist[l - 1];
02195   
02196   RR2 = iperm[l - 1];
02197   while (r != 0) {
02198     j = l;
02199     flag = 1;
02200     
02201     while (flag == 1) {
02202       i = j;
02203       j = j + j;
02204       
02205       if (j > r + 1)
02206   flag = 0;
02207       else {
02208   if (j < r + 1)
02209     if (dlist[j] > dlist[j - 1]) j = j + 1;
02210   
02211   if (dlist[j - 1] > dK) {
02212     dlist[i - 1] = dlist[j - 1];
02213     iperm[i - 1] = iperm[j - 1];
02214   }
02215   else {
02216     flag = 0;
02217   }
02218       }
02219     }
02220     dlist[i - 1] = dRR;
02221     iperm[i - 1] = RR2;
02222     
02223     if (l == 1) {
02224       dRR  = dlist [r];
02225       RR2 = iperm[r];
02226       dK = dlist[r];
02227       dlist[r] = dlist[0];
02228       iperm[r] = iperm[0];
02229       r = r - 1;
02230     }
02231     else {
02232       l   = l - 1;
02233       dRR  = dlist[l - 1];
02234       RR2  = iperm[l - 1];
02235       dK   = dlist[l - 1];
02236     }
02237   }
02238   dlist[0] = dRR;
02239   iperm[0] = RR2; 
02240 }
02241 
02242 
02243 //  This method requires the solver manager to return a string that describes itself.
02244 template<class ScalarType, class MV, class OP>
02245 std::string GCRODRSolMgr<ScalarType,MV,OP>::description() const {
02246   std::ostringstream oss;
02247   oss << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
02248   oss << "{";
02249   oss << "Ortho Type='"<<orthoType_;
02250   oss << ", Num Blocks=" <<numBlocks_;
02251   oss << ", Num Recycle Blocks=" << recycledBlocks_;
02252   oss << ", Max Restarts=" << maxRestarts_;
02253   oss << "}";
02254   return oss.str();
02255 }
02256 
02257 
02258 } // end Belos namespace
02259 
02260 #endif /* BELOS_GCRODR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines