Unverified Commit d834e02a authored by Peter Doak's avatar Peter Doak Committed by GitHub
Browse files

Merge pull request #224 from gbalduzz/const_matrix_view 

Properly construct view of constant matrix.
parents 610c5ae0 f185e9ac

include/dca/linalg/matrix_view.hpp

+85 −93
Original line number Diff line number Diff line
@@ -25,28 +25,40 @@ namespace dca { 
namespace linalg {

// dca::linalg::



template <typename ScalarType, DeviceType device_name = linalg::CPU>

template <class Scalar, DeviceType device_name = linalg::CPU>

class MatrixView {

public:

  using ValueType = ScalarType;

  using ValueType = Scalar;

  constexpr static DeviceType device = device_name;



  using Pair = std::pair<int, int>;

  MatrixView(ScalarType* data, Pair size);

  MatrixView(ScalarType* data, Pair size, int ld);

  MatrixView(ScalarType* data, int size, int ld);

  MatrixView(ScalarType* data, int size);

  template <template <typename, DeviceType> class Matrix>

  MatrixView(Matrix<ScalarType, device_name>& mat);

  template <template <typename, DeviceType> class Matrix>

  MatrixView(Matrix<ScalarType, device_name>& mat, int offset_i, int offset_j);

  template <template <typename, DeviceType> class Matrix>

  MatrixView(Matrix<ScalarType, device_name>& mat, int offset_i, int offset_j, int ni, int nj);



  template <template <typename, DeviceType> class Matrix>

  MatrixView& operator=(const Matrix<ScalarType, device_name>& rhs);



  MatrixView& operator=(const MatrixView<ScalarType, device_name>& rhs);

  MatrixView(Scalar* data, Pair size);

  MatrixView(Scalar* data, Pair size, int ld);

  MatrixView(Scalar* data, int size, int ld);

  MatrixView(Scalar* data, int size);



  template <template <class, DeviceType> class Matrix>

  MatrixView(Matrix<Scalar, device_name>& mat);

  template <template <class, DeviceType> class Matrix>

  MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_name>& mat);

  template <template <class, DeviceType> class Matrix>

  MatrixView(Matrix<Scalar, device_name>& mat, int offset_i, int offset_j);

  template <template <class, DeviceType> class Matrix>

  MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_name>& mat, int offset_i, int offset_j);

  template <template <class, DeviceType> class Matrix>

  MatrixView(Matrix<Scalar, device_name>& mat, int offset_i, int offset_j, int ni, int nj);

  template <template <class, DeviceType> class Matrix>

  MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_name>& mat, int offset_i, int offset_j,

             int ni, int nj);



  // Preconditions: Device is CPU.

  //                Assignment from Scalar2 to Scalar is defined.

  //                Size must be equal to rhs' size.

  template <template <class, DeviceType> class Matrix, class Scalar2>

  MatrixView& operator=(const Matrix<Scalar2, device_name>& rhs);



  // Same as above, necessary to override default operator=.

  MatrixView& operator=(const MatrixView& rhs);



  void print(std::ostream& out = std::cout) const;
