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