MueLu_Hierarchy_decl.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_HIERARCHY_DECL_HPP
#define MUELU_HIERARCHY_DECL_HPP
 
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_Ptr.hpp>
 
#include <Xpetra_ConfigDefs.hpp>  // global_size_t
#include <Xpetra_Matrix_fwd.hpp>
#include <Xpetra_MultiVector_fwd.hpp>
#include <Xpetra_MultiVectorFactory_fwd.hpp>
#include <Xpetra_Operator_fwd.hpp>
 
#include <Xpetra_Cloner.hpp>
#include <MueLu_SmootherCloner.hpp>
#include "MueLu_ConfigDefs.hpp"
#include "MueLu_BaseClass.hpp"
#include "MueLu_Hierarchy_fwd.hpp"
 
#include "MueLu_Types.hpp"
 
#include "MueLu_FactoryBase_fwd.hpp"
#include "MueLu_FactoryManager.hpp" // no fwd declaration because constructor of FactoryManager is used as a default parameter of Setup()
#include "MueLu_HierarchyUtils_fwd.hpp"
#include "MueLu_KeepType.hpp"
#include "MueLu_Level_fwd.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_NoFactory.hpp"
#include "MueLu_PerfUtils_fwd.hpp"
#include "MueLu_PFactory_fwd.hpp"
#include "MueLu_RFactory_fwd.hpp"
#include "MueLu_SmootherBase_fwd.hpp"
#include "MueLu_SmootherFactoryBase_fwd.hpp"
#include "MueLu_SmootherFactory_fwd.hpp"
#include "MueLu_TwoLevelFactoryBase_fwd.hpp"
#include "MueLu_Utilities_fwd.hpp"
 
namespace MueLu {
 
  enum ReturnType {
    Converged,
    Unconverged,
    Undefined
  };
 
  template <class Scalar        = Xpetra::Operator<>::scalar_type,
            class LocalOrdinal  = typename Xpetra::Operator<Scalar>::local_ordinal_type,
            class GlobalOrdinal = typename Xpetra::Operator<Scalar, LocalOrdinal>::global_ordinal_type,
            class Node          = typename Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
  class Hierarchy : public BaseClass {
#undef MUELU_HIERARCHY_SHORT
#include "MueLu_UseShortNames.hpp"
 
    typedef Teuchos::ScalarTraits<SC> STS;
    typedef typename STS::magnitudeType MagnitudeType;
 
    struct ConvData {
      ConvData()                              : maxIts_(1),       tol_(-STS::magnitude(STS::one())) { }
      ConvData(LO maxIts)                     : maxIts_(maxIts),  tol_(-STS::magnitude(STS::one())) { }
      ConvData(MagnitudeType tol)             : maxIts_(10000),   tol_(tol) { }
      ConvData(std::pair<LO,MagnitudeType> p) : maxIts_(p.first), tol_(p.second) { }
 
      LO            maxIts_;
      MagnitudeType tol_;
    };
 
  public:
 
 
 
    Hierarchy();
    Hierarchy(const std::string& label);
 
    Hierarchy(const RCP<Matrix> & A);
 
    Hierarchy(const RCP<Matrix> & A, const std::string& label);
 
    virtual ~Hierarchy() { }
 
 
 
 
    static CycleType             GetDefaultCycle()                                     { return MasterList::getDefault<std::string>("cycle type") == "V" ? VCYCLE : WCYCLE; }
    static bool                  GetDefaultImplicitTranspose()                         { return MasterList::getDefault<bool>("transpose: use implicit");  }
    static Xpetra::global_size_t GetDefaultMaxCoarseSize()                             { return MasterList::getDefault<int>("coarse: max size");   }
    static int                   GetDefaultMaxLevels()                                 { return MasterList::getDefault<int>("max levels");     }
    static bool                  GetDefaultPRrebalance()                               { return MasterList::getDefault<bool>("repartition: rebalance P and R");   }
 
    Xpetra::global_size_t        GetMaxCoarseSize() const                              { return maxCoarseSize_; }
    bool                         GetImplicitTranspose() const                          { return implicitTranspose_; }
 
