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