Constraints.hpp

#pragma once

#include <algorithm>
#include <vector>

#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>

#include <amdis/linearalgebra/Constraints.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp>

namespace AMDiS
{
  template <class T>
  struct Constraints<MTLSparseMatrix<T>>
  {
    using Matrix = MTLSparseMatrix<T>;
    using Vector = MTLVector<T>;

    template <class BitVector, class Associations>
    static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
    {
      SymmetryStructure const symmetry = mat.symmetry();
      if (symmetry == SymmetryStructure::spd ||
          symmetry == SymmetryStructure::symmetric ||
          symmetry == SymmetryStructure::hermitian)
        symmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
      else
        unsymmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
    }


    template <class Mat, class Vec, class BitVector, class Associations>
    static void symmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
    {
      error_exit("Not implemented");
    }


    template <class Value>
    struct Triplet
    {
      std::size_t row, col;
      Value value;
    };

    template <class Mat, class Vec, class BitVector, class Associations>
    static void unsymmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
    {
      std::vector<Triplet<typename Mat::value_type>> rowValues;
      rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));

      // Define the property maps
      // auto row   = mtl::mat::row_map(mat);
      auto col   = mtl::mat::col_map(mat);
      auto value = mtl::mat::value_map(mat);

      // iterate over the matrix
      std::size_t slotSize = 0;
      for (auto r : mtl::rows_of(mat)) {
        if (left[r.value()]) {
          slotSize = std::max(slotSize, std::size_t(mat.nnz_local(r.value())));
          std::size_t right = left2right.at(r.value());

          for (auto i : mtl::nz_of(r)) {
            rowValues.push_back({right,col(i),value(i)});
            value(i, 0);
          }
        }
      }

      mtl::mat::inserter<Mat, mtl::update_plus<typename Mat::value_type> > ins(mat, 2*slotSize);
      for (auto const& t : rowValues)
        ins[t.row][t.col] += t.value;

      for (std::size_t i = 0; i < mtl::size(left); ++i) {
        if (left[i]) {
          std::size_t j = left2right.at(i);
          if (setDiagonal) {
            ins[i][i] = 1;
            ins[i][j] = -1;
          }

          rhs[j] += rhs[i];
          rhs[i] = 0;

          sol[j] = sol[i];
        }
      }
    }

  };

} // end namespace AMDiS