From 95867280ab72ce6ade22f2cf7fcbc5b47b6cc95d Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Fri, 11 Mar 2016 20:00:45 +0100
Subject: solve: export aligned coefficients

---
 solve.c | 28 ++++++++++++++++++++++++----
 1 file changed, 24 insertions(+), 4 deletions(-)

(limited to 'solve.c')

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);
 
-- 
cgit v1.2.3