    void                         SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize) { maxCoarseSize_ = maxCoarseSize; }
    void                         SetPRrebalance(bool doPRrebalance)                    { doPRrebalance_ = doPRrebalance; }
    void                         SetImplicitTranspose(const bool& implicit)            { implicitTranspose_ = implicit; }
 
 
 
    template<class S2, class LO2, class GO2, class N2>
    friend class Hierarchy;
 
  private:
    int  LastLevelID()      const { return Levels_.size() - 1; }
    void DumpCurrentGraph() const;
 
  public:
 
    void AddLevel(const RCP<Level> & level);
 
    void AddNewLevel();
 
    RCP<Level> & GetLevel(const int levelID = 0);
 
    int    GetNumLevels() const;
    int    GetGlobalNumLevels() const;
 
    MagnitudeType GetRate() const { return rate_; }
 
    // This function is global
    double GetOperatorComplexity() const;
 
    // This function is global
    double GetSmootherComplexity() const;
 
    void CheckLevel(Level& level, int levelID);
 
 
    bool Setup(int coarseLevelID, const RCP<const FactoryManagerBase> fineLevelManager /* = Teuchos::null */, const RCP<const FactoryManagerBase> coarseLevelManager,
               const RCP<const FactoryManagerBase> nextLevelManager = Teuchos::null);
 
    void Setup(const FactoryManagerBase& manager = FactoryManager(), int startLevel = 0, int numDesiredLevels = GetDefaultMaxLevels());
 
    void SetupRe();
 
    void Clear(int startLevel = 0);
    void ExpertClear();
 
    CycleType GetCycle()                 const { return Cycle_;  }
 
    void      SetCycle(CycleType Cycle)        { Cycle_ = Cycle; }
 
    void SetProlongatorScalingFactor(double scalingFactor) { scalingFactor_ = scalingFactor; }
 
    ReturnType Iterate(const MultiVector& B, MultiVector& X, ConvData conv = ConvData(),
                       bool InitialGuessIsZero = false, LO startLevel = 0);
 
