TRIQS / nda

C++ library for multi-dimensional arrays
https://triqs.github.io/nda
Other
13 stars 11 forks source link

Matrix vector product returns different result whether matrix view or regular matrix is used #43

Closed dominikkiese closed 1 year ago

dominikkiese commented 1 year ago

I added the following test to tests/c++/nda_functions.cpp

TEST(Array, Permuted_view_matmul) { //NOLINT

  // generate some dummy data
  nda::array<double, 3> A = nda::rand(3, 3, 3);
  nda::vector<double> v   = nda::rand(3);

  // build permuted array
  auto B = make_regular(nda::permuted_indices_view<nda::encode(std::array<int, 3>{1, 2, 0})>(A));

  for (auto k : range(3)) {
    auto Amat = matrix<double>{A(_, _, k)}; // copy for contiguous memory
    auto Bmat = matrix_view<double>{B(k, _, _)}; // stride layout should be contiguous
    EXPECT_EQ_ARRAY(Amat, Bmat);
    EXPECT_EQ_ARRAY(Amat * v, Bmat * v);
  }
}

If we replace auto Bmat = matrix_view<double>{B(k, _, _)}; with auto Bmat = matrix<double>{B(k, _, _)}; the test passes. Note that only the second expect is triggered, the data itself seems to be the same. Currently checked out 5d9c817 on sym_grp_pr branch.

dominikkiese commented 1 year ago

The problem with the above code is the following: make_regular will copy the stride layout of the permuted array, such that Bmat will not have a stride of dimension 1 (exception triggered in debug mode). This can be avoided by replacing make_regular with the constructor for array<double, 3>.