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