MueLu_IfpackSmoother.cpp

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
#include "MueLu_ConfigDefs.hpp"
 
#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK)
#include <Ifpack.h>
#include <Ifpack_Chebyshev.h>
#include "Xpetra_MultiVectorFactory.hpp"
 
#include "MueLu_IfpackSmoother.hpp"
 
#include "MueLu_Level.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_Aggregates.hpp"
 
 
namespace MueLu {
 
  template <class Node>
  IfpackSmoother<Node>::IfpackSmoother(std::string const & type, Teuchos::ParameterList const & paramList, LO const &overlap)
    : type_(type), overlap_(overlap)
  {
    this->declareConstructionOutcome(false, "");
    SetParameterList(paramList);
  }
 
  template <class Node>
  void IfpackSmoother<Node>::SetParameterList(const Teuchos::ParameterList& paramList) {
    Factory::SetParameterList(paramList);
 
    if (SmootherPrototype::IsSetup()) {
      // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
      // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
      prec_->SetParameters(const_cast<ParameterList&>(this->GetParameterList()));
    }
  }
 
  template <class Node>
  void IfpackSmoother<Node>::SetPrecParameters(const Teuchos::ParameterList& list) const {
    ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
    paramList.setParameters(list);
 
    RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
 
    prec_->SetParameters(*precList);
 
    // We would like to have the following line here:
    //      paramList.setParameters(*precList);
    // For instance, if Ifpack sets somem parameters internally, we would like to have
    // them listed when we call this->GetParameterList()
    // But because of the way Ifpack handles the list, we cannot do that.
    // The bad scenario goes like this:
    //   * SmootherFactory calls Setup
    //   * Setup calls SetPrecParameters
    //   * We call prec_->SetParameters(*precList)
    //     This actually updates the internal parameter list  with default prec_ parameters
    //     This means that we get a parameter ("chebyshev: max eigenvalue", -1) in the list
    //   * Setup calls prec_->Compute()
    //     Here we may compute the max eigenvalue, but we get no indication of this. If we
    //     do compute it, our parameter list becomes outdated
    //   * SmootherFactory calls Apply
    //   * Apply constructs a list with a list with an entry "chebyshev: zero starting solution"
    //   * We call prec_->SetParameters(*precList)
    // The last call is the problem. At this point, we have a list with an outdated entry
    // "chebyshev: max eigenvalue", but prec_ uses this entry and replaces the computed max
    // eigenvalue with the one from the list, resulting in -1.0 eigenvalue.
    //
    // Ifpack2 does not have this problem, as it does not populate the list with new entries
  }
 
