diff options
Diffstat (limited to 'solve.c')
-rw-r--r-- | solve.c | 28 |
1 files changed, 24 insertions, 4 deletions
@@ -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); |