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

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