  template <class Node>
  void IfpackSmoother<Node>::DeclareInput(Level &currentLevel) const {
    this->Input(currentLevel, "A");
 
    if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
        type_ == "LINESMOOTHING_BANDED RELAXATION" ||
        type_ == "LINESMOOTHING_BANDEDRELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCK_RELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCK RELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCKRELAXATION") {
      this->Input(currentLevel, "CoarseNumZLayers");              // necessary for fallback criterion
      this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
    } // if (type_ == "LINESMOOTHING_BANDEDRELAXATION")
    else if (type_ == "AGGREGATE")
    {
      // Aggregate smoothing needs aggregates
      this->Input(currentLevel,"Aggregates");
    }
 
  }
 
  template <class Node>
  void IfpackSmoother<Node>::Setup(Level &currentLevel) {
    FactoryMonitor m(*this, "Setup Smoother", currentLevel);
    if (SmootherPrototype::IsSetup() == true)
      this->GetOStream(Warnings0) << "MueLu::IfpackSmoother::Setup(): Setup() has already been called" << std::endl;
 
    A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
 
    double lambdaMax = -1.0;
    if (type_ == "Chebyshev") {
      std::string maxEigString   = "chebyshev: max eigenvalue";
      std::string eigRatioString = "chebyshev: ratio eigenvalue";
 
      try {
        lambdaMax = Teuchos::getValue<Scalar>(this->GetParameter(maxEigString));
        this->GetOStream(Statistics1) << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
 
      } catch (Teuchos::Exceptions::InvalidParameterName&) {
        lambdaMax = A_->GetMaxEigenvalueEstimate();
 
        if (lambdaMax != -1.0) {
          this->GetOStream(Statistics1) << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
          this->SetParameter(maxEigString, ParameterEntry(lambdaMax));
        }
      }
 
      // Calculate the eigenvalue ratio
      const Scalar defaultEigRatio = 20;
 
      Scalar ratio = defaultEigRatio;
      try {
        ratio = Teuchos::getValue<Scalar>(this->GetParameter(eigRatioString));
 
      } catch (Teuchos::Exceptions::InvalidParameterName&) {
        this->SetParameter(eigRatioString, ParameterEntry(ratio));
      }
 
      if (currentLevel.GetLevelID()) {
        // Update ratio to be
        //   ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
        //
        // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
        RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >("A");
        size_t nRowsFine   = fineA->getGlobalNumRows();
        size_t nRowsCoarse = A_->getGlobalNumRows();
 
        ratio = std::max(ratio, as<Scalar>(nRowsFine)/nRowsCoarse);
 
        this->GetOStream(Statistics1) << eigRatioString << " (computed) = " << ratio << std::endl;
        this->SetParameter(eigRatioString, ParameterEntry(ratio));
      }
    } // if (type_ == "Chebyshev")
 
    if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
        type_ == "LINESMOOTHING_BANDED RELAXATION" ||
        type_ == "LINESMOOTHING_BANDEDRELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCK_RELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCK RELAXATION"  ||
        type_ == "LINESMOOTHING_BLOCKRELAXATION" ) {
      ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
 
      LO CoarseNumZLayers = currentLevel.Get<LO>("CoarseNumZLayers",Factory::GetFactory("CoarseNumZLayers").get());
      if (CoarseNumZLayers > 0) {
        Teuchos::ArrayRCP<LO> TVertLineIdSmoo = currentLevel.Get< Teuchos::ArrayRCP<LO> >("LineDetection_VertLineIds", Factory::GetFactory("LineDetection_VertLineIds").get());
 
        // determine number of local parts
        LO maxPart = 0;
        for(size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
          if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
        }
 
        size_t numLocalRows = A_->getLocalNumRows();
        TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
 
        if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
          myparamList.set("partitioner: type","user");
          myparamList.set("partitioner: map",&(TVertLineIdSmoo[0]));
          myparamList.set("partitioner: local parts",maxPart+1);
        } else {
          // we assume a constant number of DOFs per node
          size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
 
          // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
          Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
          for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
            for (size_t dof = 0; dof < numDofsPerNode; dof++)
              partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
          myparamList.set("partitioner: type","user");
          myparamList.set("partitioner: map",&(partitionerMap[0]));
          myparamList.set("partitioner: local parts",maxPart + 1);
        }
 
        if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
            type_ == "LINESMOOTHING_BANDED RELAXATION" ||
            type_ == "LINESMOOTHING_BANDEDRELAXATION")
          type_ = "block relaxation";
        else
          type_ = "block relaxation";
      } else {
        // line detection failed -> fallback to point-wise relaxation
        this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
        myparamList.remove("partitioner: type",false);
        myparamList.remove("partitioner: map", false);
        myparamList.remove("partitioner: local parts",false);
        type_ = "point relaxation stand-alone";
      }
 
    } // if (type_ == "LINESMOOTHING_BANDEDRELAXATION")
 
    if(type_ == "AGGREGATE") {
      SetupAggregate(currentLevel);
    }
 
    else {
      // If we're using a linear partitioner and haven't set the # local parts, set it to match the operator's block size
      ParameterList precList = this->GetParameterList();
      if(precList.isParameter("partitioner: type") && precList.get<std::string>("partitioner: type") == "linear" &&
         !precList.isParameter("partitioner: local parts")) {
        precList.set("partitioner: local parts", (int)A_->getLocalNumRows() / A_->GetFixedBlockSize());
      }
      
      
      RCP<Epetra_CrsMatrix> epA = Utilities::Op2NonConstEpetraCrs(A_);
      
      Ifpack factory;
      prec_ = rcp(factory.Create(type_, &(*epA), overlap_));
      TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(), Exceptions::RuntimeError, "Could not create an Ifpack preconditioner with type = \"" << type_ << "\"");
      SetPrecParameters();
      prec_->Compute();
    }
 
    SmootherPrototype::IsSetup(true);
 
    if (type_ == "Chebyshev" && lambdaMax == -1.0) {
      Teuchos::RCP<Ifpack_Chebyshev> chebyPrec = rcp_dynamic_cast<Ifpack_Chebyshev>(prec_);
      if (chebyPrec != Teuchos::null) {
        lambdaMax = chebyPrec->GetLambdaMax();
        A_->SetMaxEigenvalueEstimate(lambdaMax);
        this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack)" << " = " << lambdaMax << std::endl;
      }
      TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == -1.0, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
    }
 
    this->GetOStream(Statistics1) << description() << std::endl;
  }
 
 
  template <class Node>
  void IfpackSmoother<Node>::SetupAggregate(Level& currentLevel) {
 
    ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
 
    if (this->IsSetup() == true) {
      this->GetOStream(Warnings0) << "MueLu::Ifpack2moother::SetupAggregate(): Setup() has already been called" << std::endl;
      this->GetOStream(Warnings0) << "MueLu::IfpackSmoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
    }
 
    this->GetOStream(Statistics0) << "IfpackSmoother: Using Aggregate Smoothing"<<std::endl;
 
    RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates> >(currentLevel,"Aggregates");
    RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
    ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
    ArrayRCP<LO> dof_ids;
 
    // We need to unamalgamate, if the FixedBlockSize > 1
    if(A_->GetFixedBlockSize() > 1) {
      // NOTE: We're basically going to have to leave a deallocated pointer hanging out
      // in the paramList object (and inside the partitioner).  This never gets
      // use again after Compute() gets called, so this is OK, but I'm still leaving
      // this note here in case it bites us again later.
      LO blocksize = (LO) A_->GetFixedBlockSize();
      dof_ids.resize(aggregate_ids.size() * blocksize);
      for(LO i=0; i<(LO)aggregate_ids.size(); i++) {
        for(LO j=0; j<(LO)blocksize; j++)
          dof_ids[i*blocksize+j] = aggregate_ids[i];
      }
    }
    else {
      dof_ids = aggregate_ids;
    }
        
    paramList.set("partitioner: map", dof_ids.getRawPtr());
    paramList.set("partitioner: type", "user");
    paramList.set("partitioner: overlap", 0);
    paramList.set("partitioner: local parts", (int)aggregates->GetNumAggregates());
    // In case of Dirichlet nodes
    paramList.set("partitioner: keep singletons",true);
 
    RCP<Epetra_CrsMatrix> A = Utilities::Op2NonConstEpetraCrs(A_);
    type_ = "block relaxation stand-alone";
 
    Ifpack factory;
    prec_ = rcp(factory.Create(type_, &(*A), overlap_));
    TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(), Exceptions::RuntimeError, "Could not create an Ifpack preconditioner with type = \"" << type_ << "\"");
    SetPrecParameters();
 
    int rv = prec_->Compute();
    TEUCHOS_TEST_FOR_EXCEPTION(rv, Exceptions::RuntimeError, "Ifpack preconditioner with type = \"" << type_ << "\" Compute() call failed.");
 
  }
 
 
  template <class Node>
  void IfpackSmoother<Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
    TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Apply(): Setup() has not been called");
 
 
    // Forward the InitialGuessIsZero option to Ifpack
    Teuchos::ParameterList paramList;
    bool supportInitialGuess = false;
    if (type_ == "Chebyshev") {
      paramList.set("chebyshev: zero starting solution",  InitialGuessIsZero);
      supportInitialGuess = true;
 
    } else if (type_ == "point relaxation stand-alone") {
      paramList.set("relaxation: zero starting solution", InitialGuessIsZero);
      supportInitialGuess = true;
    }
 
    SetPrecParameters(paramList);
 
    // Apply
    if (InitialGuessIsZero || supportInitialGuess) {
      Epetra_MultiVector&       epX = Utilities::MV2NonConstEpetraMV(X);
      const Epetra_MultiVector& epB = Utilities::MV2EpetraMV(B);
 
      prec_->ApplyInverse(epB, epX);
 
    } else {
      RCP<MultiVector> Residual   = Utilities::Residual(*A_, X, B);
      RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
 
      Epetra_MultiVector&       epX = Utilities::MV2NonConstEpetraMV(*Correction);
      const Epetra_MultiVector& epB = Utilities::MV2EpetraMV(*Residual);
 
      prec_->ApplyInverse(epB, epX);
 
      X.update(1.0, *Correction, 1.0);
    }
  }
 
  template <class Node>
  RCP<MueLu::SmootherPrototype<double, int, int, Node> > IfpackSmoother<Node>::Copy() const {
    RCP<IfpackSmoother<Node> > smoother = rcp(new IfpackSmoother<Node>(*this) );
    smoother->SetParameterList(this->GetParameterList());
    return Teuchos::rcp_dynamic_cast<MueLu::SmootherPrototype<double, int, int, Node> >(smoother);
  }
 
  template <class Node>
  std::string IfpackSmoother<Node>::description() const {
    std::ostringstream out;
    // The check "GetVerbLevel() == Test" is to avoid
    // failures in the EasyInterface test.
    if (prec_ == Teuchos::null || this->GetVerbLevel() == InterfaceTest) {
      out << SmootherPrototype::description();
      out << "{type = " << type_ << "}";
    } else {
      out << prec_->Label();
    }
    return out.str();
  }
 
  template <class Node>
  void IfpackSmoother<Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
    MUELU_DESCRIBE;
 
    if (verbLevel & Parameters0)
      out0 << "Prec. type: " << type_ << std::endl;
 
    if (verbLevel & Parameters1) {
      out0 << "Parameter list: " << std::endl;
      Teuchos::OSTab tab2(out);
      out << this->GetParameterList();
      out0 << "Overlap: "        << overlap_ << std::endl;
    }
 
    if (verbLevel & External)
      if (prec_ != Teuchos::null) {
        Teuchos::OSTab tab2(out);
        out << *prec_ << std::endl;
      }
 
    if (verbLevel & Debug) {
      out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
           << "-" << std::endl
           << "RCP<A_>: " << A_ << std::endl
           << "RCP<prec_>: " << prec_ << std::endl;
    }
  }
 
  template <class Node>
  size_t IfpackSmoother<Node>::getNodeSmootherComplexity() const {
    // FIXME: This is a placeholder
    return Teuchos::OrdinalTraits<size_t>::invalid();
  }
  
 
} // namespace MueLu
 
// The IfpackSmoother is only templated on the Node, since it is an Epetra only object
// Therefore we do not need the full ETI instantiations as we do for the other MueLu
// objects which are instantiated on all template parameters.
#if defined(HAVE_MUELU_EPETRA)
template class MueLu::IfpackSmoother<Xpetra::EpetraNode>;
#endif
 
#endif