From 1c05a006eb2121ab45d65ff9e82a8e5370563c03 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 25 Jun 2020 16:34:58 +0200 Subject: Initial commit. --- ndarray.c | 272 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) create mode 100644 ndarray.c (limited to 'ndarray.c') 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); +} -- cgit v1.2.3