1 files changed, 24 insertions, 4 deletions
diff --git a/solve.c b/solve.c
index b9b26a2..64dcb55 100644
--- a/solve.c
+++ b/solve.c
@@ -169,6 +169,26 @@ static int brill_solve_svd(const BDContext *bd, double *mat,
     return 0;
 }
 
+static int brill_export_coeffs(BDPriv *s, const double *coeffs)
+{
+    int alignment = REQ_ALIGNMENT(*coeffs);
+    int nb_coeffs_aligned[2] = { ALIGN(s->nb_coeffs[0], alignment), ALIGN(s->nb_coeffs[0], alignment) };
+
+    int ret, i;
+
+    ret = posix_memalign((void**)&s->coeffs, REQ_ALIGNMENT(char), nb_coeffs_aligned[0] * nb_coeffs_aligned[1] * sizeof(*s->coeffs));
+    if (ret)
+        return -ret;
+
+    memset(s->coeffs, 0, nb_coeffs_aligned[0] * nb_coeffs_aligned[1] * sizeof(*s->coeffs));
+    for (i = 0; i < s->nb_coeffs[1]; i++)
+        memcpy(s->coeffs + i * nb_coeffs_aligned[0], coeffs + i * s->nb_coeffs[0], sizeof(*coeffs) * s->nb_coeffs[0]);
+
+    s->coeffs_stride = nb_coeffs_aligned[0];
+
+    return 0;
+}
+
 /*
  * Solve the equation
  * Δψ + ¼ ψ (∂²q/∂r² + ∂²q/∂z²) = 0
@@ -213,12 +233,12 @@ int bdi_solve(BDContext *bd)
         ret = brill_solve_linear(bd, mat, &coeffs, &rhs);
 
     /* export the result */
-    s->coeffs = coeffs;
-
-fail:
+    ret = brill_export_coeffs(s, coeffs);
     if (ret < 0)
-        free(coeffs);
+        goto fail;
 
+fail:
+    free(coeffs);
     free(rhs);
     free(mat);