BelosRCGSolMgr.hpp

Go to the documentation of this file.
00001 // ***********************************************************************
00002 //
00003 //                 Belos: Block Linear Solvers Package
00004 //                 Copyright (2004) Sandia Corporation
00005 //
00006 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00007 // license for use of this work by or on behalf of the U.S. Government.
00008 //
00009 // This library is free software; you can redistribute it and/or modify
00010 // it under the terms of the GNU Lesser General Public License as
00011 // published by the Free Software Foundation; either version 2.1 of the
00012 // License, or (at your option) any later version.
00013 //
00014 // This library is distributed in the hope that it will be useful, but
00015 // WITHOUT ANY WARRANTY; without even the implied warranty of
00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 // Lesser General Public License for more details.
00018 //
00019 // You should have received a copy of the GNU Lesser General Public
00020 // License along with this library; if not, write to the Free Software
00021 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00022 // USA
00023 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00024 //
00025 // ***********************************************************************
00026 // @HEADER
00027 
00028 #ifndef BELOS_RCG_SOLMGR_HPP
00029 #define BELOS_RCG_SOLMGR_HPP
00030 
00035 #include "BelosConfigDefs.hpp"
00036 #include "BelosTypes.hpp"
00037 
00038 #include "BelosLinearProblem.hpp"
00039 #include "BelosSolverManager.hpp"
00040 
00041 #include "BelosRCGIter.hpp"
00042 #include "BelosDGKSOrthoManager.hpp"
00043 #include "BelosICGSOrthoManager.hpp"
00044 #include "BelosIMGSOrthoManager.hpp"
00045 #include "BelosStatusTestMaxIters.hpp"
00046 #include "BelosStatusTestGenResNorm.hpp"
00047 #include "BelosStatusTestCombo.hpp"
00048 #include "BelosStatusTestOutput.hpp"
00049 #include "BelosOutputManager.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #include "Teuchos_LAPACK.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00063 namespace Belos {
00064   
00066 
00067 
00074   class RCGSolMgrLinearProblemFailure : public BelosError {public:
00075     RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00076     {}};
00077  
00084   class RCGSolMgrLAPACKFailure : public BelosError {public:
00085     RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00086     {}};
00087 
00094   class RCGSolMgrRecyclingFailure : public BelosError {public:
00095     RCGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
00096     {}};
00097 
00099   
00100   template<class ScalarType, class MV, class OP>
00101   class RCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00102     
00103   private:
00104     typedef MultiVecTraits<ScalarType,MV> MVT;
00105     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00106     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00107     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00108     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00109     
00110   public:
00111     
00113 
00114 
00120      RCGSolMgr();
00121 
00142     RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00143        const Teuchos::RCP<Teuchos::ParameterList> &pl );
00144     
00146     virtual ~RCGSolMgr() {};
00148     
00150 
00151     
00152     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00153       return *problem_;
00154     }
00155 
00157     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00158     
00160     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00161     
00167     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00168       return tuple(timerSolve_);
00169     }
00170 
00172     int getNumIters() const {
00173       return numIters_;
00174     }
00175 
00177     bool isLOADetected() const { return false; }
00178  
00180     
00182 
00183    
00185     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00186    
00188     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00189     
00191    
00193 
00194 
00198     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00199 
00203     void resetRecycleSpace( ) { existU_ = false; }
00205  
00207 
00208     
00226     ReturnType solve();
00227     
00229     
00232     
00234     std::string description() const;
00235     
00237     
00238   private:
00239 
00240     // Called by all constructors; Contains init instructions common to all constructors
00241     void init();
00242  
00243     //  Computes harmonic eigenpairs of projected matrix created during one cycle.
00244     //  Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude.
00245     void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F, 
00246                          const Teuchos::SerialDenseMatrix<int,ScalarType> &G, 
00247                          Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
00248  
00249     // Sort list of n floating-point numbers and return permutation vector
00250     void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
00251 
00252     // Initialize solver state storage
00253     void initializeStateStorage();
00254 
00255     // Linear problem.
00256     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00257     
00258     // Output manager.
00259     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00260     Teuchos::RCP<std::ostream> outputStream_;
00261 
00262     // Status test.
00263     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00264     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00265     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00266     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00267 
00268     // Current parameter list.
00269     Teuchos::RCP<ParameterList> params_;
00270     
00271     // Default solver values.
00272     static const MagnitudeType convtol_default_;
00273     static const int maxIters_default_;
00274     static const int numBlocks_default_;
00275     static const int recycleBlocks_default_;
00276     static const bool showMaxResNormOnly_default_;
00277     static const int verbosity_default_;
00278     static const int outputFreq_default_;
00279     static const std::string label_default_;
00280     static const Teuchos::RCP<std::ostream> outputStream_default_;
00281 
00282     // Current solver values.
00283     MagnitudeType convtol_;
00284     int maxIters_, numIters_;
00285     int numBlocks_, recycleBlocks_;
00286     bool showMaxResNormOnly_;
00287     int verbosity_, outputFreq_;
00288 
00290     // Solver State Storage
00292     // Search vectors
00293     Teuchos::RCP<MV> P_;
00294     //
00295     // A times current search direction
00296     Teuchos::RCP<MV> Ap_;
00297     //
00298     // Residual vector
00299     Teuchos::RCP<MV> r_;
00300     //
00301     // Preconditioned residual
00302     Teuchos::RCP<MV> z_;
00303     //
00304     // Flag indicating that the recycle space should be used
00305     bool existU_;
00306     //
00307     // Flag indicating that the updated recycle space has been created
00308     bool existU1_;
00309     //
00310     // Recycled subspace and its image
00311     Teuchos::RCP<MV> U_, AU_;
00312     //
00313     // Recycled subspace for next system and its image
00314     Teuchos::RCP<MV> U1_;
00315     //
00316     // Coefficients arising in RCG iteration
00317     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
00318     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
00319     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
00320     //
00321     // Solutions to local least-squares problems
00322     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
00323     //
00324     // The matrix U^T A U
00325     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
00326     //
00327     // LU factorization of U^T A U
00328     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
00329     //
00330     // Data from LU factorization of UTAU
00331     Teuchos::RCP<std::vector<int> > ipiv_;
00332     //
00333     // The matrix (AU)^T AU
00334     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
00335     //
00336     // The scalar r'*z
00337     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
00338     //
00339     // Matrices needed for calculation of harmonic Ritz eigenproblem
00340     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
00341     //
00342     // Matrices needed for updating recycle space
00343     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
00344     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
00345     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
00346     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
00347     Teuchos::RCP<MV> U1Y1_, PY2_;
00348     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
00349     ScalarType dold;
00351     
00352     // Timers.
00353     std::string label_;
00354     Teuchos::RCP<Teuchos::Time> timerSolve_;
00355 
00356     // Internal state variables.
00357     bool params_Set_;
00358   };
00359 
00360 
00361 // Default solver values.
00362 template<class ScalarType, class MV, class OP>
00363 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType RCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00364 
00365 template<class ScalarType, class MV, class OP>
00366 const int RCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00367 
00368 template<class ScalarType, class MV, class OP>
00369 const int RCGSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 25;
00370  
00371 template<class ScalarType, class MV, class OP>
00372 const int RCGSolMgr<ScalarType,MV,OP>::recycleBlocks_default_ = 3;
00373 
00374 template<class ScalarType, class MV, class OP>
00375 const bool RCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00376 
00377 template<class ScalarType, class MV, class OP>
00378 const int RCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00379 
00380 template<class ScalarType, class MV, class OP>
00381 const int RCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00382 
00383 template<class ScalarType, class MV, class OP>
00384 const std::string RCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00385 
00386 template<class ScalarType, class MV, class OP>
00387 const Teuchos::RCP<std::ostream> RCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00388 
00389 // Empty Constructor
00390 template<class ScalarType, class MV, class OP>
00391 RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr() {
00392   init();
00393 }
00394 
00395 // Basic Constructor
00396 template<class ScalarType, class MV, class OP>
00397 RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr( 
00398                  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00399                  const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00400   problem_(problem)
00401 {
00402   init();
00403   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00404 
00405   // If the parameter list pointer is null, then set the current parameters to the default parameter list.
00406   if ( !is_null(pl) ) {
00407     setParameters( pl );  
00408   }
00409 }
00410 
00411 // Common instructions executed in all constructors
00412 template<class ScalarType, class MV, class OP>
00413 void RCGSolMgr<ScalarType,MV,OP>::init()
00414 {
00415   outputStream_ = outputStream_default_;
00416   convtol_ = convtol_default_;
00417   maxIters_ = maxIters_default_;
00418   numBlocks_ = numBlocks_default_;
00419   recycleBlocks_ = recycleBlocks_default_;
00420   verbosity_ = verbosity_default_;
00421   outputFreq_= outputFreq_default_;
00422   showMaxResNormOnly_ = showMaxResNormOnly_default_;
00423   label_ = label_default_;
00424   params_Set_ = false;
00425   P_ = Teuchos::null;
00426   Ap_ = Teuchos::null;
00427   r_ = Teuchos::null;
00428   z_ = Teuchos::null;
00429   existU_ = false;
00430   existU1_ = false;
00431   U_ = Teuchos::null;
00432   AU_ = Teuchos::null;
00433   U1_ = Teuchos::null;
00434   Alpha_ = Teuchos::null;
00435   Beta_ = Teuchos::null;
00436   D_ = Teuchos::null;
00437   Delta_ = Teuchos::null;
00438   UTAU_ = Teuchos::null;
00439   LUUTAU_ = Teuchos::null;
00440   ipiv_ = Teuchos::null;
00441   AUTAU_ = Teuchos::null;
00442   rTz_old_ = Teuchos::null;
00443   F_ = Teuchos::null;
00444   G_ = Teuchos::null;
00445   Y_ = Teuchos::null;
00446   L2_ = Teuchos::null;
00447   DeltaL2_ = Teuchos::null;
00448   AU1TUDeltaL2_ = Teuchos::null;
00449   AU1TAU1_ = Teuchos::null;
00450   AU1TU1_ = Teuchos::null;
00451   AU1TAP_ = Teuchos::null;
00452   FY_ = Teuchos::null;
00453   GY_ = Teuchos::null;
00454   APTAP_ = Teuchos::null;
00455   U1Y1_ = Teuchos::null;
00456   PY2_ = Teuchos::null;
00457   AUTAP_ = Teuchos::null;
00458   AU1TU_ = Teuchos::null;
00459   dold = 0.;
00460 }
00461 
00462 template<class ScalarType, class MV, class OP>
00463 void RCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00464 {
00465   // Create the internal parameter list if ones doesn't already exist.
00466   if (params_ == Teuchos::null) {
00467     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00468   }
00469   else {
00470     params->validateParameters(*getValidParameters());
00471   }
00472 
00473   // Check for maximum number of iterations
00474   if (params->isParameter("Maximum Iterations")) {
00475     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00476 
00477     // Update parameter in our list and in status test.
00478     params_->set("Maximum Iterations", maxIters_);
00479     if (maxIterTest_!=Teuchos::null)
00480       maxIterTest_->setMaxIters( maxIters_ );
00481   }
00482 
00483   // Check for the maximum number of blocks.
00484   if (params->isParameter("Num Blocks")) {
00485     numBlocks_ = params->get("Num Blocks",numBlocks_default_);
00486     TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
00487                        "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
00488  
00489     // Update parameter in our list.
00490     params_->set("Num Blocks", numBlocks_);
00491   }
00492  
00493   // Check for the maximum number of blocks.
00494   if (params->isParameter("Num Recycled Blocks")) {
00495     recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_);
00496     TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
00497                        "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
00498  
00499     TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
00500                        "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
00501  
00502     // Update parameter in our list.
00503     params_->set("Num Recycled Blocks", recycleBlocks_);
00504   }
00505 
00506   // Check to see if the timer label changed.
00507   if (params->isParameter("Timer Label")) {
00508     std::string tempLabel = params->get("Timer Label", label_default_);
00509 
00510     // Update parameter in our list and solver timer
00511     if (tempLabel != label_) {
00512       label_ = tempLabel;
00513       params_->set("Timer Label", label_);
00514       std::string solveLabel = label_ + ": RCGSolMgr total solve time";
00515       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00516     }
00517   }
00518 
00519   // Check for a change in verbosity level
00520   if (params->isParameter("Verbosity")) {
00521     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00522       verbosity_ = params->get("Verbosity", verbosity_default_);
00523     } else {
00524       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00525     }
00526 
00527     // Update parameter in our list.
00528     params_->set("Verbosity", verbosity_);
00529     if (printer_ != Teuchos::null)
00530       printer_->setVerbosity(verbosity_);
00531   }
00532 
00533   // output stream
00534   if (params->isParameter("Output Stream")) {
00535     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00536 
00537     // Update parameter in our list.
00538     params_->set("Output Stream", outputStream_);
00539     if (printer_ != Teuchos::null)
00540       printer_->setOStream( outputStream_ );
00541   }
00542 
00543   // frequency level
00544   if (verbosity_ & Belos::StatusTestDetails) {
00545     if (params->isParameter("Output Frequency")) {
00546       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00547     }
00548 
00549     // Update parameter in out list and output status test.
00550     params_->set("Output Frequency", outputFreq_);
00551     if (outputTest_ != Teuchos::null)
00552       outputTest_->setOutputFrequency( outputFreq_ );
00553   }
00554 
00555   // Create output manager if we need to.
00556   if (printer_ == Teuchos::null) {
00557     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00558   }  
00559   
00560   // Convergence
00561   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00562   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00563 
00564   // Check for convergence tolerance
00565   if (params->isParameter("Convergence Tolerance")) {
00566     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00567 
00568     // Update parameter in our list and residual tests.
00569     params_->set("Convergence Tolerance", convtol_);
00570     if (convTest_ != Teuchos::null)
00571       convTest_->setTolerance( convtol_ );
00572   }
00573   
00574   if (params->isParameter("Show Maximum Residual Norm Only")) {
00575     showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00576 
00577     // Update parameter in our list and residual tests
00578     params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00579     if (convTest_ != Teuchos::null)
00580       convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00581   }
00582 
00583   // Create status tests if we need to.
00584 
00585   // Basic test checks maximum iterations and native residual.
00586   if (maxIterTest_ == Teuchos::null)
00587     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00588   
00589   // Implicit residual test, using the native residual to determine if convergence was achieved.
00590   if (convTest_ == Teuchos::null)
00591     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00592   
00593   if (sTest_ == Teuchos::null)
00594     sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00595   
00596   if (outputTest_ == Teuchos::null) {
00597     if (outputFreq_ > 0) {
00598       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00599                     sTest_, 
00600                     outputFreq_, 
00601                     Passed+Failed+Undefined ) ); 
00602     }
00603     else {
00604       outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_, 
00605                     sTest_, 1 ) );
00606     }
00607   }
00608 
00609   // Create the timer if we need to.
00610   if (timerSolve_ == Teuchos::null) {
00611     std::string solveLabel = label_ + ": RCGSolMgr total solve time";
00612     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00613   }
00614 
00615   // Inform the solver manager that the current parameters were set.
00616   params_Set_ = true;
00617 }
00618 
00619     
00620 template<class ScalarType, class MV, class OP>
00621 Teuchos::RCP<const Teuchos::ParameterList> 
00622 RCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00623 {
00624   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00625   
00626   // Set all the valid parameters and their default values.
00627   if(is_null(validPL)) {
00628     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00629     pl->set("Convergence Tolerance", convtol_default_,
00630       "The relative residual tolerance that needs to be achieved by the\n"
00631       "iterative solver in order for the linear system to be declared converged.");
00632     pl->set("Maximum Iterations", maxIters_default_,
00633       "The maximum number of block iterations allowed for each\n"
00634       "set of RHS solved.");
00635     pl->set("Num Blocks", numBlocks_default_,
00636       "The length of a cycle (and this max number of search vectors kept)\n");
00637     pl->set("Num Recycled Blocks", recycleBlocks_default_,
00638       "The number of vectors in the recycle subspace.");
00639     pl->set("Verbosity", verbosity_default_,
00640       "What type(s) of solver information should be outputted\n"
00641       "to the output stream.");
00642     pl->set("Output Frequency", outputFreq_default_,
00643       "How often convergence information should be outputted\n"
00644       "to the output stream.");  
00645     pl->set("Output Stream", outputStream_default_,
00646       "A reference-counted pointer to the output stream where all\n"
00647       "solver output is sent.");
00648     pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00649       "When convergence information is printed, only show the maximum\n"
00650       "relative residual norm when the block size is greater than one.");
00651     pl->set("Timer Label", label_default_,
00652       "The string to use as a prefix for the timer labels.");
00653     //  pl->set("Restart Timers", restartTimers_);
00654     validPL = pl;
00655   }
00656   return validPL;
00657 }
00658 
00659 // initializeStateStorage
00660 template<class ScalarType, class MV, class OP>
00661 void RCGSolMgr<ScalarType,MV,OP>::initializeStateStorage() {
00662 
00663     // Check if there is any multivector to clone from.
00664     Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
00665     if (rhsMV == Teuchos::null) {
00666       // Nothing to do
00667       return;
00668     }
00669     else {
00670 
00671       // Initialize the state storage
00672       TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
00673                          "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
00674 
00675       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00676       if (P_ == Teuchos::null) {
00677         P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
00678       }
00679       else {
00680         // Generate P_ by cloning itself ONLY if more space is needed.
00681         if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
00682           Teuchos::RCP<const MV> tmp = P_;
00683           P_ = MVT::Clone( *tmp, numBlocks_+2 );
00684         }
00685       }
00686 
00687       // Generate Ap_ only if it doesn't exist
00688       if (Ap_ == Teuchos::null)
00689         Ap_ = MVT::Clone( *rhsMV, 1 );
00690 
00691       // Generate r_ only if it doesn't exist
00692       if (r_ == Teuchos::null)
00693         r_ = MVT::Clone( *rhsMV, 1 );
00694 
00695       // Generate z_ only if it doesn't exist
00696       if (z_ == Teuchos::null)
00697         z_ = MVT::Clone( *rhsMV, 1 );
00698 
00699       // If the recycle space has not been initialized before, generate it using the RHS from lp_.
00700       if (U_ == Teuchos::null) {
00701         U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
00702       }
00703       else {
00704         // Generate U_ by cloning itself ONLY if more space is needed.
00705         if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
00706           Teuchos::RCP<const MV> tmp = U_;
00707           U_ = MVT::Clone( *tmp, recycleBlocks_ );
00708         }
00709       }
00710 
00711       // If the recycle space has not be initialized before, generate it using the RHS from lp_.
00712       if (AU_ == Teuchos::null) {
00713         AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
00714       }
00715       else {
00716         // Generate AU_ by cloning itself ONLY if more space is needed.
00717         if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
00718           Teuchos::RCP<const MV> tmp = AU_;
00719           AU_ = MVT::Clone( *tmp, recycleBlocks_ );
00720         }
00721       }
00722 
00723       // If the recycle space has not been initialized before, generate it using the RHS from lp_.
00724       if (U1_ == Teuchos::null) {
00725         U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
00726       }
00727       else {
00728         // Generate U1_ by cloning itself ONLY if more space is needed.
00729         if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
00730           Teuchos::RCP<const MV> tmp = U1_;
00731           U1_ = MVT::Clone( *tmp, recycleBlocks_ );
00732         }
00733       }
00734 
00735       // Generate Alpha_ only if it doesn't exist, otherwise resize it.
00736       if (Alpha_ == Teuchos::null)
00737         Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
00738       else {
00739         if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
00740           Alpha_->reshape( numBlocks_, 1 );
00741       }
00742 
00743       // Generate Beta_ only if it doesn't exist, otherwise resize it.
00744       if (Beta_ == Teuchos::null)
00745         Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
00746       else {
00747         if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
00748           Beta_->reshape( numBlocks_ + 1, 1 );
00749       }
00750 
00751       // Generate D_ only if it doesn't exist, otherwise resize it.
00752       if (D_ == Teuchos::null)
00753         D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
00754       else {
00755         if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
00756           D_->reshape( numBlocks_, 1 );
00757       }
00758 
00759       // Generate Delta_ only if it doesn't exist, otherwise resize it.
00760       if (Delta_ == Teuchos::null)
00761         Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
00762       else {
00763         if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
00764           Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
00765       }
00766 
00767       // Generate UTAU_ only if it doesn't exist, otherwise resize it.
00768       if (UTAU_ == Teuchos::null)
00769         UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00770       else {
00771         if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
00772           UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
00773       }
00774 
00775       // Generate LUUTAU_ only if it doesn't exist, otherwise resize it.
00776       if (LUUTAU_ == Teuchos::null)
00777         LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00778       else {
00779         if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
00780           LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
00781       }
00782 
00783       // Generate ipiv_ only if it doesn't exist, otherwise resize it.
00784       if (ipiv_ == Teuchos::null)
00785         ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) );
00786       else {
00787         if ( ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it
00788           ipiv_->resize(recycleBlocks_);
00789       }
00790 
00791       // Generate AUTAU_ only if it doesn't exist, otherwise resize it.
00792       if (AUTAU_ == Teuchos::null)
00793         AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00794       else {
00795         if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
00796           AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
00797       }
00798 
00799       // Generate rTz_old_ only if it doesn't exist
00800       if (rTz_old_ == Teuchos::null)
00801         rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
00802       else {
00803         if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
00804           rTz_old_->reshape( 1, 1 );
00805       }
00806 
00807       // Generate F_ only if it doesn't exist
00808       if (F_ == Teuchos::null)
00809         F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
00810       else {
00811         if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
00812           F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
00813       }
00814 
00815       // Generate G_ only if it doesn't exist
00816       if (G_ == Teuchos::null)
00817         G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
00818       else {
00819         if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
00820           G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
00821       }
00822 
00823       // Generate Y_ only if it doesn't exist
00824       if (Y_ == Teuchos::null)
00825         Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
00826       else {
00827         if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
00828           Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
00829       }
00830 
00831       // Generate L2_ only if it doesn't exist
00832       if (L2_ == Teuchos::null)
00833         L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
00834       else {
00835         if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
00836           L2_->reshape( numBlocks_+1, numBlocks_ );
00837       }
00838 
00839       // Generate DeltaL2_ only if it doesn't exist
00840       if (DeltaL2_ == Teuchos::null)
00841         DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
00842       else {
00843         if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
00844           DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
00845       }
00846 
00847       // Generate AU1TUDeltaL2_ only if it doesn't exist
00848       if (AU1TUDeltaL2_ == Teuchos::null)
00849         AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
00850       else {
00851         if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
00852           AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
00853       }
00854 
00855       // Generate AU1TAU1_ only if it doesn't exist
00856       if (AU1TAU1_ == Teuchos::null)
00857         AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00858       else {
00859         if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
00860           AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
00861       }
00862 
00863       // Generate GY_ only if it doesn't exist
00864       if (GY_ == Teuchos::null)
00865         GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
00866       else {
00867         if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
00868           GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
00869       }
00870 
00871       // Generate AU1TU1_ only if it doesn't exist
00872       if (AU1TU1_ == Teuchos::null)
00873         AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00874       else {
00875         if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
00876           AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
00877       }
00878 
00879       // Generate FY_ only if it doesn't exist
00880       if (FY_ == Teuchos::null)
00881         FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
00882       else {
00883         if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
00884           FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
00885       }
00886 
00887       // Generate AU1TAP_ only if it doesn't exist
00888       if (AU1TAP_ == Teuchos::null)
00889         AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
00890       else {
00891         if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
00892           AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
00893       }
00894 
00895       // Generate APTAP_ only if it doesn't exist
00896       if (APTAP_ == Teuchos::null)
00897         APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
00898       else {
00899         if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
00900           APTAP_->reshape( numBlocks_, numBlocks_ );
00901       }
00902 
00903       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00904       if (U1Y1_ == Teuchos::null) {
00905         U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
00906       }
00907       else {
00908         // Generate U1Y1_ by cloning itself ONLY if more space is needed.
00909         if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
00910           Teuchos::RCP<const MV> tmp = U1Y1_;
00911           U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
00912         }
00913       }
00914 
00915       // If the subspace has not been initialized before, generate it using the RHS from lp_.
00916       if (PY2_ == Teuchos::null) {
00917         PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
00918       }
00919       else {
00920         // Generate PY2_ by cloning itself ONLY if more space is needed.
00921         if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
00922           Teuchos::RCP<const MV> tmp = PY2_;
00923           PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
00924         }
00925       }
00926 
00927       // Generate AUTAP_ only if it doesn't exist
00928       if (AUTAP_ == Teuchos::null)
00929         AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
00930       else {
00931         if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
00932           AUTAP_->reshape( recycleBlocks_, numBlocks_ );
00933       }
00934 
00935       // Generate AU1TU_ only if it doesn't exist
00936       if (AU1TU_ == Teuchos::null)
00937         AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
00938       else {
00939         if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
00940           AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
00941       }
00942 
00943 
00944     }
00945 }
00946 
00947 template<class ScalarType, class MV, class OP>
00948 ReturnType RCGSolMgr<ScalarType,MV,OP>::solve() {
00949 
00950   Teuchos::LAPACK<int,ScalarType> lapack;
00951   std::vector<int> index(recycleBlocks_);
00952   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00953   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00954 
00955   // Count of number of cycles performed on current rhs
00956   int cycle = 0;
00957 
00958   // Set the current parameters if they were not set before.
00959   // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
00960   // then didn't set any parameters using setParameters().
00961   if (!params_Set_) {
00962     setParameters(Teuchos::parameterList(*getValidParameters()));
00963   }
00964 
00965   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
00966                      "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
00967   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
00968                      "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00969 
00970   // Create indices for the linear systems to be solved.
00971   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00972   std::vector<int> currIdx(1);
00973   currIdx[0] = 0;
00974 
00975   // Inform the linear problem of the current linear system to solve.
00976   problem_->setLSIndex( currIdx );
00977 
00978   // Check the number of blocks and change them if necessary.
00979   int dim = MVT::GetVecLength( *(problem_->getRHS()) );
00980   if (numBlocks_ > dim) {
00981     numBlocks_ = dim;
00982     params_->set("Num Blocks", numBlocks_);
00983     printer_->stream(Warnings) <<
00984       "Warning! Requested Krylov subspace dimension is larger that operator dimension!" << endl <<
00985       " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << endl;
00986   }
00987 
00988   // Initialize storage for all state variables
00989   initializeStateStorage();
00990 
00991   // Parameter list
00992   Teuchos::ParameterList plist;
00993   plist.set("Num Blocks",numBlocks_);
00994   plist.set("Recycled Blocks",recycleBlocks_);
00995   
00996   // Reset the status test.  
00997   outputTest_->reset();
00998 
00999   // Assume convergence is achieved, then let any failed convergence set this to false.
01000   bool isConverged = true;  
01001 
01002   // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
01003   if (existU_) {
01004     index.resize(recycleBlocks_);
01005     for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01006     Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_,  index );
01007     index.resize(recycleBlocks_);
01008     for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01009     Teuchos::RCP<MV> AUtmp = MVT::CloneView( *AU_, index );
01010     // Initialize AU
01011     problem_->applyOp( *Utmp, *AUtmp );
01012     // Initialize UTAU
01013     Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
01014     MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
01015     // Initialize AUTAU  ( AUTAU = AU'*(M\AU) )
01016     Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
01017     if ( problem_->getLeftPrec() != Teuchos::null ) {
01018       index.resize(recycleBlocks_);
01019       for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01020       index.resize(recycleBlocks_);
01021       for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01022       Teuchos::RCP<MV> LeftPCAU = MVT::CloneView( *U1_, index ); // use U1 as temp storage
01023       problem_->applyLeftPrec( *AUtmp, *LeftPCAU );
01024       MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
01025     } else {
01026       MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
01027     }
01028   }
01029 
01031   // RCG solver
01032 
01033   Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
01034   rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
01035 
01036   // Enter solve() iterations
01037   {
01038     Teuchos::TimeMonitor slvtimer(*timerSolve_);
01039 
01040     while ( numRHS2Solve > 0 ) {
01041 
01042       // Debugging output to tell use if recycle space exists and will be used
01043       if (printer_->isVerbosity( Debug ) ) {
01044         if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
01045         else printer_->print( Debug, "No recycle space exists." );
01046       }
01047 
01048       // Reset the number of iterations.
01049       rcg_iter->resetNumIters();
01050 
01051       // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
01052       rcg_iter->setSize( recycleBlocks_, numBlocks_ );
01053 
01054       // Reset the number of calls that the status test output knows about.
01055       outputTest_->resetNumCalls();
01056 
01057       // indicate that updated recycle space has not yet been generated for this linear system
01058       existU1_ = false;
01059 
01060       // reset cycle count
01061       cycle = 0;
01062 
01063       // Get the current residual 
01064       problem_->computeCurrResVec( &*r_ );
01065 
01066       // If U exists, find best soln over this space first
01067       if (existU_) {
01068         // Solve linear system UTAU * y = (U'*r)
01069         Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
01070         index.resize(recycleBlocks_);
01071         for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01072         Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_,  index );
01073         MVT::MvTransMv( one, *Utmp, *r_, Utr );
01074         Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
01075         Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
01076         LUUTAUtmp.assign(UTAUtmp);
01077         int info = 0;
01078         lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
01079         TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
01080                            "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
01081 
01082         // Update solution (x = x + U*y)
01083         MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
01084  
01085         // Update residual ( r = r - AU*y )
01086         index.resize(recycleBlocks_);
01087         for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01088         Teuchos::RCP<MV> AUtmp  = MVT::CloneView( *AU_, index );
01089         MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
01090       }
01091 
01092       if ( problem_->getLeftPrec() != Teuchos::null ) {
01093         problem_->applyLeftPrec( *r_, *z_ );
01094       } else {
01095         z_ = r_;
01096       }
01097    
01098       // rTz_old = r'*z
01099       MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
01100 
01101       if ( existU_ ) {
01102         // mu = UTAU\(AU'*z);
01103         Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
01104         index.resize(recycleBlocks_);
01105         for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01106         Teuchos::RCP<MV> AUtmp  = MVT::CloneView( *AU_, index );
01107         MVT::MvTransMv( one, *AUtmp, *z_, mu );
01108         Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
01109         char TRANS = 'N';
01110         int info;
01111         lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
01112         TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
01113                            "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
01114         // p  = z - U*mu;
01115         index.resize(recycleBlocks_);
01116         for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01117         Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_, index );
01118         index.resize( 1 );
01119         index[0] = 0;
01120         Teuchos::RCP<MV> Ptmp  = MVT::CloneView( *P_, index );
01121         MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
01122         MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
01123       } else {
01124         // p = z;
01125         index.resize( 1 );
01126         index[0] = 0;
01127         Teuchos::RCP<MV> Ptmp  = MVT::CloneView( *P_, index );
01128         MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
01129       }
01130 
01131       // Set the new state and initialize the solver.
01132       RCGIterState<ScalarType,MV> newstate;
01133 
01134       // Create RCP views here 
01135       index.resize( numBlocks_+1 );
01136       for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
01137       newstate.P  = MVT::CloneView( *P_,  index );
01138       index.resize( recycleBlocks_ );
01139       for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01140       newstate.U  = MVT::CloneView( *U_,  index );
01141       index.resize( recycleBlocks_ );
01142       for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01143       newstate.AU  = MVT::CloneView( *AU_,  index );
01144       newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) ); 
01145       newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) ); 
01146       newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) ); 
01147       newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 
01148       newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) ); 
01149       // assign the rest of the values in the struct
01150       newstate.curDim = 1; // We have initialized the first search vector
01151       newstate.Ap = Ap_;
01152       newstate.r = r_;
01153       newstate.z = z_;
01154       newstate.existU = existU_;
01155       newstate.ipiv = ipiv_;
01156       newstate.rTz_old = rTz_old_;
01157 
01158       rcg_iter->initialize(newstate);
01159 
01160       while(1) {
01161 
01162   // tell rcg_iter to iterate
01163   try {
01164     rcg_iter->iterate();
01165 
01167     //
01168     // check convergence first
01169     //
01171     if ( convTest_->getStatus() == Passed ) {
01172       // We have convergence
01173             break; // break from while(1){rcg_iter->iterate()}
01174     }
01176     //
01177     // check for maximum iterations
01178     //
01180     else if ( maxIterTest_->getStatus() == Passed ) {
01181       // we don't have convergence
01182       isConverged = false;
01183             break; // break from while(1){rcg_iter->iterate()}
01184     }
01186     //
01187     // check if cycle complete; update for next cycle
01188     //
01190     else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
01191             // index into P_ of last search vector generated this cycle
01192             int lastp = -1;
01193             // index into Beta_ of last entry generated this cycle
01194             int lastBeta = -1;
01195             if (recycleBlocks_ > 0) {
01196               if (!existU_) {
01197                 if (cycle == 0) { // No U, no U1
01198 
01199                    Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
01200                    Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
01201                    Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
01202                    Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
01203                    Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
01204                    Ftmp.putScalar(zero);
01205                    Gtmp.putScalar(zero);
01206                    for (int ii=0;ii<numBlocks_;ii++) {
01207                      Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
01208                      if (ii > 0) {
01209                        Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01210                        Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01211                      }
01212                      Ftmp(ii,ii) = Dtmp(ii,0); 
01213                    }
01214 
01215                    // compute harmonic Ritz vectors
01216                    Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
01217                    getHarmonicVecs(Ftmp,Gtmp,Ytmp);
01218 
01219                    // U1 = [P(:,1:end-1)*Y];
01220                    index.resize( numBlocks_ );
01221                    for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
01222                    Teuchos::RCP<MV> Ptmp = MVT::CloneView( *P_,  index );
01223                    index.resize( recycleBlocks_ );
01224                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01225                    Teuchos::RCP<MV> U1tmp  = MVT::CloneView( *U1_,  index );
01226                    MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
01227 
01228                    // Precompute some variables for next cycle
01229 
01230                    // AU1TAU1     = Y'*G*Y;
01231                    Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
01232                    GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
01233                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
01234                    AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
01235 
01236                    // AU1TU1      = Y'*F*Y;
01237                    Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
01238                    FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
01239                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
01240                    AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
01241 
01242                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
01243                    // Must reinitialize AU1TAP; can become dense later
01244                    AU1TAPtmp.putScalar(zero);
01245                    // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
01246                    ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
01247                    for (int ii=0; ii<recycleBlocks_; ++ii) {
01248                       AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
01249                    }
01250  
01251                    // indicate that updated recycle space now defined
01252                    existU1_ = true;
01253 
01254                    // Indicate the size of the P, Beta structures generated this cycle
01255                    lastp = numBlocks_;
01256                    lastBeta = numBlocks_-1;
01257 
01258                 } // if (cycle == 0)
01259                 else { // No U, but U1 guaranteed to exist now
01260 
01261                    // Finish computation of subblocks
01262                    // AU1TAP = AU1TAP * D(1);
01263                    Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
01264                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
01265                    AU1TAPtmp.scale(Dtmp(0,0));
01266 
01267                    Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
01268                    Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
01269                    Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
01270                    APTAPtmp.putScalar(zero);
01271                    for (int ii=0; ii<numBlocks_; ii++) {
01272                      APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
01273                      if (ii > 0) {
01274                        APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01275                        APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01276                      }
01277                    }
01278 
01279                    // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
01280                    Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01281                    Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
01282                    Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01283                    Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01284                    Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01285                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
01286                    F11.assign(AU1TU1tmp);
01287                    F12.putScalar(zero);
01288                    F21.putScalar(zero);
01289                    F22.putScalar(zero);
01290                    for(int ii=0;ii<numBlocks_;ii++) {
01291                      F22(ii,ii) = Dtmp(ii,0);
01292                    }
01293 
01294                    // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
01295                    Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01296                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
01297                    Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
01298                    Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01299                    Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01300                    Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01301                    G11.assign(AU1TAU1tmp);
01302                    G12.assign(AU1TAPtmp);
01303                    // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
01304                    for (int ii=0;ii<recycleBlocks_;++ii)
01305                      for (int jj=0;jj<numBlocks_;++jj)
01306                        G21(jj,ii) = G12(ii,jj);
01307                    G22.assign(APTAPtmp);
01308 
01309                    // compute harmonic Ritz vectors
01310                    Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
01311                    getHarmonicVecs(Ftmp,Gtmp,Ytmp);
01312 
01313                    // U1 = [U1 P(:,2:end-1)]*Y;
01314                    index.resize( numBlocks_ );
01315                    for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
01316                    Teuchos::RCP<MV> Ptmp = MVT::CloneView( *P_,  index );
01317                    index.resize( recycleBlocks_ );
01318                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01319                    Teuchos::RCP<MV> PY2tmp = MVT::CloneView( *PY2_,  index );
01320                    Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01321                    index.resize( recycleBlocks_ );
01322                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01323                    Teuchos::RCP<MV> U1tmp  = MVT::CloneView( *U1_,  index );
01324                    index.resize( recycleBlocks_ );
01325                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01326                    Teuchos::RCP<MV> U1Y1tmp  = MVT::CloneView( *U1Y1_,  index );
01327                    Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
01328                    MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
01329                    MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
01330                    MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
01331 
01332                    // Precompute some variables for next cycle
01333  
01334                    // AU1TAU1     = Y'*G*Y;
01335                    Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01336                    GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
01337                    AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
01338 
01339                    // AU1TAP      = zeros(k,m);
01340                    // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
01341                    AU1TAPtmp.putScalar(zero);
01342                    ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
01343                    for (int ii=0; ii<recycleBlocks_; ++ii) {
01344                       AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
01345                    }
01346 
01347                    // AU1TU1      = Y'*F*Y;
01348                    Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01349                    FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
01350                    //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
01351                    AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
01352  
01353                    // Indicate the size of the P, Beta structures generated this cycle
01354                    lastp = numBlocks_+1;
01355                    lastBeta = numBlocks_;
01356 
01357                 } // if (cycle != 1)
01358               } // if (!existU_)
01359               else { // U exists
01360                 if (cycle == 0) { // No U1, but U exists
01361                    Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
01362                    Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
01363                    Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
01364                    Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
01365                    APTAPtmp.putScalar(zero);
01366                    for (int ii=0; ii<numBlocks_; ii++) {
01367                      APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
01368                      if (ii > 0) {
01369                        APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01370                        APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01371                      }
01372                    }
01373 
01374                    Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
01375                    L2tmp.putScalar(zero);
01376                    for(int ii=0;ii<numBlocks_;ii++) {
01377                      L2tmp(ii,ii)   = 1./Alphatmp(ii,0);
01378                      L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
01379                    }
01380 
01381                    // AUTAP = UTAU*Delta*L2;
01382                    Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
01383                    Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
01384                    Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
01385                    Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
01386                    DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
01387                    AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
01388 
01389                    // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
01390                    Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01391                    Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
01392                    Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01393                    Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01394                    Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01395                    F11.assign(UTAUtmp);
01396                    F12.putScalar(zero);
01397                    F21.putScalar(zero);
01398                    F22.putScalar(zero);
01399                    for(int ii=0;ii<numBlocks_;ii++) {
01400                      F22(ii,ii) = Dtmp(ii,0);
01401                    }
01402  
01403                    // G = [AUTAU AUTAP; AUTAP' APTAP];
01404                    Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01405                    Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
01406                    Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01407                    Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01408                    Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01409                    Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
01410                    G11.assign(AUTAUtmp);
01411                    G12.assign(AUTAPtmp);
01412                    // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
01413                    for (int ii=0;ii<recycleBlocks_;++ii)
01414                      for (int jj=0;jj<numBlocks_;++jj)
01415                        G21(jj,ii) = G12(ii,jj);
01416                    G22.assign(APTAPtmp);
01417 
01418                    // compute harmonic Ritz vectors
01419                    Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
01420                    getHarmonicVecs(Ftmp,Gtmp,Ytmp);
01421 
01422                    // U1 = [U P(:,1:end-1)]*Y;
01423                    index.resize( recycleBlocks_ );
01424                    for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
01425                    Teuchos::RCP<MV> Utmp = MVT::CloneView( *U_,  index );
01426                    index.resize( numBlocks_ );
01427                    for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
01428                    Teuchos::RCP<MV> Ptmp = MVT::CloneView( *P_,  index );
01429                    index.resize( recycleBlocks_ );
01430                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01431                    Teuchos::RCP<MV> PY2tmp = MVT::CloneView( *PY2_,  index );
01432                    Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01433                    index.resize( recycleBlocks_ );
01434                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01435                    Teuchos::RCP<MV> UY1tmp  = MVT::CloneView( *U1Y1_,  index );
01436                    Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
01437                    index.resize( recycleBlocks_ );
01438                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01439                    Teuchos::RCP<MV> U1tmp  = MVT::CloneView( *U1_,  index );
01440                    MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
01441                    MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
01442                    MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
01443 
01444                    // Precompute some variables for next cycle
01445  
01446                    // AU1TAU1     = Y'*G*Y;
01447                    Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01448                    GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
01449                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
01450                    AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
01451  
01452                    // AU1TU1      = Y'*F*Y;
01453                    Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01454                    FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
01455                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
01456                    AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
01457 
01458                    // AU1TU   = UTAU;
01459                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
01460                    AU1TUtmp.assign(UTAUtmp);
01461 
01462                    // dold    = D(end);
01463                    dold = Dtmp(numBlocks_-1,0);
01464 
01465                    // indicate that updated recycle space now defined
01466                    existU1_ = true; 
01467 
01468                    // Indicate the size of the P, Beta structures generated this cycle
01469                    lastp = numBlocks_;
01470                    lastBeta = numBlocks_-1;
01471                 }
01472                 else { // Have U and U1
01473                    Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
01474                    Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
01475                    Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
01476                    Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
01477                    for (int ii=0; ii<numBlocks_; ii++) {
01478                      APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
01479                      if (ii > 0) {
01480                        APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01481                        APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
01482                      }
01483                    }
01484 
01485                    Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
01486                    for(int ii=0;ii<numBlocks_;ii++) {
01487                      L2tmp(ii,ii)   = 1./Alphatmp(ii,0);
01488                      L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
01489                    }
01490 
01491                    // M(end,1) = dold*(-Beta(1)/Alpha(1));
01492                    // AU1TAP = Y'*[AU1TU*Delta*L2; M];
01493                    Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
01494                    Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
01495                    DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
01496                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
01497                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
01498                    AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
01499                    Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
01500                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
01501                    AU1TAPtmp.putScalar(zero);
01502                    AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
01503                    Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01504                    ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
01505                    for(int ii=0;ii<recycleBlocks_;ii++) {
01506                      AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
01507                    }
01508 
01509                    // AU1TU = Y1'*AU1TU
01510                    Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
01511                    Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
01512                    AU1TUtmp.assign(Y1TAU1TU);
01513 
01514                    // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
01515                    Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01516                    Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
01517                    Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01518                    Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01519                    Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01520                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
01521                    F11.assign(AU1TU1tmp);
01522                    for(int ii=0;ii<numBlocks_;ii++) {
01523                      F22(ii,ii) = Dtmp(ii,0);
01524                    }
01525  
01526                    // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
01527                    Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
01528                    Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
01529                    Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
01530                    Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
01531                    Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
01532                    Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
01533                    G11.assign(AU1TAU1tmp);
01534                    G12.assign(AU1TAPtmp);
01535                    // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
01536                    for (int ii=0;ii<recycleBlocks_;++ii)
01537                      for (int jj=0;jj<numBlocks_;++jj)
01538                        G21(jj,ii) = G12(ii,jj);
01539                    G22.assign(APTAPtmp);
01540 
01541                    // compute harmonic Ritz vectors
01542                    Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
01543                    getHarmonicVecs(Ftmp,Gtmp,Ytmp);
01544 
01545                    // U1 = [U1 P(:,2:end-1)]*Y;
01546                    index.resize( numBlocks_ );
01547                    for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
01548                    Teuchos::RCP<MV> Ptmp = MVT::CloneView( *P_,  index );
01549                    index.resize( recycleBlocks_ );
01550                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01551                    Teuchos::RCP<MV> PY2tmp = MVT::CloneView( *PY2_,  index );
01552                    index.resize( recycleBlocks_ );
01553                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01554                    Teuchos::RCP<MV> U1tmp  = MVT::CloneView( *U1_,  index );
01555                    index.resize( recycleBlocks_ );
01556                    for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01557                    Teuchos::RCP<MV> U1Y1tmp  = MVT::CloneView( *U1Y1_,  index );
01558                    MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
01559                    MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
01560                    MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
01561  
01562                    // Precompute some variables for next cycle
01563 
01564                    // AU1TAU1     = Y'*G*Y;
01565                    Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01566                    GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
01567                    AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
01568  
01569                    // AU1TU1      = Y'*F*Y;
01570                    Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
01571                    FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
01572                    AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
01573 
01574                    // dold    = D(end);
01575                    dold = Dtmp(numBlocks_-1,0);
01576 
01577                    // Indicate the size of the P, Beta structures generated this cycle
01578                    lastp = numBlocks_+1;
01579                    lastBeta = numBlocks_;
01580 
01581                 }
01582               }
01583             } // if (recycleBlocks_ > 0)
01584 
01585             // Cleanup after end of cycle
01586 
01587             // P = P(:,end-1:end);
01588             index.resize( 2 );
01589             for (int ii=0; ii<2; ++ii) { index[ii] = ii; }
01590             Teuchos::RCP<MV> Ptmp1 = MVT::CloneView( *P_,  index );
01591             index.resize( 2 );
01592             index[0] = lastp-1; index[1] = lastp;
01593             Teuchos::RCP<MV> Ptmp2 = MVT::CloneView( *P_,  index );
01594             index.resize( 2 );
01595             for (int ii=0; ii<2; ++ii) { index[ii] = ii; }
01596             MVT::SetBlock(*Ptmp2,index,*Ptmp1);
01597 
01598             // Beta = Beta(end);
01599             (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
01600  
01601             // Delta = Delta(:,end);
01602             if (existU_) { // Delta only initialized if U exists
01603               Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
01604               Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
01605               mu1.assign(mu2);
01606             }
01607 
01608             // Now reinitialize state variables for next cycle
01609             newstate.P = Teuchos::null;
01610             index.resize( numBlocks_+1 );
01611             for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
01612             newstate.P  = MVT::CloneView( *P_,  index );
01613 
01614             newstate.Beta = Teuchos::null;
01615             newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
01616 
01617             newstate.Delta = Teuchos::null;
01618             newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
01619 
01620             newstate.curDim = 1; // We have initialized the first search vector
01621 
01622             // Pass to iteration object
01623             rcg_iter->initialize(newstate);
01624 
01625             // increment cycle count
01626             cycle = cycle + 1;
01627 
01628     }
01630     //
01631     // we returned from iterate(), but none of our status tests Passed.
01632     // something is wrong, and it is probably our fault.
01633     //
01635     else {
01636       TEST_FOR_EXCEPTION(true,std::logic_error,
01637              "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
01638     }
01639   }
01640   catch (const std::exception &e) {
01641     printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration " 
01642            << rcg_iter->getNumIters() << std::endl 
01643            << e.what() << std::endl;
01644     throw;
01645   }
01646       }
01647 
01648       // Inform the linear problem that we are finished with this block linear system.
01649       problem_->setCurrLS();
01650       
01651       // Update indices for the linear systems to be solved.
01652       numRHS2Solve -= 1;
01653       if ( numRHS2Solve > 0 ) {
01654         currIdx[0]++;
01655   // Set the next indices.
01656   problem_->setLSIndex( currIdx );
01657       } 
01658       else {
01659         currIdx.resize( numRHS2Solve );
01660       }
01661 
01662       // Update the recycle space for the next linear system
01663       if (existU1_) { // be sure updated recycle space was created
01664         // U = U1
01665         index.resize(recycleBlocks_);
01666         for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01667         MVT::SetBlock(*U1_,index,*U_);
01668         // Set flag indicating recycle space is now defined
01669         existU_ = true;
01670         if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
01671           // Free pointers in newstate
01672           newstate.P = Teuchos::null;
01673           newstate.Ap = Teuchos::null;
01674           newstate.r = Teuchos::null;
01675           newstate.z = Teuchos::null;
01676           newstate.U = Teuchos::null;
01677           newstate.AU = Teuchos::null;
01678           newstate.Alpha = Teuchos::null;
01679           newstate.Beta = Teuchos::null;
01680           newstate.D = Teuchos::null;
01681           newstate.Delta = Teuchos::null;
01682           newstate.LUUTAU = Teuchos::null;
01683           newstate.ipiv = Teuchos::null;
01684           newstate.rTz_old = Teuchos::null;
01685 
01686           // Reinitialize AU, UTAU, AUTAU
01687           index.resize(recycleBlocks_);
01688           for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01689           Teuchos::RCP<MV> Utmp  = MVT::CloneView( *U_,  index );
01690           index.resize(recycleBlocks_);
01691           for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01692           Teuchos::RCP<MV> AUtmp = MVT::CloneView( *AU_, index );
01693           // Initialize AU
01694           problem_->applyOp( *Utmp, *AUtmp );
01695           // Initialize UTAU
01696           Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
01697           MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
01698           // Initialize AUTAU  ( AUTAU = AU'*(M\AU) )
01699           Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
01700           if ( problem_->getLeftPrec() != Teuchos::null ) {
01701             index.resize(recycleBlocks_);
01702             for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
01703             index.resize(recycleBlocks_);
01704             for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
01705             Teuchos::RCP<MV> LeftPCAU = MVT::CloneView( *U1_, index ); // use U1 as temp storage
01706             problem_->applyLeftPrec( *AUtmp, *LeftPCAU );
01707             MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
01708           } else {
01709             MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
01710           }
01711         } // if (numRHS2Solve > 0)
01712 
01713       } // if (existU1) 
01714     }// while ( numRHS2Solve > 0 )
01715     
01716   }
01717 
01718   // print final summary
01719   sTest_->print( printer_->stream(FinalSummary) );
01720  
01721   // print timing information
01722   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01723 
01724   // get iteration information for this solve
01725   numIters_ = maxIterTest_->getNumIters();
01726   
01727   if (!isConverged) {
01728     return Unconverged; // return from RCGSolMgr::solve() 
01729   }
01730   return Converged; // return from RCGSolMgr::solve() 
01731 }
01732 
01733 //  Compute the harmonic eigenpairs of the projected, dense system.
01734 template<class ScalarType, class MV, class OP>
01735 void RCGSolMgr<ScalarType,MV,OP>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F,
01736                                                   const Teuchos::SerialDenseMatrix<int,ScalarType>& G,
01737                                                         Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
01738   // order of F,G 
01739   int n = F.numCols(); 
01740 
01741   // The LAPACK interface
01742   Teuchos::LAPACK<int,ScalarType> lapack;
01743  
01744   // Magnitude of harmonic Ritz values
01745   std::vector<MagnitudeType> w(n);
01746  
01747   // Sorted order of harmonic Ritz values
01748   std::vector<int> iperm(n);
01749  
01750   // Compute k smallest harmonic Ritz pairs
01751   // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
01752   int itype = 1; // solve A*x = (lambda)*B*x
01753   char jobz='V'; // compute eigenvalues and eigenvectors
01754   char uplo='U'; // since F,G symmetric, reference only their upper triangular data
01755   std::vector<ScalarType> work(1);
01756   int lwork = -1;
01757   int info = 0;
01758   // since SYGV destroys workspace, create copies of F,G
01759   Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_, F_->numRows(), F_->numCols() );
01760   Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_, G_->numRows(), G_->numCols() );
01761 
01762   // query for optimal workspace size
01763   lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 
01764   TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
01765                      "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
01766   lwork = (int)work[0];
01767   work.resize(lwork);
01768   lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 
01769   TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
01770                      "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
01771 
01772 
01773   // Construct magnitude of each harmonic Ritz value
01774   this->sort(w,n,iperm);
01775  
01776   // Select recycledBlocks_ smallest eigenvectors
01777   for( int i=0; i<recycleBlocks_; i++ ) {
01778     for( int j=0; j<n; j++ ) {
01779       Y(j,i) = G2(j,iperm[i]);
01780     }
01781   }
01782  
01783 }
01784 
01785 // This method sorts list of n floating-point numbers and return permutation vector
01786 template<class ScalarType, class MV, class OP>
01787 void RCGSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
01788 {
01789   int l, r, j, i, flag;
01790   int    RR2;
01791   double dRR, dK;
01792  
01793   // Initialize the permutation vector.
01794   for(j=0;j<n;j++)
01795     iperm[j] = j;
01796  
01797   if (n <= 1) return;
01798  
01799   l    = n / 2 + 1;
01800   r    = n - 1;
01801   l    = l - 1;
01802   dRR  = dlist[l - 1];
01803   dK   = dlist[l - 1];
01804  
01805   RR2 = iperm[l - 1];
01806   while (r != 0) {
01807     j = l;
01808     flag = 1;
01809  
01810     while (flag == 1) {
01811       i = j;
01812       j = j + j;
01813  
01814       if (j > r + 1)
01815         flag = 0;
01816       else {
01817         if (j < r + 1)
01818           if (dlist[j] > dlist[j - 1]) j = j + 1;
01819  
01820         if (dlist[j - 1] > dK) {
01821           dlist[i - 1] = dlist[j - 1];
01822           iperm[i - 1] = iperm[j - 1];
01823         }
01824         else {
01825           flag = 0;
01826         }
01827       }
01828     }
01829     dlist[i - 1] = dRR;
01830     iperm[i - 1] = RR2;
01831     if (l == 1) {
01832       dRR  = dlist [r];
01833       RR2 = iperm[r];
01834       dK = dlist[r];
01835       dlist[r] = dlist[0];
01836       iperm[r] = iperm[0];
01837       r = r - 1;
01838     }
01839     else {
01840       l   = l - 1;
01841       dRR  = dlist[l - 1];
01842       RR2  = iperm[l - 1];
01843       dK   = dlist[l - 1];
01844     }
01845   }
01846   dlist[0] = dRR;
01847   iperm[0] = RR2;
01848 }
01849 
01850 //  This method requires the solver manager to return a std::string that describes itself.
01851 template<class ScalarType, class MV, class OP>
01852 std::string RCGSolMgr<ScalarType,MV,OP>::description() const
01853 {
01854   std::ostringstream oss;
01855   oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01856   return oss.str();
01857 }
01858   
01859 } // end Belos namespace
01860 
01861 #endif /* BELOS_RCG_SOLMGR_HPP */

Generated on Wed May 12 21:30:09 2010 for Belos by  doxygen 1.4.7