BlockDiagonalPreconditioner.hpp

#pragma once
 
#include <memory>
#include <string>
#include <utility>
#include <vector>
 
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>
 
#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/mtl/BlockPreconditioner.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp>
 
namespace AMDiS
{
  template <class M, class V, class Basis>
  class BlockDiagonalPreconditioner
      : public BlockPreconditioner<M, V, Basis>
  {
    using Super = BlockPreconditioner<M, V, Basis>;
    using Self = BlockDiagonalPreconditioner;
 
    using X = V;
    using Y = V;
    using PreconBase = PreconditionerInterface<M,X,Y>;
 
  public:
    struct Creator : CreatorInterfaceName<PreconBase>
    {
      using Interface = PreconBase;
 
      Creator(std::shared_ptr<Basis> basis)
        : basis_(std::move(basis))
      {}
 
      std::unique_ptr<Interface> createWithString(std::string prefix) override
      {
        return std::make_unique<Self>(std::move(prefix), basis_);
      }
 
      std::shared_ptr<Basis> basis_;
    };
 
  public:
    BlockDiagonalPreconditioner(std::string prefix, std::shared_ptr<Basis> basis)
      : Super(std::move(basis))
      , prefix_(std::move(prefix))
    {}
 
    void init(M const& A) override
    {
      assert(rows() == cols());
      Super::initBlocks(A, true);
      initSubPrecons();
    }
 
    void initSubPrecons()
    {
      // create preconditioners for the diagonal blocks
      subPrecon_.resize(rows());
      for (std::size_t i = 0; i < subPrecon_.size(); ++i) {
        std::string subPreconName = "default";
        std::string subPreconPrefix = prefix_ + "->sub precon[" + std::to_string(i) + "]";
        Parameters::get(subPreconPrefix, subPreconName);
 
        auto* subCreator = CreatorMap<PreconBase>::get(subPreconName, subPreconPrefix);
        subPrecon_[i] = named(subCreator)->createWithString(subPreconPrefix);
        subPrecon_[i]->init(Super::subMatrix(i,i));
      }
    }
 
    void finish() override
    {
      finishSubPrecons();
    }
 
    void finishSubPrecons()
    {
      for (std::size_t i = 0; i < subPrecon_.size(); ++i)
        subPrecon_[i]->finish();
    }
 
 
    void solve(X const& x, Y& y) const override
    {
      y.checked_change_resource(x);
 
      for (std::size_t i = 0; i < rows(); ++i)
      {
        X const x_i(Super::subVector(x, Super::rowIndices(i)));
        Y y_i(resource(x_i));
 
        subPrecon_[i]->solve(x_i, y_i);
        Super::subVector(y, Super::colIndices(i)) = y_i;
      }
    }
 
 
    void adjoint_solve(X const& x, Y& y) const override
    {
      y.checked_change_resource(x);
 
      for (std::size_t i = 0; i < cols(); ++i)
      {
        X const x_i(Super::subVector(x, Super::colIndices(i)));
        Y y_i(resource(x_i));
 
        subPrecon_[i]->adjoint_solve(x_i, y_i);
        Super::subVector(y, Super::rowIndices(i)) = y_i;
      }
    }
 
  protected:
    using Super::rows;
    using Super::cols;
 
  protected:
    std::string prefix_;
    std::vector<std::shared_ptr<PreconBase>> subPrecon_;
  };
 
 
  template <class T = double, class Basis>
  void registerBlockDiagonalPreconditioner(std::shared_ptr<Basis> const& basis, std::string name = "block_diag")
  {
    using M = typename MTLSparseMatrix<T>::BaseMatrix;
    using V = typename MTLVector<T>::BaseVector;
    using Creator = typename BlockDiagonalPreconditioner<M, V, Basis>::Creator;
    CreatorMap<typename Creator::Interface>::addCreator(name, new Creator(basis));
  }
 
} // end namespace AMDiS