MueLu_NullspacePresmoothFactory_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef MUELU_NULLSPACEPRESMOOTHFACTORY_DEF_HPP
#define MUELU_NULLSPACEPRESMOOTHFACTORY_DEF_HPP

#include <Xpetra_Matrix.hpp>
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>

#include "MueLu_NullspacePresmoothFactory_decl.hpp"
#include "MueLu_Utilities.hpp"

#include "MueLu_Monitor.hpp"

namespace MueLu {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> NullspacePresmoothFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());

  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory for A");
  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory for the nonsmoothed nullspace");

  return validParamList;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void NullspacePresmoothFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &currentLevel) const {
  Input(currentLevel, "Nullspace");

  if (currentLevel.GetLevelID() == 0)
    Input(currentLevel, "A");
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void NullspacePresmoothFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
  FactoryMonitor m(*this, "Nullspace presmooth factory", currentLevel);

  RCP<MultiVector> newB;
  if (currentLevel.GetLevelID() == 0) {
    RCP<Matrix> A      = Get<RCP<Matrix> >(currentLevel, "A");
    RCP<MultiVector> B = Get<RCP<MultiVector> >(currentLevel, "Nullspace");
    newB               = MultiVectorFactory::Build(B->getMap(), B->getNumVectors());

    Teuchos::ArrayRCP<SC> D = Utilities::GetMatrixDiagonal_arcp(*A);

    SC damping = 4. / 3;
    damping /= Utilities::PowerMethod(*A, true, (LO)10, 1e-4);

    A->apply(*B, *newB, Teuchos::NO_TRANS);

    size_t numVec  = newB->getNumVectors();
    LO numElements = newB->getLocalLength();
    for (size_t j = 0; j < numVec; j++) {
      Teuchos::ArrayRCP<const SC> Bj = B->getData(j);
      Teuchos::ArrayRCP<SC> newBj    = newB->getDataNonConst(j);

      for (LO i = 0; i < numElements; i++)
        newBj[i] = Bj[i] - damping * newBj[i] / D[i];
    }
  } else {
    newB = Get<RCP<MultiVector> >(currentLevel, "Nullspace");
  }

  // provide "Nullspace" variable on current level
  Set(currentLevel, "Nullspace", newB);

}  // Build

}  // namespace MueLu

#endif  // MUELU_NULLSPACEPRESMOOTHFACTORY_DEF_HPP