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