Belos Package Browser (Single Doxygen Collection) Development
BelosPCPGSolMgr.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_PCPG_SOLMGR_HPP
00043 #define BELOS_PCPG_SOLMGR_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosPCPGIter.hpp"
00056 
00057 #include "BelosDGKSOrthoManager.hpp"
00058 #include "BelosICGSOrthoManager.hpp"
00059 #include "BelosIMGSOrthoManager.hpp"
00060 #include "BelosStatusTestMaxIters.hpp"
00061 #include "BelosStatusTestGenResNorm.hpp"
00062 #include "BelosStatusTestCombo.hpp"
00063 #include "BelosStatusTestOutputFactory.hpp"
00064 #include "BelosOutputManager.hpp"
00065 #include "Teuchos_BLAS.hpp"
00066 #include "Teuchos_LAPACK.hpp"
00067 #include "Teuchos_TimeMonitor.hpp"
00068 
00069 //#include <vector>              //getPermThatSorts
00070 //#include <algorithm>           //getPermThatSorts
00071 
00072 
00089 namespace Belos {
00090   
00092 
00093   
00100   class PCPGSolMgrLinearProblemFailure : public BelosError {public:
00101     PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00102     {}};
00103   
00109   class PCPGSolMgrOrthoFailure : public BelosError {public:
00110     PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00111     {}};
00112   
00119   class PCPGSolMgrLAPACKFailure : public BelosError {public:
00120     PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00121     {}};
00122 
00129   class PCPGSolMgrRecyclingFailure : public BelosError {public:
00130     PCPGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
00131     {}};
00132 
00134 
00135 
00136   template<class ScalarType, class MV, class OP>
00137   class PCPGSolMgr : public SolverManager<ScalarType,MV,OP> {
00138     
00139   private:
00140     typedef MultiVecTraits<ScalarType,MV> MVT;
00141     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00142     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00143     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00144     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00145     
00146   public:
00147     
00149 
00150    
00157      PCPGSolMgr();
00158  
00194     PCPGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00195           const Teuchos::RCP<Teuchos::ParameterList> &pl );
00196     
00198     virtual ~PCPGSolMgr() {};
00200     
00202 
00203     
00206     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00207       return *problem_;
00208     }
00209 
00212     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00213 
00216     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00217  
00223     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00224       return tuple(timerSolve_);
00225     }
00226 
00228     int getNumIters() const {
00229       return numIters_;
00230     }
00231  
00234     bool isLOADetected() const { return false; }
00235  
00237     
00239 
00240    
00242     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00243    
00245     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00246     
00248    
00250 
00251 
00255     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00257  
00259 
00260     
00278     ReturnType solve();
00279     
00281     
00284     
00286     std::string description() const;
00287     
00289     
00290   private:
00291 
00292     // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary.  Given 
00293     // the seed space U and C = A U, find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU.
00294     int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
00295 
00296     // Linear problem.
00297     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00298     
00299     // Output manager.
00300     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00301     Teuchos::RCP<std::ostream> outputStream_;
00302 
00303     // Status test.
00304     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00305     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00306     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00307     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00308 
00309     // Orthogonalization manager.
00310     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00311     
00312     // Current parameter list.
00313     Teuchos::RCP<ParameterList> params_;
00314 
00315     // Default solver values.
00316     static const MagnitudeType convtol_default_;
00317     static const MagnitudeType orthoKappa_default_;
00318     static const int maxIters_default_;
00319     static const int deflatedBlocks_default_;
00320     static const int savedBlocks_default_;
00321     static const int verbosity_default_;
00322     static const int outputStyle_default_;
00323     static const int outputFreq_default_;
00324     static const std::string label_default_;
00325     static const std::string orthoType_default_;
00326     static const Teuchos::RCP<std::ostream> outputStream_default_;
00327 
00328     // Current solver values.
00329     MagnitudeType convtol_, orthoKappa_;
00330     int numIters_, maxIters_, deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
00331     std::string orthoType_; 
00332 
00333     // Recycled subspace, its image and the residual
00334     Teuchos::RCP<MV> U_, C_, R_;
00335 
00336     // Actual dimension of current recycling subspace (<= savedBlocks_ )
00337     int dimU_; 
00338 
00339     // Timers.
00340     std::string label_;
00341     Teuchos::RCP<Teuchos::Time> timerSolve_;
00342 
00343     // Internal state variables.
00344     bool isSet_;
00345   };
00346 
00347 
00348 // Default solver values.
00349 template<class ScalarType, class MV, class OP>
00350 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00351 
00352 template<class ScalarType, class MV, class OP>
00353 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00354 
00355 template<class ScalarType, class MV, class OP>
00356 const int PCPGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00357 
00358 template<class ScalarType, class MV, class OP>
00359 const int PCPGSolMgr<ScalarType,MV,OP>::deflatedBlocks_default_ = 2;
00360 
00361 template<class ScalarType, class MV, class OP>
00362 const int PCPGSolMgr<ScalarType,MV,OP>::savedBlocks_default_ = 16;
00363 
00364 template<class ScalarType, class MV, class OP>
00365 const int PCPGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00366 
00367 template<class ScalarType, class MV, class OP>
00368 const int PCPGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00369 
00370 template<class ScalarType, class MV, class OP>
00371 const int PCPGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00372 
00373 template<class ScalarType, class MV, class OP>
00374 const std::string PCPGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00375 
00376 template<class ScalarType, class MV, class OP>
00377 const std::string PCPGSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00378 
00379 template<class ScalarType, class MV, class OP>
00380 const Teuchos::RCP<std::ostream> PCPGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00381 
00382 
00383 // Empty Constructor
00384 template<class ScalarType, class MV, class OP>
00385 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr() :
00386   outputStream_(outputStream_default_),
00387   convtol_(convtol_default_),
00388   orthoKappa_(orthoKappa_default_),
00389   maxIters_(maxIters_default_),
00390   deflatedBlocks_(deflatedBlocks_default_),
00391   savedBlocks_(savedBlocks_default_),
00392   verbosity_(verbosity_default_),
00393   outputStyle_(outputStyle_default_),
00394   outputFreq_(outputFreq_default_),
00395   orthoType_(orthoType_default_),
00396   label_(label_default_),
00397   isSet_(false)
00398 {}
00399 
00400 
00401 // Basic Constructor
00402 template<class ScalarType, class MV, class OP>
00403 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr( 
00404                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00405                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00406   problem_(problem),
00407   outputStream_(outputStream_default_),
00408   convtol_(convtol_default_),
00409   orthoKappa_(orthoKappa_default_),
00410   maxIters_(maxIters_default_),
00411   deflatedBlocks_(deflatedBlocks_default_),
00412   savedBlocks_(savedBlocks_default_),
00413   verbosity_(verbosity_default_),
00414   outputStyle_(outputStyle_default_),
00415   outputFreq_(outputFreq_default_),
00416   orthoType_(orthoType_default_),
00417   label_(label_default_),
00418   isSet_(false)
00419 {
00420   TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00421 
00422   if (!is_null(pl)) {
00423     // Set the parameters using the list that was passed in.
00424     setParameters( pl );  
00425   }
00426 }
00427 
00428 
00429 template<class ScalarType, class MV, class OP>
00430 void PCPGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00431 {
00432   // Create the internal parameter list if ones doesn't already exist.
00433   if (params_ == Teuchos::null) {
00434     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00435   }
00436   else {
00437     params->validateParameters(*getValidParameters());
00438   }
00439 
00440   // Check for maximum number of iterations
00441   if (params->isParameter("Maximum Iterations")) {
00442     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00443 
00444     // Update parameter in our list and in status test.
00445     params_->set("Maximum Iterations", maxIters_);
00446     if (maxIterTest_!=Teuchos::null)
00447       maxIterTest_->setMaxIters( maxIters_ );
00448   }
00449 
00450   // Check for the maximum numbers of saved and deflated blocks.
00451   if (params->isParameter("Num Saved Blocks")) {
00452     savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
00453     TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
00454            "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
00455 
00456     // savedBlocks > number of matrix rows and columns, not known in parameters.
00457     //TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
00458     //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
00459 
00460     // Update parameter in our list.
00461     params_->set("Num Saved Blocks", savedBlocks_);
00462   }
00463   if (params->isParameter("Num Deflated Blocks")) {
00464     deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
00465     TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
00466            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
00467 
00468     TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
00469            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
00470 
00471     // Update parameter in our list.
00472     params_->set("Num Deflated Blocks", deflatedBlocks_);
00473   }
00474 
00475   // Check to see if the timer label changed.
00476   if (params->isParameter("Timer Label")) {
00477     std::string tempLabel = params->get("Timer Label", label_default_);
00478 
00479     // Update parameter in our list and solver timer
00480     if (tempLabel != label_) {
00481       label_ = tempLabel;
00482       params_->set("Timer Label", label_);
00483       std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00484       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00485     }
00486   }
00487 
00488   // Check if the orthogonalization changed.
00489   if (params->isParameter("Orthogonalization")) {
00490     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00491     TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00492       std::invalid_argument,
00493       "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00494     if (tempOrthoType != orthoType_) {
00495       orthoType_ = tempOrthoType;
00496       // Create orthogonalization manager
00497       if (orthoType_=="DGKS") {
00498   if (orthoKappa_ <= 0) {
00499     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00500   }
00501   else {
00502     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00503     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00504   }
00505       }
00506       else if (orthoType_=="ICGS") {
00507   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00508       } 
00509       else if (orthoType_=="IMGS") {
00510   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00511       } 
00512     }  
00513   }
00514 
00515   // Check which orthogonalization constant to use.
00516   if (params->isParameter("Orthogonalization Constant")) {
00517     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00518 
00519     // Update parameter in our list.
00520     params_->set("Orthogonalization Constant",orthoKappa_);
00521     if (orthoType_=="DGKS") {
00522       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00523   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00524       }
00525     } 
00526   }
00527 
00528   // Check for a change in verbosity level
00529   if (params->isParameter("Verbosity")) {
00530     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00531       verbosity_ = params->get("Verbosity", verbosity_default_);
00532     } else {
00533       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00534     }
00535 
00536     // Update parameter in our list.
00537     params_->set("Verbosity", verbosity_);
00538     if (printer_ != Teuchos::null)
00539       printer_->setVerbosity(verbosity_);
00540   }
00541 
00542   // Check for a change in output style
00543   if (params->isParameter("Output Style")) {
00544     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00545       outputStyle_ = params->get("Output Style", outputStyle_default_);
00546     } else {
00547       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00548     }
00549 
00550     // Reconstruct the convergence test if the explicit residual test is not being used.
00551     params_->set("Output Style", outputStyle_);
00552     outputTest_ = Teuchos::null;
00553   }
00554 
00555   // output stream
00556   if (params->isParameter("Output Stream")) {
00557     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00558 
00559     // Update parameter in our list.
00560     params_->set("Output Stream", outputStream_);
00561     if (printer_ != Teuchos::null)
00562       printer_->setOStream( outputStream_ );
00563   }
00564 
00565   // frequency level
00566   if (verbosity_ & Belos::StatusTestDetails) {
00567     if (params->isParameter("Output Frequency")) {
00568       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00569     }
00570 
00571     // Update parameter in out list and output status test.
00572     params_->set("Output Frequency", outputFreq_);
00573     if (outputTest_ != Teuchos::null)
00574       outputTest_->setOutputFrequency( outputFreq_ );
00575   }
00576 
00577   // Create output manager if we need to.
00578   if (printer_ == Teuchos::null) {
00579     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00580   }  
00581   
00582   // Convergence
00583   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00584   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00585 
00586   // Check for convergence tolerance
00587   if (params->isParameter("Convergence Tolerance")) {
00588     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00589 
00590     // Update parameter in our list and residual tests.
00591     params_->set("Convergence Tolerance", convtol_);
00592     if (convTest_ != Teuchos::null)
00593       convTest_->setTolerance( convtol_ );
00594   }
00595 
00596   // Create status tests if we need to.
00597 
00598   // Basic test checks maximum iterations and native residual.
00599   if (maxIterTest_ == Teuchos::null)
00600     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00601 
00602   if (convTest_ == Teuchos::null)
00603     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00604 
00605   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00606   
00607   // Create the status test output class.
00608   // This class manages and formats the output from the status test.
00609   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00610   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00611 
00612   // Set the solver string for the output test
00613   std::string solverDesc = " PCPG ";
00614   outputTest_->setSolverDesc( solverDesc );
00615 
00616 
00617   // Create orthogonalization manager if we need to.
00618   if (ortho_ == Teuchos::null) {
00619     if (orthoType_=="DGKS") {
00620       if (orthoKappa_ <= 0) {
00621   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00622       }
00623       else {
00624   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00625   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00626       }
00627     }
00628     else if (orthoType_=="ICGS") {
00629       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00630     } 
00631     else if (orthoType_=="IMGS") {
00632       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00633     } 
00634     else {
00635       TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00636        "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
00637     }  
00638   }
00639 
00640   // Create the timer if we need to.
00641   if (timerSolve_ == Teuchos::null) {
00642     std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00643     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00644   }
00645 
00646   // Inform the solver manager that the current parameters were set.
00647   isSet_ = true;
00648 }
00649 
00650     
00651 template<class ScalarType, class MV, class OP>
00652 Teuchos::RCP<const Teuchos::ParameterList>
00653 PCPGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00654 {
00655   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00656   if (is_null(validPL)) {
00657     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00658     // Set all the valid parameters and their default values.
00659     pl->set("Convergence Tolerance", convtol_default_,
00660       "The relative residual tolerance that needs to be achieved by the\n"
00661       "iterative solver in order for the linear system to be declared converged.");
00662     pl->set("Maximum Iterations", maxIters_default_,
00663       "The maximum number of iterations allowed for each\n"
00664       "set of RHS solved.");
00665     pl->set("Num Deflated Blocks", deflatedBlocks_default_,
00666       "The maximum number of vectors in the seed subspace." );
00667     pl->set("Num Saved Blocks", savedBlocks_default_,
00668       "The maximum number of vectors saved from old Krylov subspaces." );
00669     pl->set("Verbosity", verbosity_default_,
00670       "What type(s) of solver information should be outputted\n"
00671       "to the output stream.");
00672     pl->set("Output Style", outputStyle_default_,
00673       "What style is used for the solver information outputted\n"
00674       "to the output stream.");
00675     pl->set("Output Frequency", outputFreq_default_,
00676       "How often convergence information should be outputted\n"
00677       "to the output stream.");  
00678     pl->set("Output Stream", outputStream_default_,
00679       "A reference-counted pointer to the output stream where all\n"
00680       "solver output is sent.");
00681     pl->set("Timer Label", label_default_,
00682       "The string to use as a prefix for the timer labels.");
00683     //  pl->set("Restart Timers", restartTimers_);
00684     pl->set("Orthogonalization", orthoType_default_,
00685       "The type of orthogonalization to use: DGKS, ICGS, IMGS");
00686     pl->set("Orthogonalization Constant",orthoKappa_default_,
00687       "The constant used by DGKS orthogonalization to determine\n"
00688       "whether another step of classical Gram-Schmidt is necessary.");
00689     validPL = pl;
00690   }
00691   return validPL;
00692 }
00693 
00694   
00695 // solve()
00696 template<class ScalarType, class MV, class OP>
00697 ReturnType PCPGSolMgr<ScalarType,MV,OP>::solve() {
00698 
00699   // Set the current parameters if are not set already.
00700   if (!isSet_) { setParameters( params_ ); }
00701 
00702   Teuchos::BLAS<int,ScalarType> blas;
00703   Teuchos::LAPACK<int,ScalarType> lapack;
00704   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00705   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00706   
00707   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PCPGSolMgrLinearProblemFailure,
00708                      "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
00709 
00710   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PCPGSolMgrLinearProblemFailure,
00711                      "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00712 
00713   // Create indices for the linear systems to be solved.
00714   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00715   std::vector<int> currIdx(1);
00716   currIdx[0] = 0;
00717 
00718    bool debug = false;
00719 
00720   // Inform the linear problem of the current linear system to solve.
00721   problem_->setLSIndex( currIdx ); // block size == 1
00722 
00723   // Assume convergence is achieved, then let any failed convergence set this to false.
00724   bool isConverged = true;  
00725 
00727   // PCPG iteration parameter list
00728   Teuchos::ParameterList plist;
00729   plist.set("Saved Blocks", savedBlocks_);
00730   plist.set("Block Size", 1);
00731   plist.set("Keep Diagonal", true);
00732   plist.set("Initialize Diagonal", true);
00733   
00735   // PCPG solver
00736   
00737   Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
00738   pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00739   // Number of iterations required to generate initial recycle space (if needed)
00740 
00741   // Enter solve() iterations
00742   {
00743     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00744     while ( numRHS2Solve > 0 ) {  // test for quick return
00745       
00746       // Reset the status test.
00747       outputTest_->reset();
00748 
00749       // Create the first block in the current Krylov basis (residual).
00750       if (R_ == Teuchos::null)
00751         R_ = MVT::Clone( *(problem_->getRHS()), 1 ); 
00752 
00753       problem_->computeCurrResVec( &*R_ );
00754 
00755       
00756       // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00757       // TODO: ensure hypothesis right here ... I have to think about use cases.
00758 
00759       if( U_ != Teuchos::null ){
00760         // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00761 
00762         // possibly over solved equation ...  I want residual norms
00763         // relative to the initial residual, not what I am about to compute.
00764         Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
00765         std::vector<MagnitudeType> rnorm0(1);
00766         MVT::MvNorm( *R_, rnorm0 ); // rnorm0  = norm(R_);
00767 
00768         // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
00769         std::cout  << "Solver Manager:  dimU_ = " << dimU_ << std::endl;
00770         Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
00771 
00772         Teuchos::RCP<const MV> Uactive, Cactive;
00773         std::vector<int> active_columns( dimU_ );
00774         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
00775         Uactive = MVT::CloneView(*U_, active_columns);
00776         Cactive = MVT::CloneView(*C_, active_columns);
00777 
00778         if( debug ){
00779           std::cout << " Solver Manager : check duality of seed basis " << std::endl;
00780           Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
00781           MVT::MvTransMv( one, *Uactive, *Cactive, H );
00782           H.print( std::cout );
00783         }
00784         
00785         MVT::MvTransMv( one, *Uactive, *R_, Z );
00786         Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
00787         MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );  // UZ
00788         MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += tmp;
00789         MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );  // CZ
00790         MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= tmp;
00791         std::vector<MagnitudeType> rnorm(1);
00792         MVT::MvNorm( *R_, rnorm );
00793         if( rnorm[0] < rnorm0[0] * .001 ){  //reorthogonalize
00794           MVT::MvTransMv( one, *Uactive, *R_, Z );
00795           MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
00796           MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += UZ;
00797           MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
00798           MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= CZ;
00799         }
00800         Uactive = Teuchos::null;
00801         Cactive = Teuchos::null;
00802         tempU = Teuchos::null;
00803       }
00804       else { 
00805         dimU_ = 0;
00806       }
00807 
00808 
00809       // Set the new state and initialize the solver.
00810       PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
00811 
00812       pcpgState.R = R_;
00813       if( U_ != Teuchos::null ) pcpgState.U = U_;
00814       if( C_ != Teuchos::null ) pcpgState.C = C_;
00815       if( dimU_ > 0 ) pcpgState.curDim = dimU_;
00816       pcpg_iter->initialize(pcpgState);
00817 
00818       // treat initialize() exceptions here?  how to use try-catch-throw? DMD
00819 
00820       // Get the current number of deflated blocks with the PCPG iteration
00821       dimU_ = pcpgState.curDim;
00822       if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
00823       pcpg_iter->resetNumIters();
00824 
00825       if( dimU_ > savedBlocks_ )
00826         std::cout << "Error: dimU_  = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl; 
00827 
00828       while(1) { // dummy loop for break
00829 
00830         // tell pcpg_iter to iterate
00831         try {
00832           if( debug ) printf("********** Calling iterate...\n");
00833           pcpg_iter->iterate();
00834 
00836           //
00837           // check convergence first
00838           //
00840           if ( convTest_->getStatus() == Passed ) {
00841             // we have convergence
00842             break;  // break from while(1){pcpg_iter->iterate()}
00843           }
00845           //
00846           // check for maximum iterations
00847           //
00849           else if ( maxIterTest_->getStatus() == Passed ) {
00850             // we don't have convergence
00851             isConverged = false;
00852             break;  // break from while(1){pcpg_iter->iterate()}
00853           }
00854           else {
00855 
00857           //
00858           // we returned from iterate(), but none of our status tests Passed.
00859           // Something is wrong, and it is probably the developers fault.
00860           //
00862 
00863             TEST_FOR_EXCEPTION(true,std::logic_error,
00864                                "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
00865           } // end if
00866         } // end try
00867         catch (const PCPGIterOrthoFailure &e) {
00868 
00869           // Check to see if the most recent solution yielded convergence.
00870           sTest_->checkStatus( &*pcpg_iter );
00871           if (convTest_->getStatus() != Passed)
00872             isConverged = false;
00873           break;
00874         }
00875         catch (const std::exception &e) {
00876           printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
00877                                    << pcpg_iter->getNumIters() << std::endl
00878                                    << e.what() << std::endl;
00879           throw;
00880         }
00881       } // end of while(1)
00882 
00883       // Update the linear problem.
00884       Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
00885       problem_->updateSolution( update, true );
00886 
00887       // Inform the linear problem that we are finished with this block linear system.
00888       problem_->setCurrLS();
00889 
00890       // Get the state.   How did pcpgState die?
00891       PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
00892 
00893       dimU_ = oldState.curDim;
00894       int q = oldState.prevUdim;
00895 
00896       std::cout << "SolverManager: dimU_ " << dimU_ << "   prevUdim= " << q << std::endl;
00897 
00898       if( q > deflatedBlocks_ ) 
00899         std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
00900 
00901       int rank; 
00902       if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
00903         //Given the seed space U and C = A U for some symmetric positive definite A,
00904         //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
00905 
00906         //oldState.D->print( std::cout ); D = diag( C'*U )
00907 
00908         U_ = oldState.U; //MVT::MvPrint( *U_, std::cout ); 
00909         C_ = oldState.C; //MVT::MvPrint( *C_, std::cout ); 
00910         rank = ARRQR(dimU_,q, *oldState.D );
00911         if( rank < dimU_ ) {
00912                 std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
00913         }  
00914         dimU_ = rank;
00915 
00916       } // Now U_ and C_ = AU are dual bases.
00917 
00918       if( dimU_ > deflatedBlocks_ ){
00919 
00920         if( !deflatedBlocks_ ){
00921            U_ = Teuchos::null;
00922            C_ = Teuchos::null;
00923            dimU_ = deflatedBlocks_;
00924            break;
00925         }
00926 
00927         bool Harmonic = false; // (Harmonic) Ritz vectors
00928      
00929         Teuchos::RCP<MV> Uorth;
00930 
00931         std::vector<int> active_cols( dimU_ ); 
00932         for (int i=0; i < dimU_; ++i) active_cols[i] = i;
00933 
00934         if( Harmonic ){
00935           Uorth = MVT::CloneCopy(*C_, active_cols); 
00936         }
00937         else{
00938           Uorth = MVT::CloneCopy(*U_, active_cols); 
00939         }
00940 
00941         // Explicitly construct Q and R factors 
00942         Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
00943         rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
00944         Uorth = Teuchos::null;
00945         // TODO:  During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.  
00946         // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
00947 
00948         // throw an error if U is both A-orthonormal and rank deficient
00949         TEST_FOR_EXCEPTION(rank < dimU_,PCPGSolMgrOrthoFailure,
00950                            "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
00951 
00952 
00953         // R VT' = Ur S,   
00954         Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
00955         Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
00956         int lwork = 5*dimU_;                           // minimal, extra computation < 67*dimU_
00957         int info = 0;  // Hermite
00958         int lrwork = 1;
00959         if( problem_->isHermitian() ) lrwork = dimU_;
00960         std::vector<ScalarType> work(lwork); // 
00961         std::vector<ScalarType> Svec(dimU_); // 
00962         std::vector<ScalarType> rwork(lrwork); 
00963         lapack.GESVD('N', 'O',
00964                    R.numRows(),R.numCols(),R.values(), R.numRows(),
00965                    &Svec[0],
00966                    Ur.values(),1,
00967                    VT.values(),1, // Output: VT stored in R
00968                    &work[0], lwork,
00969                    &rwork[0], &info);
00970 
00971         TEST_FOR_EXCEPTION(info != 0, PCPGSolMgrLAPACKFailure,
00972            "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
00973 
00974         if( work[0] !=  67. * dimU_ )
00975            std::cout << " SVD " << dimU_ <<  " lwork " << work[0]  << std::endl;
00976         for( int i=0; i< dimU_; i++)
00977            std::cout << i << " " << Svec[i] << std::endl;
00978            
00979         Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
00980 
00981         int startRow = 0, startCol = 0;
00982         if( Harmonic )
00983           startCol = dimU_ - deflatedBlocks_;
00984 
00985         Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
00986                                                      wholeV,
00987                                                      wholeV.numRows(),
00988                                                      deflatedBlocks_,
00989                                          startRow,
00990                                          startCol);
00991         std::vector<int> active_columns( dimU_ );
00992         std::vector<int> def_cols( deflatedBlocks_ );
00993         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
00994         for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
00995 
00996         Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
00997         Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
00998         MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); //  U:= U*V
00999         Ucopy   = Teuchos::null;
01000         Uactive = Teuchos::null;
01001         Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
01002         Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
01003         MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); //  C:= C*V
01004         Ccopy  = Teuchos::null;
01005         Cactive = Teuchos::null;
01006         dimU_ = deflatedBlocks_;
01007       }
01008       printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
01009 
01010       // Inform the linear problem that we are finished with this block linear system.
01011       problem_->setCurrLS();
01012 
01013       // Update indices for the linear systems to be solved.
01014       numRHS2Solve -= 1;
01015       if ( numRHS2Solve > 0 ) {
01016   currIdx[0]++;
01017 
01018         // Set the next indices.
01019         problem_->setLSIndex( currIdx );
01020       }
01021       else {
01022         currIdx.resize( numRHS2Solve );
01023       }
01024     }// while ( numRHS2Solve > 0 )
01025   }
01026     
01027   // print final summary
01028   sTest_->print( printer_->stream(FinalSummary) );
01029   
01030   // print timing information
01031   Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01032 
01033   // get iteration information for this solve
01034   numIters_ = maxIterTest_->getNumIters();
01035  
01036   if (!isConverged) {
01037     return Unconverged; // return from PCPGSolMgr::solve() 
01038   }
01039   return Converged; // return from PCPGSolMgr::solve() 
01040 }
01041 
01042 // A-orthogonalize the Seed Space
01043 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
01044 // that are not rank revealing, and are not designed for PCPG in other ways too.
01045 template<class ScalarType, class MV, class OP>
01046 int PCPGSolMgr<ScalarType,MV,OP>::ARRQR(int p, int q, const Teuchos::SerialDenseMatrix<int,ScalarType>& D)
01047 {
01048   using Teuchos::RCP;
01049   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01050   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01051 
01052   // Allocate memory for scalars.
01053   Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
01054   Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
01055   Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
01056   std::vector<int> curind(1);
01057   std::vector<int> ipiv(p - q); // RRQR Pivot indices
01058   std::vector<ScalarType> Pivots(p); // RRQR Pivots
01059   int i, imax, j, k, l;
01060   ScalarType rteps = 1.5e-8;
01061 
01062   // Scale such that diag( U'C) = I
01063   for( int i = q ; i < p ; i++ ){
01064     ipiv[i-q] = i;
01065     curind[0] = i;
01066     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01067     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01068     anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
01069     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
01070     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
01071     Pivots[i]  = one;
01072   }
01073 
01074   for( i = q ; i < p ; i++ ){  
01075     if( q < i && i < p-1 ){ // Find the largest pivot
01076       imax = i;
01077       l = ipiv[imax-q];
01078       for( j = i+1 ; j < p ; j++ ){  
01079          k = ipiv[j-q];  
01080          if( Pivots[k] > Pivots[l] ){
01081            imax = j;  
01082            l = k;
01083          }
01084       } // end for
01085       if( imax > i ){ 
01086           l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
01087           ipiv[imax-q] = ipiv[i-q];
01088           ipiv[i-q] = l;
01089       }
01090     } // largest pivot found
01091     int k = ipiv[i-q];
01092  
01093     if( Pivots[k]  > 1.5625e-2 ){
01094       anorm(0,0) =  Pivots[k]; // A-norm of u
01095     }
01096     else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
01097       curind[0] = k;
01098       RCP<const MV> P = MVT::CloneView(*U_,curind);
01099       RCP<const MV> AP = MVT::CloneView(*C_,curind);
01100       MVT::MvTransMv( one, *P, *AP, anorm );
01101       anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
01102     }
01103     if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
01104        /*
01105        C(:,k) = A*U(:,k);  % Change C
01106        fixC = U(:, ipiv(1:i-1) )'*C(:,k);
01107        U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
01108        C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
01109        anorm = sqrt( U(:,k)'*C(:,k) );
01110        */
01111        std::cout << "ARRQR: Bad case not implemented" << std::endl;
01112     }
01113     if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
01114        std::cout << "ARRQR : deficient case not implemented " << std::endl;
01115        //U = U(:, ipiv(1:i-1) );
01116        //C = C(:, ipiv(1:i-1) );
01117        p = q + i; 
01118        // update curDim_ in State
01119        break; 
01120     }
01121     curind[0] = k;
01122     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01123     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01124     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
01125     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
01126     P = Teuchos::null;
01127     AP = Teuchos::null;
01128     Pivots[k] = one;                 // delete,  for diagonostic purposes
01129     P = MVT::CloneViewNonConst(*U_,curind);  // U(:,k)
01130     AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
01131     for( j = i+1 ; j < p ; j++ ){
01132       l = ipiv[j-q];   // ahhh
01133       curind[0] = l;
01134       RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault,  j=i+1=5
01135       MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
01136       MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
01137       Q = Teuchos::null;
01138       RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
01139       MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
01140       AQ = Teuchos::null;
01141       gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
01142       if( gamma(0,0) > 0){
01143         Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
01144       }
01145       else {
01146         Pivots[l] = zero; //rank deficiency revealed
01147       }
01148     }
01149   }
01150   return p;
01151 }
01152 
01153 //  The method returns a string describing the solver manager.
01154 template<class ScalarType, class MV, class OP>
01155 std::string PCPGSolMgr<ScalarType,MV,OP>::description() const
01156 {
01157   std::ostringstream oss;
01158   oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01159   oss << "{";
01160   oss << "Ortho Type='"<<orthoType_;
01161   oss << "}";
01162   return oss.str();
01163 }
01164   
01165 } // end Belos namespace
01166 
01167 #endif /* BELOS_PCPG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines