Initial commit.

6 files changed, 715 insertions, 0 deletions

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..0edb539
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,19 @@
+CFLAGS = -std=c99 -D_XOPEN_SOURCE=700 -fPIC -g
+LDFLAGS = -shared -lm
+
+TARGET = libmg2d.so
+
+OBJECTS =               \
+    ell_relax.o         \
+    log.o               \
+
+
+all: $(TARGET)
+
+$(TARGET): $(OBJECTS)
+	ld ${LDFLAGS} -o $@ $(OBJECTS)
+
+clean:
+	-rm $(OBJECTS) $(TARGET)
+
+.PHONY: clean test
diff --git a/common.h b/common.h
new file mode 100644
index 0000000..d030ea4
--- /dev/null
+++ b/common.h
@@ -0,0 +1,37 @@
+/*
+ * Copyright 2017 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 MG2D_COMMON_H
+#define MG2D_COMMON_H
+
+#define SQR(x) ((x) * (x))
+#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#define MIN(x, y) ((x) > (y) ? (y) : (x))
+#define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
+
+#include <stdlib.h>
+#include <stdint.h>
+#include <sys/time.h>
+static inline int64_t gettime(void)
+{
+    struct timeval tv;
+    gettimeofday(&tv, NULL);
+    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
+}
+
+#endif /* MG2D_COMMON_H */
diff --git a/ell_relax.c b/ell_relax.c
new file mode 100644
index 0000000..f51c671
--- /dev/null
+++ b/ell_relax.c
@@ -0,0 +1,328 @@
+/*
+ * Copyright 2018 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 <float.h>
+#include <errno.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "common.h"
+#include "ell_relax.h"
+#include "log.h"
+
+#define FD_STENCIL_MAX 1
+
+struct EllRelaxInternal {
+    ptrdiff_t stride;
+    double *u_base;
+    double *rhs_base;
+    double *residual_base;
+    double *diff_coeffs_base[ELL_RELAX_DIFF_COEFF_NB];
+};
+
+static const struct {
+    unsigned int stride_idx;
+    unsigned int is_upper;
+} boundary_def[] = {
+    [ELL_RELAX_BOUNDARY_00] = {
+        .stride_idx = 0,
+        .is_upper   = 0
+    },
+    [ELL_RELAX_BOUNDARY_01] = {
+        .stride_idx = 0,
+        .is_upper   = 1,
+    },
+    [ELL_RELAX_BOUNDARY_10] = {
+        .stride_idx = 1,
+        .is_upper   = 0,
+    },
+    [ELL_RELAX_BOUNDARY_11] = {
+        .stride_idx = 1,
+        .is_upper   = 1,
+    },
+};
+
+static inline void
+derivatives_calc(EllRelaxContext *ctx, double *dst, ptrdiff_t idx, ptrdiff_t stride)
+{
+    const double *u = ctx->u;
+    const double dx = ctx->step[0];
+    const double dy = ctx->step[1];
+
+    dst[ELL_RELAX_DIFF_COEFF_00] = u[idx];
+    dst[ELL_RELAX_DIFF_COEFF_10] = (u[idx + 1]      - u[idx - 1])      / (2.0 * dx);
+    dst[ELL_RELAX_DIFF_COEFF_01] = (u[idx + stride] - u[idx - stride]) / (2.0 * dy);
+
+    dst[ELL_RELAX_DIFF_COEFF_20] = (u[idx + 1]      - 2.0 * u[idx] + u[idx - 1])      / (SQR(dx));
+    dst[ELL_RELAX_DIFF_COEFF_02] = (u[idx + stride] - 2.0 * u[idx] + u[idx - stride]) / (SQR(dy));
+
+    dst[ELL_RELAX_DIFF_COEFF_11] = (u[idx + 1 + stride] - u[idx + stride - 1] - u[idx - stride + 1] + u[idx - stride - 1]) / (4.0 * dx * dy);
+}
+
+static void residual_calc(EllRelaxContext *ctx)
+{
+    EllRelaxInternal *priv = ctx->priv;
+
+    for (int idx1 = 0; idx1 < ctx->domain_size[1]; idx1++)
+        for (int idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
+            ptrdiff_t idx = idx1 * priv->stride + idx0;
+
+            double u[ELL_RELAX_DIFF_COEFF_NB];
+            double res;
+
+            derivatives_calc(ctx, u, idx, priv->stride);
+
+            res = -ctx->rhs[idx];
+            for (int i = 0; i < ARRAY_ELEMS(u); i++)
+                res += u[i] * ctx->diff_coeffs[i][idx];
+
+            ctx->residual[idx] = res;
+        }
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
+        const ptrdiff_t strides[2] = { 1, priv->stride };
+        const int si = boundary_def[i].stride_idx;
+        const ptrdiff_t stride_boundary = strides[si];
+        const ptrdiff_t stride_offset   = strides[!si];
+        const size_t size_boundary      = ctx->domain_size[si];
+        const size_t size_offset        = ctx->domain_size[!si];
+
+        double *res = ctx->residual + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset);
+
+        if (ctx->boundaries[i].type == ELL_RELAX_BC_TYPE_FIXVAL) {
+            for (size_t j = 0; j < size_boundary; j++) {
+                *res = 0.0;
+                res += stride_boundary;
+            }
+        }
+    }
+}
+
+static void boundaries_apply_fixval(double *dst, ptrdiff_t dst_stride,
+                                    const double *src,
+                                    size_t boundary_size, ptrdiff_t boundary_stride,
+                                    int is_upper)
+{
+    for (size_t i = 0; i < boundary_size; i++) {
+        *dst = *src;
+        dst += dst_stride;
+        src++;
+    }
+}
+
+static void boundaries_apply_fixdiff(double *dst, ptrdiff_t dst_stride,
+                                    const double *src,
+                                    size_t boundary_size, ptrdiff_t boundary_stride,
+                                    int is_upper)
+{
+    if (!is_upper)
+        boundary_stride *= -1;
+
+    for (size_t i = 0; i < boundary_size; i++) {
+        for (int i = 1; i <= FD_STENCIL_MAX; i++)
+            dst[boundary_stride * i] = dst[-boundary_stride * i];
+
+        dst += dst_stride;
+    }
+}
+
+static void boundaries_apply(EllRelaxContext *ctx)
+{
+    EllRelaxInternal *priv = ctx->priv;
+    const ptrdiff_t strides[2] = { 1, priv->stride };
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
+        const int si = boundary_def[i].stride_idx;
+        const ptrdiff_t stride_boundary = strides[si];
+        const ptrdiff_t stride_offset   = strides[!si];
+        const size_t size_boundary      = ctx->domain_size[si];
+        const size_t size_offset        = ctx->domain_size[!si];
+
+        double *dst = ctx->u + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset);
+
+        switch (ctx->boundaries[i].type) {
+        case ELL_RELAX_BC_TYPE_FIXVAL:
+            boundaries_apply_fixval(dst, stride_boundary, ctx->boundaries[i].val, size_boundary,
+                                    stride_offset, boundary_def[i].is_upper);
+            break;
+        case ELL_RELAX_BC_TYPE_FIXDIFF:
+            boundaries_apply_fixdiff(dst, stride_boundary, ctx->boundaries[i].val, size_boundary,
+                                     stride_offset, boundary_def[i].is_upper);
+            break;
+        }
+    }
+}
+
+int mg2di_ell_relax_step(EllRelaxContext *ctx)
+{
+    EllRelaxInternal *priv = ctx->priv;
+    const double cfl_fac = (1.0 / 4.0) * ctx->step[0] * ctx->step[1];
+
+    for (int idx1 = 0; idx1 < ctx->domain_size[1]; idx1++)
+        for (int idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) {
+            ptrdiff_t idx = idx1 * priv->stride + idx0;
+
+            ctx->u[idx] += cfl_fac * ctx->residual[idx];
+        }
+
+    boundaries_apply(ctx);
+    residual_calc(ctx);
+
+    return 0;
+}
+
+int mg2di_ell_relax_init(EllRelaxContext *ctx)
+{
+    EllRelaxInternal *priv = ctx->priv;
+    int ret;
+
+    if (ctx->fd_stencil == 0 || ctx->fd_stencil > FD_STENCIL_MAX) {
+        mg2di_log(&ctx->logger, 0, "Invalid finite difference stencil: %zd\n",
+                  ctx->fd_stencil);
+        return -EINVAL;
+    }
+
+    if (ctx->step[0] <= DBL_EPSILON || ctx->step[1] <= DBL_EPSILON) {
+        mg2di_log(&ctx->logger, 0, "Spatial step size too small\n");
+        return -EINVAL;
+    }
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
+        if (ctx->boundaries[i].type == ELL_RELAX_BC_TYPE_FIXDIFF) {
+           for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++)
+               if (ctx->boundaries[i].val[k] != 0.0) {
+                   mg2di_log(&ctx->logger, 0, "Only zero boundary derivative supported\n");
+                   return -ENOSYS;
+               }
+        }
+    }
+
+    boundaries_apply(ctx);
+    residual_calc(ctx);
+
+    return 0;
+}
+
+static int ell_relax_arrays_alloc(EllRelaxContext *ctx, const size_t domain_size[2])
+{
+    EllRelaxInternal *priv = ctx->priv;
+
+    const size_t ghosts = FD_STENCIL_MAX;
+    const size_t size_padded[2] = {
+        domain_size[0] + 2 * ghosts,
+        domain_size[1] + 2 * ghosts,
+    };
+    const size_t stride = size_padded[0];
+    const size_t start_offset = ghosts * stride + ghosts;
+    const size_t arr_size = size_padded[0] * size_padded[1];
+    int ret;
+
+    ret = posix_memalign((void**)&priv->u_base, 32, sizeof(*priv->u_base) * arr_size);
+    if (ret != 0)
+        return -ret;
+    ctx->u        = priv->u_base + start_offset;
+    ctx->u_stride = stride;
+
+    ret = posix_memalign((void**)&priv->rhs_base, 32, sizeof(*priv->rhs_base) * arr_size);
+    if (ret != 0)
+        return -ret;
+    ctx->rhs        = priv->rhs_base + start_offset;
+    ctx->rhs_stride = stride;
+
+    ret = posix_memalign((void**)&priv->residual_base, 32, sizeof(*priv->residual_base) * arr_size);
+    if (ret != 0)
+        return -ret;
+    memset(priv->residual_base, 0, sizeof(*priv->residual_base) * arr_size);
+    ctx->residual        = priv->residual_base + start_offset;
+    ctx->residual_stride = stride;
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) {
+        ret = posix_memalign((void**)&priv->diff_coeffs_base[i], 32,
+                             sizeof(*priv->diff_coeffs_base[i]) * arr_size);
+        if (ret != 0)
+            return -ret;
+        ctx->diff_coeffs[i] = priv->diff_coeffs_base[i] + start_offset;
+    }
+    ctx->diff_coeffs_stride = stride;
+
+    priv->stride = stride;
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) {
+        const size_t boundary_size = domain_size[boundary_def[i].stride_idx];
+        ret = posix_memalign((void**)&ctx->boundaries[i].val, 32,
+                             sizeof(*ctx->boundaries[i].val) * boundary_size);
+        if (ret != 0)
+            return -ret;
+        ctx->boundaries[i].val_len = boundary_size;
+    }
+
+    return 0;
+}
+
+EllRelaxContext *mg2di_ell_relax_alloc(size_t domain_size[2])
+{
+    EllRelaxContext *ctx;
+    EllRelaxInternal *priv;
+    int ret;
+
+    ctx = calloc(1, sizeof(*ctx));
+    if (!ctx)
+        return NULL;
+
+    ctx->priv = calloc(1, sizeof(*ctx->priv));
+    if (!ctx->priv)
+        goto fail;
+    priv = ctx->priv;
+
+    if (!domain_size[0] || !domain_size[1] ||
+        domain_size[0] > SIZE_MAX / domain_size[1])
+        goto fail;
+
+    ret = ell_relax_arrays_alloc(ctx, domain_size);
+    if (ret < 0)
+        goto fail;
+
+    ctx->domain_size[0] = domain_size[0];
+    ctx->domain_size[1] = domain_size[1];
+
+    return ctx;
+fail:
+    mg2di_ell_relax_free(&ctx);
+    return NULL;
+}
+
+void mg2di_ell_relax_free(EllRelaxContext **pctx)
+{
+    EllRelaxContext *ctx = *pctx;
+
+    if (!ctx)
+        return;
+
+    free(ctx->priv->u_base);
+    free(ctx->priv->rhs_base);
+    free(ctx->priv->residual_base);
+    for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs_base); i++)
+        free(ctx->priv->diff_coeffs_base[i]);
+    free(ctx->priv);
+
+    for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++)
+        free(ctx->boundaries[i].val);
+
+    free(ctx);
+    *pctx = NULL;
+}
diff --git a/ell_relax.h b/ell_relax.h
new file mode 100644
index 0000000..f24b1ba
--- /dev/null
+++ b/ell_relax.h
@@ -0,0 +1,258 @@
+/*
+ * Relaxation solver for a 2nd order 2D linear PDE.
+ * Copyright 2018 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 MG2D_ELL_RELAX_H
+#define MG2D_ELL_RELAX_H
+
+/**
+ * The problem being solved is a linear partial differential
+ * equation
+ *
+ *   ∑ C_{ab} ∂_a　∂_b u  + ∑ C_{a} ∂_a u + C u = rhs
+ *  a,b                     a
+ *
+ *  where
+ *      * a and b identify spatial directions and run from 0 to 1
+ *      * u = u(x_a) is the unknown function
+ *      * C_{ab}, C_{a} and C are the coefficients in front of
+ *        the corresponding derivative of unknown function
+ *      * rhs is the right-hand side of the equation
+ * C_{*} and rhs are all (known) functions of space and define the
+ * equation to be solved.
+ */
+
+#include <stddef.h>
+#include <stdint.h>
+
+#include "log.h"
+
+/**
+ * Type of the boundary condition on a given boundary.
+ */
+enum EllRelaxBCType {
+    /**
+     * The value of the unknown function is fixed on the boundary.
+     */
+    ELL_RELAX_BC_TYPE_FIXVAL,
+    /**
+     * The normal derivative of the unkown function is fixed on the boundary.
+     *
+     * TODO: the only supported value of the derivative is currently zero, i.e.
+     * periodic boundary conditions.
+     */
+    ELL_RELAX_BC_TYPE_FIXDIFF,
+};
+
+/**
+ * Location of the boundary.
+ */
+enum EllRelaxBoundaryLoc {
+    /**
+     * Lower coord0, lower coord1.
+     */
+    ELL_RELAX_BOUNDARY_00,
+    /**
+     * Lower coord0, upper coord1.
+     */
+    ELL_RELAX_BOUNDARY_01,
+    /**
+     * Upper coord0, lower coord1.
+     */
+    ELL_RELAX_BOUNDARY_10,
+    /**
+     * Upper coord0, upper coord1.
+     */
+    ELL_RELAX_BOUNDARY_11,
+};
+
+/**
+ * Derivative specifier.
+ */
+enum EllRelaxDiffCoeff {
+    /**
+     * Zeroth derivative, i.e. function value.
+     */
+    ELL_RELAX_DIFF_COEFF_00,
+    /**
+     * First derivative wrt coord1.
+     */
+    ELL_RELAX_DIFF_COEFF_01,
+    /**
+     * First derivative wrt coord0.
+     */
+    ELL_RELAX_DIFF_COEFF_10,
+    /**
+     * Second derivative, once wrt coord0 and once wrt coord1
+     */
+    ELL_RELAX_DIFF_COEFF_11,
+    /**
+     * Second derivative wrt coord1
+     */
+    ELL_RELAX_DIFF_COEFF_02,
+    /**
+     * Second derivative wrt coord0
+     */
+    ELL_RELAX_DIFF_COEFF_20,
+    /**
+     * Total number of allowed derivatives.
+     */
+    ELL_RELAX_DIFF_COEFF_NB,
+};
+
+typedef struct EllRelaxInternal EllRelaxInternal;
+
+/**
+ * Boundary condition definition.
+ */
+typedef struct EllRelaxBoundary {
+    /**
+     * Type of the boundary condition.
+     */
+    enum EllRelaxBCType type;
+    /**
+     * For type = ELL_RELAX_BC_TYPE_FIXVAL:
+     *     Values of the unknown function on the boundary.
+     *
+     * Ignored otherwise.
+     */
+    double *val;
+    /**
+     * Number of elements in val.
+     */
+    size_t  val_len;
+} EllRelaxBoundary;
+
+typedef struct EllRelaxContext {
+    /**
+     * Solver private data, not to be accessed in any way by the caller.
+     */
+    EllRelaxInternal *priv;
+
+    /**
+     * The logging context, to be filled by the caller before mg2di_ell_relax_init().
+     */
+    MG2DLogger logger;
+
+    /**
+     * Size of the solver grid, set by mg2di_ell_relax_alloc().
+     * Read-only for the caller.
+     */
+    size_t domain_size[2];
+
+    /**
+     * Distance between the neighbouring grid points.
+     * Must be set by the caller before mg2di_ell_relax_init().
+     */
+    double step[2];
+
+    /**
+     * Order of the finite difference operators used for approximating derivatives.
+     * Must be set by the caller before mg2di_ell_relax_init().
+     */
+    size_t fd_stencil;
+
+    /**
+     * Boundary specification, indexed by EllRelaxBoundaryLoc.
+     * To be filled by the caller before mg2di_ell_relax_init().
+     */
+    EllRelaxBoundary boundaries[4];
+
+    /**
+     * Values of the unknown function.
+     *
+     * Allocated by the solver in mg2di_ell_relax_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_ell_relax_init() to set the
+     * initial guess.
+     * Afterwards updated in mg2di_ell_relax_step().
+     */
+    double   *u;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t u_stride;
+
+    /**
+     * Values of the right-hand side.
+     *
+     * Allocated by the solver in mg2di_ell_relax_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_ell_relax_init().
+     */
+    double    *rhs;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t  rhs_stride;
+
+    /**
+     * Values of the right-hand side.
+     *
+     * Allocated by the solver in mg2di_ell_relax_alloc(), owned by the solver.
+     * Read-only for the caller. Initialized after mg2di_ell_relax_init(),
+     * afterwards updated in mg2di_ell_relax_step().
+     */
+    double   *residual;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t residual_stride;
+
+    /**
+     * Coefficients C_{*} that define the differential equation.
+     *
+     * Allocated by the solver in mg2di_ell_relax_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_ell_relax_init().
+     */
+    double   *diff_coeffs[ELL_RELAX_DIFF_COEFF_NB];
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t diff_coeffs_stride;
+} EllRelaxContext;
+
+/**
+ * Allocate the solver for the given domain size.
+ *
+ * @param domain_size number of grid points in each direction.
+ *
+ * @return The solver context on success, NULL on failure.
+ */
+EllRelaxContext *mg2di_ell_relax_alloc(size_t domain_size[2]);
+/**
+ * Initialize the solver for use, after all the required fields are filled by
+ * the caller.
+ *
+ * This function may be called multiple times to re-initialize the solver after
+ * certain parameters (e.g. the right-hand side or the equation coefficients)
+ * change.
+ *
+ * @return 0 on success, a negative error code on failure.
+ */
+int              mg2di_ell_relax_init(EllRelaxContext *ctx);
+/**
+ * Free the solver and write NULL to the provided pointer.
+ */
+void             mg2di_ell_relax_free(EllRelaxContext **ctx);
+
+/**
+ * Perform a single relaxation step.
+ *
+ * @return 0 on success, a negative error code on failure.
+ */
+int mg2di_ell_relax_step(EllRelaxContext *ctx);
+
+#endif /* MG2D_ELL_RELAX_H */
diff --git a/log.c b/log.c
new file mode 100644
index 0000000..a940ba6
--- /dev/null
+++ b/log.c
@@ -0,0 +1,40 @@
+/*
+ * Copyright 2014-2017 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/>.
+ */
+
+/**
+ * @file
+ * logging code
+ */
+
+#include <stdarg.h>
+#include <stdio.h>
+
+#include "log.h"
+
+void mg2di_log_default_callback(MG2DLogger *log, int level,
+                                const char *fmt, va_list vl)
+{
+    vfprintf(stderr, fmt, vl);
+}
+
+void mg2di_log(MG2DLogger *log, int level, const char *fmt, ...)
+{
+    va_list vl;
+    va_start(vl, fmt);
+    log->log(log, level, fmt, vl);
+    va_end(vl);
+}
diff --git a/log.h b/log.h
new file mode 100644
index 0000000..99fa9d1
--- /dev/null
+++ b/log.h
@@ -0,0 +1,33 @@
+/*
+ * Logging infrastructure
+ * Copyright 2017 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 MG2D_LOG_H
+#define MG2D_LOG_H
+
+#include <stdarg.h>
+
+typedef struct MG2DLogger {
+    void *opaque;
+    void (*log)(struct MG2DLogger *log, int level, const char *fmt, va_list vl);
+} MG2DLogger;
+
+void mg2di_log(MG2DLogger *log, int level, const char *fmt, ...);
+void mg2di_log_default_callback(MG2DLogger *log, int level,
+                                const char *fmt, va_list vl);
+
+#endif // MG2D_LOG_H