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