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