Initial commit.

4 files changed, 380 insertions, 0 deletions

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..4a7040b
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,28 @@
+TARGET = libndarray.so
+
+CFLAGS = -std=c99 -D_XOPEN_SOURCE=700 -fPIC -g -I.
+TARGET_LDFLAGS = -Wl,--version-script=libndarray.v -shared -lm -pthread
+TEST_LIBS = -lm -lpthread
+CC = cc
+
+OBJS =                       \
+       ndarray.o             \
+
+all: $(TARGET)
+
+$(TARGET): $(OBJS)
+	$(CC) ${TARGET_LDFLAGS} -o $@ $(OBJS)
+
+%.o: %.c
+	$(CC) $(CFLAGS) -MMD -MF $(@:.o=.d) -MT $@ -c -o $@ $<
+
+%.o: %.asm
+	nasm -f elf -m amd64 -M $< > $(@:.o=.d)
+	nasm -f elf -m amd64 -o $@ $<
+
+clean:
+	-rm -f *.o *.d *.pyc $(TARGET)
+
+-include $(OBJS:.o=.d)
+
+.PHONY: clean
diff --git a/libndarray.v b/libndarray.v
new file mode 100644
index 0000000..90335b2
--- /dev/null
+++ b/libndarray.v
@@ -0,0 +1,4 @@
+LIBNDARRAY_0 {
+    global: ndarray_*;
+    local: *;
+};
diff --git a/ndarray.c b/ndarray.c
new file mode 100644
index 0000000..d4d4a18
--- /dev/null
+++ b/ndarray.c
@@ -0,0 +1,272 @@
+/*
+ * N-dimensional strided array
+ * Copyright 2018-2020 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <errno.h>
+#include <limits.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "ndarray.h"
+
+struct NDArrayInternal {
+    size_t     alloc_size;
+    void      *data;
+
+    size_t    *shape;
+    ptrdiff_t *stride;
+};
+
+int ndarray_alloc(NDArray **result, unsigned int dims,
+                  const size_t * const size, unsigned int flags)
+{
+    NDArray *a;
+    int ret;
+
+    if (!dims || dims >= INT_MAX)
+        return -EINVAL;
+
+    a = calloc(1, sizeof(*a));
+    if (!a)
+        return -ENOMEM;
+
+    a->priv = calloc(1, sizeof(*a->priv));
+    if (!a->priv) {
+        free(a);
+        return -ENOMEM;
+    }
+
+    a->priv->shape = calloc(dims, sizeof(*a->priv->shape));
+    if (!a->priv->shape) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memcpy(a->priv->shape, size, dims * sizeof(*a->priv->shape));
+
+    a->priv->stride = calloc(dims, sizeof(*a->priv->stride));
+    if (!a->priv->stride) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+
+    a->priv->stride[dims - 1] = 1;
+    for (int i = (int)dims - 2; i >= 0; i--)
+        a->priv->stride[i] = size[i + 1] * a->priv->stride[i + 1];
+
+    a->priv->alloc_size = a->priv->stride[0] * size[0] * sizeof(*a->data);
+    ret = posix_memalign((void**)&a->priv->data, 64, a->priv->alloc_size);
+    if (ret != 0) {
+        ret = -ret;
+        goto fail;
+    }
+
+    if (flags & NDARRAY_ALLOC_ZERO)
+        memset(a->priv->data, 0, a->priv->alloc_size);
+
+    a->dims   = dims;
+
+    a->data   = a->priv->data;
+    a->shape  = a->priv->shape;
+    a->stride = a->priv->stride;
+
+    *result = a;
+
+    return 0;
+fail:
+    ndarray_free(&a);
+    return ret;
+}
+
+int ndarray_slice(NDArray **result, NDArray *src,
+                  const NDASlice * const slice)
+{
+    NDArray *a;
+    ptrdiff_t offset;
+    int ret;
+
+    a = calloc(1, sizeof(*a));
+    if (!a)
+        return -ENOMEM;
+
+    a->priv = calloc(1, sizeof(*a->priv));
+    if (!a->priv) {
+        free(a);
+        return -ENOMEM;
+    }
+
+    a->priv->shape = calloc(src->dims, sizeof(*a->priv->shape));
+    if (!a->priv->shape) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memcpy(a->priv->shape, src->priv->shape, src->dims * sizeof(*a->priv->shape));
+
+    a->priv->stride = calloc(src->dims, sizeof(*a->priv->stride));
+    if (!a->priv->stride) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memcpy(a->priv->stride, src->priv->stride, src->dims * sizeof(*a->priv->shape));
+
+    offset = 0;
+    for (unsigned int i = 0; i < src->dims; i++) {
+        const NDASlice *s = &slice[i];
+        ptrdiff_t start, end, step;
+
+        if (s->start == PTRDIFF_MAX) {
+            start = 0;
+            end   = src->shape[i];
+            step  = 1;
+        } else {
+            start = s->start;
+            end   = s->end;
+            step  = s->step;
+
+            if (step < 1) {
+                ret = step ? -ENOSYS : -EINVAL;
+                goto fail;
+            }
+
+            if (start < 0)
+                start = src->shape[i] + start;
+            if (end <= 0)
+                end = src->shape[i] + end;
+
+            if (start >= end || end > src->shape[i]) {
+                ret = -EINVAL;
+                goto fail;
+            }
+        }
+
+        offset += src->stride[i] * start;
+
+        a->priv->stride[i] *= step;
+        a->priv->shape[i]   = (end - start) / step;
+    }
+
+    a->dims   = src->dims;
+
+    a->data   = src->data + offset;
+    a->shape  = a->priv->shape;
+    a->stride = a->priv->stride;
+
+    *result = a;
+    return 0;
+
+fail:
+    ndarray_free(&a);
+    return ret;
+}
+
+int ndarray_wrap(NDArray **result, unsigned int dims,
+                 const size_t * const size, double *data,
+                 const ptrdiff_t * const stride)
+{
+    NDArray *a;
+    int ret;
+
+    if (!dims || dims >= INT_MAX)
+        return -EINVAL;
+
+    a = calloc(1, sizeof(*a));
+    if (!a)
+        return -ENOMEM;
+
+    a->priv = calloc(1, sizeof(*a->priv));
+    if (!a->priv) {
+        free(a);
+        return -ENOMEM;
+    }
+
+    a->priv->shape = calloc(dims, sizeof(*a->priv->shape));
+    if (!a->priv->shape) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memcpy(a->priv->shape, size, dims * sizeof(*a->priv->shape));
+
+    a->priv->stride = calloc(dims, sizeof(*a->priv->stride));
+    if (!a->priv->stride) {
+        ret = -ENOMEM;
+        goto fail;
+    }
+    memcpy(a->priv->stride, stride, dims * sizeof(*a->priv->shape));
+
+    a->dims   = dims;
+
+    a->data   = data;
+    a->shape  = a->priv->shape;
+    a->stride = a->priv->stride;
+
+    *result = a;
+
+    return 0;
+fail:
+    ndarray_free(&a);
+    return ret;
+}
+
+void ndarray_free(NDArray **pa)
+{
+    NDArray *a = *pa;
+
+    if (!a)
+        return;
+
+    if (a->priv) {
+        free(a->priv->data);
+        free(a->priv->shape);
+        free(a->priv->stride);
+    }
+    free(a->priv);
+
+    free(a);
+    *pa = NULL;
+}
+
+static int copy_axis(NDArray *dst, const NDArray *src, unsigned int axis,
+                     ptrdiff_t offset_dst, ptrdiff_t offset_src)
+{
+    if (dst->shape[axis] != src->shape[axis])
+        return -EINVAL;
+
+    if (axis == dst->dims - 1) {
+        if (dst->stride[axis] == 1 && src->stride[axis] == 1)
+            memcpy(dst->data + offset_dst, src->data + offset_src, sizeof(*dst->data) * dst->shape[axis]);
+        else {
+            for (size_t idx = 0; idx < dst->shape[axis]; idx++)
+                dst->data[offset_dst + idx * dst->stride[axis]] = src->data[offset_src + idx * src->stride[axis]];
+        }
+        return 0;
+    }
+
+    for (size_t idx = 0; idx < dst->shape[axis]; idx++)
+        copy_axis(dst, src, axis + 1, offset_dst + idx * dst->stride[axis], offset_src + idx * src->stride[axis]);
+
+    return 0;
+}
+
+int ndarray_copy(NDArray *dst, const NDArray *src)
+{
+    const unsigned int dims = src->dims;
+
+    if (dims != dst->dims)
+        return -EINVAL;
+
+    return copy_axis(dst, src, 0, 0, 0);
+}
diff --git a/ndarray.h b/ndarray.h
new file mode 100644
index 0000000..16ba59e
--- /dev/null
+++ b/ndarray.h
@@ -0,0 +1,76 @@
+/*
+ * N-dimensional strided array
+ * Copyright 2018-2020 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef NDARRAY_NDARRAY_H
+#define NDARRAY_NDARRAY_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+/**
+ * TODO:
+ * custom allocators
+ * destructors
+ * other data types than double
+ */
+
+typedef struct NDASlice {
+    ptrdiff_t start;
+    ptrdiff_t end;
+    ptrdiff_t step;
+} NDASlice;
+
+#define NDASLICE(start, end, step) ((NDASlice){ start, end, step })
+#define NDASLICE_NULL ((NDASlice){ PTRDIFF_MAX })
+
+typedef struct NDArrayInternal NDArrayInternal;
+
+typedef struct NDArray {
+    NDArrayInternal *priv;
+
+    unsigned int     dims;
+
+    double          *data;
+    const size_t    *shape;
+    const ptrdiff_t *stride;
+} NDArray;
+
+#define NDA_IDX1D(arr, x)       ((arr)->stride[0] * (x))
+#define NDA_IDX2D(arr, x, y)    ((arr)->stride[0] * (y) + (arr)->stride[1] * (x))
+#define NDA_IDX3D(arr, x, y, z) ((arr)->stride[0] * (z) + (arr)->stride[1] * (y) + (arr)->stride[0] * (x))
+
+#define NDA_PTR1D(arr, x)       ((arr)->data + NDA_IDX1D(arr, x))
+#define NDA_PTR2D(arr, x, y)    ((arr)->data + NDA_IDX2D(arr, x, y))
+#define NDA_PTR3D(arr, x, y, z) ((arr)->data + NDA_IDX2D(arr, x, y, z))
+
+#define NDARRAY_ALLOC_ZERO (1 << 0)
+
+int ndarray_alloc(NDArray **arr, unsigned int dims,
+                  const size_t * const size, unsigned int flags);
+void ndarray_free(NDArray **arr);
+
+int ndarray_wrap(NDArray **arr, unsigned int dims,
+                 const size_t * const size, double *data,
+                 const ptrdiff_t * const strides);
+
+int ndarray_slice(NDArray **dst, NDArray *src,
+                  const NDASlice * const slice);
+
+int ndarray_copy(NDArray *dst, const NDArray *src);
+
+#endif // NDARRAY_NDARRAY_H