src

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