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