@@ -60,18 +72,18 @@ public: 
  __DEVICE__ __HOST__ std::pair<int, int> size() const {

    return size_;

  }

  __DEVICE__ __HOST__ ScalarType* ptr() {

  __DEVICE__ __HOST__ Scalar* ptr() {

    return ptr_;

  }

  __DEVICE__ __HOST__ const ScalarType* ptr() const {

  __DEVICE__ __HOST__ const Scalar* ptr() const {

    return ptr_;

  }

  __DEVICE__ __HOST__ ScalarType* ptr(int i, int j) {

  __DEVICE__ __HOST__ Scalar* ptr(int i, int j) {

    assert(0 <= i && i <= size_.first);

    assert(0 <= j && j <= size_.second);

    return ptr_ + leadingDimension() * j + i;

  }

  __DEVICE__ __HOST__ const ScalarType* ptr(int i, int j) const {

  __DEVICE__ __HOST__ const Scalar* ptr(int i, int j) const {

    assert(0 <= i && i <= size_.first);

    assert(0 <= j && j <= size_.second);

    return ptr_ + leadingDimension() * j + i;
@@ -85,66 +97,87 @@ public: 
  bool is_square() const {

    return size_.first == size_.second;

  }

  __DEVICE__ __HOST__ ScalarType& operator()(int i, int j) {

  __DEVICE__ __HOST__ Scalar& operator()(int i, int j) {

    assert(0 <= i && i < size_.first);

    assert(0 <= j && j < size_.second);

    return ptr_[i + j * ldm_];

  }

  __DEVICE__ __HOST__ const ScalarType& operator()(int i, int j) const {

  __DEVICE__ __HOST__ const Scalar& operator()(int i, int j) const {

    assert(0 <= i && i < size_.first);

    assert(0 <= j && j < size_.second);

    return ptr_[i + j * ldm_];

  }



private:

  ScalarType* const ptr_;

  Scalar* const ptr_;

  const int ldm_;

  const std::pair<int, int> size_;

};



template <typename ScalarType, DeviceType device_t>

MatrixView<ScalarType, device_t>::MatrixView(ScalarType* const data, const Pair size)

template <class Scalar, DeviceType device_t>

MatrixView<Scalar, device_t>::MatrixView(Scalar* const data, const Pair size)

    : MatrixView(data, size, size.first) {}



template <typename ScalarType, DeviceType device_t>

MatrixView<ScalarType, device_t>::MatrixView(ScalarType* const data, const Pair size, const int ld)

template <class Scalar, DeviceType device_t>

MatrixView<Scalar, device_t>::MatrixView(Scalar* const data, const Pair size, const int ld)

    : ptr_(data), ldm_(ld), size_(size) {}



template <typename ScalarType, DeviceType device_t>

MatrixView<ScalarType, device_t>::MatrixView(ScalarType* const data, const int size, const int ld)

template <class Scalar, DeviceType device_t>

MatrixView<Scalar, device_t>::MatrixView(Scalar* const data, const int size, const int ld)

    : ptr_(data), ldm_(ld), size_(std::make_pair(size, size)) {}



template <typename ScalarType, DeviceType device_t>

MatrixView<ScalarType, device_t>::MatrixView(ScalarType* const data, const int size)

template <class Scalar, DeviceType device_t>

MatrixView<Scalar, device_t>::MatrixView(Scalar* const data, const int size)

    : ptr_(data), ldm_(size), size_(std::make_pair(size, size)) {}



template <typename ScalarType, DeviceType device_t>

template <template <typename, DeviceType> class Matrix>

MatrixView<ScalarType, device_t>::MatrixView(Matrix<ScalarType, device_t>& mat)

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(Matrix<Scalar, device_t>& mat)

    : ptr_(mat.ptr()), ldm_(mat.leadingDimension()), size_(mat.size()) {}



template <typename ScalarType, DeviceType device_t>

template <template <typename, DeviceType> class Matrix>

MatrixView<ScalarType, device_t>::MatrixView(Matrix<ScalarType, device_t>& mat, int offset_i,

                                             int offset_j)

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_t>& mat)

    : ptr_(mat.ptr()), ldm_(mat.leadingDimension()), size_(mat.size()) {}



template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(Matrix<Scalar, device_t>& mat, int offset_i, int offset_j)

    : MatrixView(mat, offset_i, offset_j, mat.nrRows() - offset_i, mat.nrCols() - offset_j) {

  assert(offset_i < mat.nrCols());

  assert(offset_j < mat.nrRows());

}

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_t>& mat,

                                         int offset_i, int offset_j)

    : MatrixView(mat, offset_i, offset_j, mat.nrRows() - offset_i, mat.nrCols() - offset_j) {

  assert(offset_i < mat.nrCols());

  assert(offset_j < mat.nrRows());

}



template <typename ScalarType, DeviceType device_t>

template <template <typename, DeviceType> class Matrix>

MatrixView<ScalarType, device_t>::MatrixView(Matrix<ScalarType, device_t>& mat, int offset_i,

                                             int offset_j, int ni, int nj)

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(Matrix<Scalar, device_t>& mat, int offset_i, int offset_j,

                                         int ni, int nj)

    : ptr_(mat.ptr(offset_i, offset_j)), ldm_(mat.leadingDimension()), size_(std::make_pair(ni, nj)) {

  assert(ni + offset_i <= mat.nrRows());

  assert(nj + offset_j <= mat.nrCols());

}

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix>

MatrixView<Scalar, device_t>::MatrixView(const Matrix<std::remove_cv_t<Scalar>, device_t>& mat,

                                         int offset_i, int offset_j, int ni, int nj)

    : ptr_(mat.ptr(offset_i, offset_j)), ldm_(mat.leadingDimension()), size_(std::make_pair(ni, nj)) {

  assert(ni + offset_i <= mat.nrRows());

  assert(nj + offset_j <= mat.nrCols());

}



template <typename ScalarType, DeviceType device_t>

template <template <typename, DeviceType> class Matrix>

MatrixView<ScalarType, device_t>& MatrixView<ScalarType, device_t>::operator=(

    const Matrix<ScalarType, device_t>& rhs) {

template <class Scalar, DeviceType device_t>

template <template <class, DeviceType> class Matrix, class Scalar2>

MatrixView<Scalar, device_t>& MatrixView<Scalar, device_t>::operator=(

    const Matrix<Scalar2, device_t>& rhs) {

  static_assert(device_t == CPU, "Copy implemented only on CPU");

  if (nrCols() != rhs.nrCols() || nrRows() != rhs.nrRows()) {

    throw(std::invalid_argument("Matrix size mismatch."));

  }
@@ -155,9 +188,9 @@ MatrixView<ScalarType, device_t>& MatrixView<ScalarType, device_t>::operator=( 
  return *this;

}



template <typename ScalarType, DeviceType device_t>

MatrixView<ScalarType, device_t>& MatrixView<ScalarType, device_t>::operator=(

    const MatrixView<ScalarType, device_t>& rhs) {

template <class Scalar, DeviceType device_t>

MatrixView<Scalar, device_t>& MatrixView<Scalar, device_t>::operator=(const MatrixView& rhs) {

  static_assert(device_t == CPU, "Copy implemented only on CPU");

  if (nrCols() != rhs.nrCols() || nrRows() != rhs.nrRows()) {

    throw(std::invalid_argument("Matrix size mismatch."));

  }
@@ -168,8 +201,8 @@ MatrixView<ScalarType, device_t>& MatrixView<ScalarType, device_t>::operator=( 
  return *this;

}



template <typename ScalarType, DeviceType device_t>

void MatrixView<ScalarType, device_t>::print(std::ostream& out) const {

template <class Scalar, DeviceType device_t>

void MatrixView<Scalar, device_t>::print(std::ostream& out) const {

  out << "\tMatrix view:\n";

  out << "Size: \t" << size_.first << ", " << size_.second << "\n";
@@ -181,47 +214,6 @@ void MatrixView<ScalarType, device_t>::print(std::ostream& out) const { 
  out << "\n" << std::endl;

}



// Methods for returning a pointer to constant matrix view from non const arguments

template <typename ScalarType, DeviceType device_t>

auto inline makeConstantView(const ScalarType* data, const std::pair<int, int> size, const int ld) {

  return std::make_unique<const MatrixView<ScalarType, device_t>>(const_cast<ScalarType*>(data),

                                                                  size, ld);

}



template <typename ScalarType, DeviceType device_t>

auto inline makeConstantView(ScalarType* const data, const int size, const int ld) {

  return makeConstantView<ScalarType, device_t>(data, std::make_pair(size, size), ld);

}



template <typename ScalarType, DeviceType device_t>

auto inline makeConstantView(ScalarType* const data, const int size) {

  return makeConstantView<ScalarType, device_t>(data, std::make_pair(size, size), size);

}



template <template <typename, DeviceType> class Matrix, typename ScalarType, DeviceType device_t>

auto inline makeConstantView(const Matrix<ScalarType, device_t>& mat) {

  return makeConstantView<ScalarType, device_t>(mat.ptr(), mat.size(), mat.leadingDimension());

}



template <template <typename, DeviceType> class Matrix, typename ScalarType, DeviceType device_t>

auto inline makeConstantView(const Matrix<ScalarType, device_t>& mat, const int offset_i,

                             const int offset_j) {

  assert(offset_i < mat.nrCols());

  assert(offset_j < mat.nrRows());

  return makeConstantView<ScalarType, device_t>(

      mat.ptr(offset_i, offset_j), std::make_pair(mat.nrRows() - offset_i, mat.nrCols() - offset_j),

      mat.leadingDimension());

}



template <template <typename, DeviceType> class Matrix, typename ScalarType, DeviceType device_t>

auto inline makeConstantView(const Matrix<ScalarType, device_t>& mat, const int offset_i,

                             const int offset_j, const int ni, const int nj) {

  assert(ni + offset_i <= mat.nrRows());

  assert(nj + offset_j <= mat.nrCols());

  return makeConstantView<ScalarType, device_t>(mat.ptr(offset_i, offset_j), std::make_pair(ni, nj),

                                                mat.leadingDimension());

}



}  // namespace linalg

}  // namespace dca

include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu.hpp

+8 −7
Original line number Diff line number Diff line
@@ -45,6 +45,7 @@ public: 
  using Matrix = typename BaseClass::Matrix;

  using MatrixPair = typename BaseClass::MatrixPair;

  using MatrixView = typename linalg::MatrixView<Real, linalg::CPU>;

  using ConstView = typename linalg::MatrixView<const Real, linalg::CPU>;



public:

  CtintWalker(const Parameters& pars_ref, const Data& /*data*/, Rng& rng_ref, int id = 0);
@@ -74,7 +75,7 @@ private: 
  Real removalProbability();

  void applyRemoval();



  virtual void smallInverse(const MatrixView& in, MatrixView& out, int s);

  virtual void smallInverse(const ConstView& in, MatrixView& out, int s);

  virtual void smallInverse(MatrixView& in_out, int s);



protected:
@@ -306,7 +307,7 @@ void CtintWalker<linalg::CPU, Parameters, Real>::applyInsertion(const MatrixPair 
      continue;

    // update M matrix.

    const auto& R = Rp[s];

    const auto S = linalg::makeConstantView(Sp[s]);

    const ConstView S(Sp[s]);

    auto& M = M_[s];

    const auto& M_Q = M_Q_[s];

    const int m_size = M.nrCols();
@@ -314,7 +315,7 @@ void CtintWalker<linalg::CPU, Parameters, Real>::applyInsertion(const MatrixPair 
    if (not m_size) {

      M.resizeNoCopy(delta);

      auto M_view = MatrixView(M);

      smallInverse(*S, M_view, s);

      smallInverse(S, M_view, s);

      continue;

    }
@@ -325,7 +326,7 @@ void CtintWalker<linalg::CPU, Parameters, Real>::applyInsertion(const MatrixPair 
    M.resize(m_size + delta);

    //  S_tilde = S^-1.

    MatrixView S_tilde(M, m_size, m_size, delta, delta);

    smallInverse(*S, S_tilde, s);

    smallInverse(S, S_tilde, s);



    // R_tilde = - S * R * M

    MatrixView R_tilde(M, m_size, 0, delta, m_size);
@@ -411,7 +412,7 @@ void CtintWalker<linalg::CPU, Parameters, Real>::moveRemovalToEnd() { 
}



template <class Parameters, typename Real>

void CtintWalker<linalg::CPU, Parameters, Real>::smallInverse(const MatrixView& in, MatrixView& out,

void CtintWalker<linalg::CPU, Parameters, Real>::smallInverse(const ConstView& in, MatrixView& out,

                                                              const int s) {

  details::smallInverse(in, out, det_ratio_[s], ipiv_, v_work_);

}

include/dca/phys/dca_step/cluster_solver/ctint/walker/tools/walker_methods.hpp

+1 −1
Original line number Diff line number Diff line
@@ -89,7 +89,7 @@ template <class MatrixA, class MatrixB, typename Real> 
inline void smallInverse(const MatrixA& in, MatrixB& out, const Real det,

                         linalg::Vector<int, linalg::CPU>& ipiv,

                         linalg::Vector<Real, linalg::CPU>& work) {

  static_assert(std::is_same<typename MatrixA::ValueType, Real>::value,

  static_assert(std::is_same<std::remove_cv_t<typename MatrixA::ValueType>, Real>::value,

                "Scalar type MatrixA mismatch.");

  static_assert(std::is_same<typename MatrixB::ValueType, Real>::value,

                "Scalar type MatrixB mismatch.");

test/unit/linalg/matrix_view_test.cpp

+32 −46
Original line number Diff line number Diff line
@@ -28,6 +28,16 @@ TEST(MatrixViewTest, Constructors) { 
    EXPECT_EQ(mat.ptr(), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

  {

    // Construct from const

    const dca::linalg::Matrix<float, dca::linalg::CPU> mat(std::make_pair(5, 4));

    dca::linalg::MatrixView<const float, dca::linalg::CPU> view(mat);

    // view(0, 0) = 4.; Does not compile.

    EXPECT_EQ(mat.nrRows(), view.nrRows());

    EXPECT_EQ(mat.nrCols(), view.nrCols());

    EXPECT_EQ(mat.ptr(), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

  {

    dca::linalg::Matrix<double, dca::linalg::CPU> mat(size);

    const int delta_i(1), delta_j(2);
@@ -37,6 +47,16 @@ TEST(MatrixViewTest, Constructors) { 
    EXPECT_EQ(mat.ptr(delta_i, delta_j), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

  {

    // Construct from const

    const dca::linalg::Matrix<float, dca::linalg::CPU> mat(std::make_pair(10, 12));

    const int delta_i(4), delta_j(6);

    dca::linalg::MatrixView<const float, dca::linalg::CPU> view(mat, delta_i, delta_j);

    EXPECT_EQ(mat.nrRows() - delta_i, view.nrRows());

    EXPECT_EQ(mat.nrCols() - delta_j, view.nrCols());

    EXPECT_EQ(mat.ptr(delta_i, delta_j), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

  {

    dca::linalg::Matrix<double, dca::linalg::CPU> mat(size);

    const int delta_i(0), delta_j(3);
@@ -47,6 +67,17 @@ TEST(MatrixViewTest, Constructors) { 
    EXPECT_EQ(mat.ptr(delta_i, delta_j), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

  {

    // Construct from const

    const dca::linalg::Matrix<int, dca::linalg::CPU> mat(std::make_pair(2, 13));

    const int delta_i(1), delta_j(5);

    const int ni(1), nj(0);

    dca::linalg::MatrixView<const int, dca::linalg::CPU> view(mat, delta_i, delta_j, ni, nj);

    EXPECT_EQ(ni, view.nrRows());

    EXPECT_EQ(nj, view.nrCols());

    EXPECT_EQ(mat.ptr(delta_i, delta_j), view.ptr());

    EXPECT_EQ(mat.leadingDimension(), view.leadingDimension());

  }

}



TEST(MatrixViewTest, ReadWrite) {
@@ -68,7 +99,7 @@ TEST(MatrixViewTest, ReadWrite) { 


TEST(MatrixViewTest, Assignment) {

  dca::linalg::Matrix<int, dca::linalg::CPU> mat(4);

  auto init_func = [](int i, int j) { return i >= 2 ? 1 : 0; };

  auto init_func = [](int i, int /*j*/) { return i >= 2 ? 1 : 0; };

  testing::setMatrixElements(mat, init_func);



  dca::linalg::MatrixView<int, dca::linalg::CPU> upper_left(mat, 0, 0, 2, 2);
@@ -83,48 +114,3 @@ TEST(MatrixViewTest, Assignment) { 
  // Invalid assignment:

  EXPECT_THROW(upper_left = another_size, std::invalid_argument);

}
 


TEST(MatrixViewTest, MakeConstantView) {

  dca::linalg::Matrix<double, dca::linalg::CPU> mat(std::make_pair(4, 2));

  auto init_func = [](int i, int j) { return j + 10. * i; };

  testing::setMatrixElements(mat, init_func);

  {

    const dca::linalg::Matrix<double, dca::linalg::CPU> const_mat(mat);

    auto const_view_ptr = dca::linalg::makeConstantView(const_mat);

    const auto& const_view = *const_view_ptr;

    EXPECT_EQ(const_mat.nrRows(), const_view.nrRows());

    EXPECT_EQ(const_mat.nrCols(), const_view.nrCols());

    EXPECT_EQ(const_mat.ptr(), const_view.ptr());

    EXPECT_EQ(const_mat.leadingDimension(), const_view.leadingDimension());



    for (int j = 0; j < const_view.nrCols(); ++j)

      for (int i = 0; i < const_view.nrCols(); ++i)

        EXPECT_EQ(const_mat(i, j), const_view(i, j));

  }

  {

    const dca::linalg::Matrix<double, dca::linalg::CPU> const_mat(mat);

    auto const_view_ptr = dca::linalg::makeConstantView(const_mat, 1, 0);

    const auto& const_view = *const_view_ptr;

    EXPECT_EQ(const_mat.nrRows() - 1, const_view.nrRows());

    EXPECT_EQ(const_mat.nrCols(), const_view.nrCols());

    EXPECT_EQ(const_mat.ptr(1, 0), const_view.ptr());

    EXPECT_EQ(const_mat.leadingDimension(), const_view.leadingDimension());



    for (int j = 0; j < const_view.nrCols(); ++j)

      for (int i = 0; i < const_view.nrCols(); ++i)

        EXPECT_EQ(const_mat(i + 1, j), const_view(i, j));

  }

  {

    const dca::linalg::Matrix<double, dca::linalg::CPU> const_mat(mat);

    auto const_view_ptr = dca::linalg::makeConstantView(const_mat, 0, 1, 2, 1);

    const auto& const_view = *const_view_ptr;

    EXPECT_EQ(2, const_view.nrRows());

    EXPECT_EQ(1, const_view.nrCols());

    EXPECT_EQ(const_mat.ptr(0, 1), const_view.ptr());

    EXPECT_EQ(const_mat.leadingDimension(), const_view.leadingDimension());



    for (int j = 0; j < const_view.nrCols(); ++j)

      for (int i = 0; i < const_view.nrCols(); ++i)

        EXPECT_EQ(const_mat(i, j + 1), const_view(i, j));

  }

}