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