Bug report: incorrect dot product value for non-contiguous arrays

[pavel-po]

Description

nda::blas::dot gives incorrect results when nda arrays are non-contiguous in memory.

Steps to Reproduce

A patch to the existing linear algebra test demonstrates the issue by comparison of the dot products of non-contiguous arrays with dot products of contiguous array copies.

diff --git a/test/c++/nda_linear_algebra.cpp b/test/c++/nda_linear_algebra.cpp
index 3404c2a..837250e 100644
--- a/test/c++/nda_linear_algebra.cpp
+++ b/test/c++/nda_linear_algebra.cpp
@@ -60,6 +60,28 @@ TEST(Vector, Dot2) { //NOLINT
   EXPECT_COMPLEX_NEAR(nda::blas::dotc(v, v), 1);
 }

+TEST(Vector, Dot3) { //NOLINT
+
+  // Test for non-contiguous arrays
+  size_t n = 3;
+  nda::array<std::complex<double>, 4> v_4d(n,n,n,n);
+  for(size_t i = 0; i < n; i++)
+  for(size_t j = 0; j < n; j++)
+  for(size_t k = 0; k < n; k++)
+  for(size_t l = 0; l < n; l++)
+      v_4d(i,j,k,l) = i+j+k+l;
+
+
+  auto v_3d_noncont = v_4d(nda::range::all, 1, nda::range::all, nda::range::all);
+  auto v_3d_noncont_1d = nda::reshaped_view(v_3d_noncont, std::array<long,1>{n*n*n});
+
+  nda::array<std::complex<double>, 3> v_3d_copy = nda::make_regular(v_3d_noncont);
+  auto v_3d_copy_1d = nda::reshaped_view(v_3d_copy, std::array<long,1>{n*n*n});
+
+  EXPECT_COMPLEX_NEAR(nda::blas::dot(v_3d_noncont_1d, v_3d_noncont_1d), nda::blas::dot(v_3d_copy_1d, v_3d_copy_1d));
+
+}
+
 // ==============================================================

 template <typename T, typename L1, typename L2, typename L3>

Versions

The issue is reproduced on Rusty with the following NDA version used in BeyondDFT:

#define NDA_VERSION "1.1.0"
#define NDA_VERSION_MAJOR "1"
#define NDA_VERSION_MINOR "1"
#define NDA_VERSION_PATCH "0"
#define NDA_GIT_HASH "af14bd53403a55a977e9345da4c2867d72cfc31c"

[Thoemi09]

Thank you @pavel-po for opening this issue!

I think the problem is already the reshaping of the array, i.e.

auto v_3d_noncont_1d = nda::reshaped_view(v_3d_noncont, std::array<long,1>{n*n*n});

should not be allowed since v_3d_noncont is not contiguous in memory.

If you compile the code in debug mode, i.e. with -DCMAKE_BUILD_TYPE=Debug, then you should in fact get an error that the reshaping fails (see nda/layout_transforms.hpp).

[Wentzell]

Yes, as detailed also in the nda::reshape documentation the input needs to be contiguous. The runtime checks are currently enabled only in debug mode as @Thoemi09 pointed out. Closing this.