    void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="");
 
 
 
 
    void Keep(const std::string & ename, const FactoryBase* factory = NoFactory::get());
 
    void Delete(const std::string& ename, const FactoryBase* factory = NoFactory::get());
 
    void AddKeepFlag(const std::string & ename, const FactoryBase* factory = NoFactory::get(), KeepType keep = MueLu::Keep);
 
    void RemoveKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep = MueLu::All);
 
 
 
 
    std::string description() const; 
 
    void describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel = Default) const;
    void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
 
    // Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones
    void print(std::ostream& out = std::cout, const VerbLevel verbLevel = (MueLu::Parameters | MueLu::Statistics0)) const;
 
    void IsPreconditioner(const bool flag);
 
 
    void EnableGraphDumping(const std::string& filename, int levelID = 1) {
      isDumpingEnabled_ = true;
      dumpLevel_ = levelID;
      dumpFile_  = filename;
    }
 
    template<class Node2>
    Teuchos::RCP< Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node2> >
    clone (const RCP<Node2> &node2) const;
 
    void setlib(Xpetra::UnderlyingLib inlib) { lib_ = inlib; }
    Xpetra::UnderlyingLib lib() { return lib_; }
 
    // force recreation of cached description_ next time description() is called: 
    void ResetDescription() {
      description_ = "";
    }
 
    void AllocateLevelMultiVectors(int numvecs);
    void DeleteLevelMultiVectors();
 
  protected:
    const RCP<const FactoryManagerBase>& GetLevelManager(const int levelID) const {
      return levelManagers_[levelID];
    }
 
  private:
    Hierarchy(const Hierarchy &h);
 
    Array<RCP<Level> > Levels_;
 
    // We replace coordinates GIDs to make them consistent with matrix GIDs,
    // even if user does not do that.  Ideally, though, we should completely
    // remove any notion of coordinate GIDs, and deal only with LIDs, assuming
    // that they are consistent with matrix block IDs
    void ReplaceCoordinateMap(Level& level);
 
    // Minimum size of a matrix on any level. If we fall below that, we stop
    // the coarsening
    Xpetra::global_size_t maxCoarseSize_;
 
    // Potential speed up of the setup by skipping R construction, and using
    // transpose matrix-matrix product for RAP
    bool implicitTranspose_;
 
    // Potential speed up of the setup by skipping rebalancing of P and R, and
    // doing extra import during solve
    bool doPRrebalance_;
 
    // Hierarchy may be used in a standalone mode, or as a preconditioner
    bool isPreconditioner_;
 
    // V- or W-cycle
    CycleType Cycle_;
 
    // Scaling factor to be applied to coarse grid correction.
    double scalingFactor_;
 
    // Epetra/Tpetra mode
    Xpetra::UnderlyingLib lib_;
 
    // cache description to avoid recreating in each call to description() - use ResetDescription() to force recreation in Setup, SetupRe, etc.
    mutable std::string description_ = ""; // mutable so that we can lazily initialize in description(), which is declared const
 
    // If enabled, we dump the graph on a specified level into a specified file
    bool isDumpingEnabled_;
    int  dumpLevel_;
    std::string dumpFile_;
 
    MagnitudeType rate_;
 
    // Level managers used during the Setup
    Array<RCP<const FactoryManagerBase> > levelManagers_;
 
    // Caching (Multi)Vectors used in Hierarchy::Iterate()
    int sizeOfAllocatedLevelMultiVectors_;
    Array<RCP<MultiVector> > residual_, coarseRhs_, coarseX_, coarseImport_, coarseExport_, correction_;
    
 
  }; //class Hierarchy
 
  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  template<typename Node2>
  Teuchos::RCP<Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node2> >
  Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  clone (const Teuchos::RCP<Node2> &node2) const {
    typedef Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node2>           New_H_Type;
    typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2>      CloneOperator;
    typedef MueLu::SmootherBase<Scalar, LocalOrdinal, GlobalOrdinal, Node2> CloneSmoother;
 
    Teuchos::RCP<New_H_Type> new_h = Teuchos::rcp(new New_H_Type());
    new_h->Levels_.resize(this->GetNumLevels());
    new_h->maxCoarseSize_     = maxCoarseSize_;
    new_h->implicitTranspose_ = implicitTranspose_;
    new_h->isPreconditioner_  = isPreconditioner_;
    new_h->isDumpingEnabled_  = isDumpingEnabled_;
    new_h->dumpLevel_         = dumpLevel_;
    new_h->dumpFile_          = dumpFile_;
 
    RCP<SmootherBase>  Pre, Post;
    RCP<CloneSmoother> clonePre, clonePost;
    RCP<CloneOperator> cloneA, cloneR, cloneP;
    RCP<Operator>      A, R, P;
    for (int i = 0; i < GetNumLevels(); i++) {
      RCP<Level> level      = this->Levels_[i];
      RCP<Level> clonelevel = rcp(new Level());
 
      if (level->IsAvailable("A")) {
        A      = level->template Get<RCP<Operator> >("A");
        cloneA = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*A, node2);
        clonelevel->template Set<RCP<CloneOperator> >("A", cloneA);
      }
      if (level->IsAvailable("R")){
        R      = level->template Get<RCP<Operator> >("R");
        cloneR = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*R, node2);
        clonelevel->template Set<RCP<CloneOperator> >("R", cloneR);
      }
      if (level->IsAvailable("P")){
        P      = level->template Get<RCP<Operator> >("P");
        cloneP = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*P,  node2);
        clonelevel->template Set<RCP<CloneOperator> >("P", cloneP);
      }
      if (level->IsAvailable("PreSmoother")){
        Pre      = level->template Get<RCP<SmootherBase> >("PreSmoother");
        clonePre = MueLu::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2> (Pre, cloneA, node2);
        clonelevel->template Set<RCP<CloneSmoother> >("PreSmoother", clonePre);
      }
      if (level->IsAvailable("PostSmoother")){
        Post      = level->template Get<RCP<SmootherBase> >("PostSmoother");
        clonePost = MueLu::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2> (Post, cloneA, node2);
        clonelevel-> template Set<RCP<CloneSmoother> >("PostSmoother", clonePost);
      }
      new_h->Levels_[i] = clonelevel;
    }
 
    return new_h;
  }
 
} //namespace MueLu
 
#define MUELU_HIERARCHY_SHORT
#endif // MUELU_HIERARCHY_DECL_HPP