From 1c05a006eb2121ab45d65ff9e82a8e5370563c03 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 25 Jun 2020 16:34:58 +0200 Subject: Initial commit. --- Makefile | 28 ++++++ libndarray.v | 4 + ndarray.c | 272 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ndarray.h | 76 +++++++++++++++++ 4 files changed, 380 insertions(+) create mode 100644 Makefile create mode 100644 libndarray.v create mode 100644 ndarray.c create mode 100644 ndarray.h 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 + * + * 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 . + */ + +#include +#include +#include +#include +#include + +#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 + * + * 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 . + */ + +#ifndef NDARRAY_NDARRAY_H +#define NDARRAY_NDARRAY_H + +#include +#include + +/** + * 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 -- cgit v1.2.3