basis: disentangle the API from the global solver API

6 files changed, 145 insertions, 42 deletions

diff --git a/basis.c b/basis.c
index 0735752..face638 100644
--- a/basis.c
+++ b/basis.c
@@ -21,20 +21,47 @@
  */
 
 #include <math.h>
+#include <stdlib.h>
 
+#include "basis.h"
 #include "internal.h"
 
-static double sb_even_eval(double coord, int idx, double sf)
+
+/* a set of basis functions */
+typedef struct BasisSet {
+    /* evaluate the idx-th basis function at the specified point*/
+    double (*eval)      (const BasisSetContext *s, double coord, int idx);
+    /* evaluate the first derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff1)(const BasisSetContext *s, double coord, int idx);
+    /* evaluate the second derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff2)(const BasisSetContext *s, double coord, int idx);
+
+    void (*eval_multi)(const BasisSetContext *s, double coord,
+                       double **out, int order_start, int order_end);
+    /**
+     * Get the idx-th collocation point for the specified order.
+     * idx runs from 0 to order - 1 (inclusive)
+     */
+    double (*colloc_point)(const BasisSetContext *s, int order, int idx);
+} BasisSet;
+
+struct BasisSetContext {
+    const struct BasisSet *bs;
+    double sf;
+};
+
+static double sb_even_eval(const BasisSetContext *s, double coord, int idx)
 {
-    double val = atan2(sf, coord);
+    double val = atan2(s->sf, coord);
 
     idx *= 2;   // even only
 
     return sin((idx + 1) * val);
 }
 
-static double sb_even_eval_diff1(double coord, int idx, double sf)
+static double sb_even_eval_diff1(const BasisSetContext *s, double coord, int idx)
 {
+    const double sf = s->sf;
     double val = atan2(sf, coord);
 
     idx *= 2;   // even only
@@ -42,8 +69,9 @@ static double sb_even_eval_diff1(double coord, int idx, double sf)
     return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord));
 }
 
-static double sb_even_eval_diff2(double coord, int idx, double sf)
+static double sb_even_eval_diff2(const BasisSetContext *s, double coord, int idx)
 {
+    const double sf = s->sf;
     double val = atan2(sf, coord);
 
     idx *= 2;   // even only
@@ -51,14 +79,14 @@ static double sb_even_eval_diff2(double coord, int idx, double sf)
     return sf * (idx + 1) * (2 * coord * cos((idx + 1) * val) - sf * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord));
 }
 
-static double sb_even_colloc_point(int order, int idx, double sf)
+static double sb_even_colloc_point(const BasisSetContext *s, int order, int idx)
 {
     double t;
 
     idx = order - idx - 1;
 
     t = (idx + 2) * M_PI / (2 * order + 2);
-    return sf / tan(t);
+    return s->sf / tan(t);
 }
 
 /*
@@ -66,9 +94,54 @@ static double sb_even_colloc_point(int order, int idx, double sf)
  * SB(x, n) = sin((n + 1) arccot(|x| / L))
  * They are symmetric wrt origin and decay as odd powers of x in infinity.
  */
-const BasisSet bdi_sb_even_basis = {
+static const BasisSet sb_even_basis = {
     .eval         = sb_even_eval,
     .eval_diff1   = sb_even_eval_diff1,
     .eval_diff2   = sb_even_eval_diff2,
     .colloc_point = sb_even_colloc_point,
 };
+
+double bdi_basis_eval(BasisSetContext *s, enum BSEvalType type, double coord, int order)
+{
+    double (*eval)(const BasisSetContext *s, double coord, int order);
+
+    switch (type) {
+    case BS_EVAL_TYPE_VALUE: eval = s->bs->eval;       break;
+    case BS_EVAL_TYPE_DIFF1: eval = s->bs->eval_diff1; break;
+    case BS_EVAL_TYPE_DIFF2: eval = s->bs->eval_diff2; break;
+    }
+
+    return eval(s, coord, order);
+}
+
+double bdi_basis_colloc_point(BasisSetContext *s, int idx, int order)
+{
+    return s->bs->colloc_point(s, order, idx);
+}
+
+void bdi_basis_free(BasisSetContext **ps)
+{
+    BasisSetContext *s = *ps;
+    free(s);
+    *ps = NULL;
+}
+
+BasisSetContext *bdi_basis_init(enum BasisFamily family, double sf)
+{
+    BasisSetContext *s = calloc(1, sizeof(*s));
+
+    if (!s)
+        return NULL;
+
+    switch (family) {
+    case BASIS_FAMILY_SB_EVEN:
+        s->bs = &sb_even_basis;
+        s->sf = sf;
+        break;
+    default:
+        free(s);
+        return NULL;
+    }
+
+    return s;
+}
diff --git a/basis.h b/basis.h
new file mode 100644
index 0000000..6b6f3c4
--- /dev/null
+++ b/basis.h
@@ -0,0 +1,40 @@
+/*
+ * Copyright 2015 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 BRILL_DATA_BASIS_H
+#define BRILL_DATA_BASIS_H
+
+enum BSEvalType {
+    BS_EVAL_TYPE_VALUE,
+    BS_EVAL_TYPE_DIFF1,
+    BS_EVAL_TYPE_DIFF2,
+};
+
+enum BasisFamily {
+    BASIS_FAMILY_SB_EVEN,
+};
+
+typedef struct BasisSetContext BasisSetContext;
+
+BasisSetContext *bdi_basis_init(enum BasisFamily family, double sf);
+void bdi_basis_free(BasisSetContext **s);
+
+double bdi_basis_eval(BasisSetContext *s, enum BSEvalType type,
+                      double coord, int order);
+double bdi_basis_colloc_point(BasisSetContext *s, int idx, int order);
+
+#endif /* BRILL_DATA_BASIS_H */
diff --git a/eval.c b/eval.c
index 98fd0a7..c8b4b95 100644
--- a/eval.c
+++ b/eval.c
@@ -36,8 +36,7 @@ int bd_eval_psi(const BDContext *bd,
                 double *psi, unsigned int psi_stride)
 {
     BDPriv *s = bd->priv;
-    const double *sf = bd->basis_scale_factor;
-    double (*eval)(double coord, int idx, double sf);
+    enum BSEvalType type;
 
     double *basis_val_rho = NULL, *basis_val_z = NULL;
 
@@ -62,26 +61,26 @@ int bd_eval_psi(const BDContext *bd,
     }
 
     switch (diff_order[0]) {
-    case 0: eval = s->basis[0]->eval;       break;
-    case 1: eval = s->basis[0]->eval_diff1; break;
-    case 2: eval = s->basis[0]->eval_diff2; break;
+    case 0: type = BS_EVAL_TYPE_VALUE; break;
+    case 1: type = BS_EVAL_TYPE_DIFF1; break;
+    case 2: type = BS_EVAL_TYPE_DIFF2; break;
     }
 
     for (int i = 0; i < nb_coords_rho; i++) {
         double rrho = rho[i];
         for (int j = 0; j < s->nb_coeffs[0]; j++)
-            basis_val_rho[i * s->nb_coeffs[0] + j] = eval(rrho, j, sf[0]);
+            basis_val_rho[i * s->nb_coeffs[0] + j] = bdi_basis_eval(s->basis[0], type, rrho, j);
     }
 
     switch (diff_order[1]) {
-    case 0: eval = s->basis[1]->eval;       break;
-    case 1: eval = s->basis[1]->eval_diff1; break;
-    case 2: eval = s->basis[1]->eval_diff2; break;
+    case 0: type = BS_EVAL_TYPE_VALUE; break;
+    case 1: type = BS_EVAL_TYPE_DIFF1; break;
+    case 2: type = BS_EVAL_TYPE_DIFF2; break;
     }
     for (int i = 0; i < nb_coords_z; i++) {
         double zz   = z[i];
         for (int j = 0; j < s->nb_coeffs[1]; j++)
-            basis_val_z[i * s->nb_coeffs[1] + j] = eval(zz, j, sf[1]);
+            basis_val_z[i * s->nb_coeffs[1] + j] = bdi_basis_eval(s->basis[1], type, zz, j);
     }
 
     for (int i = 0; i < nb_coords_z; i++) {
diff --git a/init.c b/init.c
index 5f9e3e6..fa15325 100644
--- a/init.c
+++ b/init.c
@@ -22,6 +22,7 @@
 #include <stdio.h>
 
 #include "brill_data.h"
+#include "basis.h"
 #include "internal.h"
 #include "qfunc.h"
 
@@ -35,8 +36,10 @@ static int brill_init_check_options(BDContext *bd)
     if (!s->qfunc)
         return ret;
 
-    s->basis[0] = &bdi_sb_even_basis;
-    s->basis[1] = &bdi_sb_even_basis;
+    s->basis[0] = bdi_basis_init(BASIS_FAMILY_SB_EVEN, bd->basis_scale_factor[0]);
+    s->basis[1] = bdi_basis_init(BASIS_FAMILY_SB_EVEN, bd->basis_scale_factor[0]);
+    if (!s->basis[0] || !s->basis[1])
+        return -ENOMEM;
 
     s->nb_colloc_points[0] = bd->nb_coeffs[0];
     s->nb_colloc_points[1] = bd->nb_coeffs[1];
@@ -108,6 +111,9 @@ void bd_context_free(BDContext **pbd)
 
     s = bd->priv;
 
+    bdi_basis_free(&s->basis[0]);
+    bdi_basis_free(&s->basis[1]);
+
     bdi_qfunc_free(&s->qfunc);
 
     free(s->coeffs);
diff --git a/internal.h b/internal.h
index ae9932a..78b3741 100644
--- a/internal.h
+++ b/internal.h
@@ -22,6 +22,7 @@
 
 #include "brill_data.h"
 
+#include "basis.h"
 #include "qfunc.h"
 
 #define MAX(x, y) ((x) > (y) ? (x) : (y))
@@ -37,23 +38,9 @@
  */
 #define EPS 1E-08
 
-/* a set of basis functions */
-typedef struct BasisSet {
-    /* evaluate the idx-th basis function at the specified point*/
-    double (*eval)      (double coord, int idx, double sf);
-    /* evaluate the first derivative of the idx-th basis function at the specified point*/
-    double (*eval_diff1)(double coord, int idx, double sf);
-    /* evaluate the second derivative of the idx-th basis function at the specified point*/
-    double (*eval_diff2)(double coord, int idx, double sf);
-    /**
-     * Get the idx-th collocation point for the specified order.
-     * idx runs from 0 to order - 1 (inclusive)
-     */
-    double (*colloc_point)(int order, int idx, double sf);
-} BasisSet;
 
 typedef struct BDPriv {
-    const BasisSet *basis[2];
+    BasisSetContext *basis[2];
     QFuncContext    *qfunc;
 
     int nb_colloc_points[2];
@@ -66,7 +53,6 @@ typedef struct BDPriv {
 #define NB_COEFFS(s)        (s->nb_coeffs[0]        * s->nb_coeffs[1])
 #define NB_COLLOC_POINTS(s) (s->nb_colloc_points[0] * s->nb_colloc_points[1])
 
-extern const BasisSet bdi_sb_even_basis;
 
 
 int bdi_solve(BDContext *bd);
diff --git a/solve.c b/solve.c
index 5803d87..fe47832 100644
--- a/solve.c
+++ b/solve.c
@@ -38,7 +38,6 @@ static int brill_construct_matrix(const BDContext *bd, double *mat,
                                   double *rhs)
 {
     BDPriv *s = bd->priv;
-    const double *sf = bd->basis_scale_factor;
 
     double *basis_val, *basis_dval, *basis_d2val;
     int idx_coeff_rho, idx_coeff_z, idx_grid_rho, idx_grid_z;
@@ -56,11 +55,11 @@ static int brill_construct_matrix(const BDContext *bd, double *mat,
             }
         }
         for (int j = 0; j < s->nb_colloc_points[i]; j++) {
-            double coord = s->basis[i]->colloc_point(s->nb_coeffs[i], j, sf[i]);
+            double coord = bdi_basis_colloc_point(s->basis[i], j, s->nb_coeffs[i]);
             for (int k = 0; k < s->nb_coeffs[i]; k++) {
-                basis[i][0][j * s->nb_coeffs[i] + k] = s->basis[i]->eval      (coord, k, sf[i]);
-                basis[i][1][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff1(coord, k, sf[i]);
-                basis[i][2][j * s->nb_coeffs[i] + k] = s->basis[i]->eval_diff2(coord, k, sf[i]);
+                basis[i][0][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_VALUE, coord, k);
+                basis[i][1][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_DIFF1, coord, k);
+                basis[i][2][j * s->nb_coeffs[i] + k] = bdi_basis_eval(s->basis[i], BS_EVAL_TYPE_DIFF2, coord, k);
             }
         }
     }
@@ -73,10 +72,10 @@ static int brill_construct_matrix(const BDContext *bd, double *mat,
 #define D2BASIS_Z   (basis[1][2][idx_grid_z   * s->nb_coeffs[1] + idx_coeff_z])
 
     for (idx_grid_z = 0; idx_grid_z < s->nb_colloc_points[1]; idx_grid_z++) {
-        double z_val = s->basis[1]->colloc_point(s->nb_coeffs[1], idx_grid_z, sf[1]);
+        double z_val = bdi_basis_colloc_point(s->basis[1], idx_grid_z, s->nb_coeffs[1]);
 
         for (idx_grid_rho = 0; idx_grid_rho < s->nb_colloc_points[0]; idx_grid_rho++) {
-            double x_val = s->basis[0]->colloc_point(s->nb_coeffs[0], idx_grid_rho, sf[0]);
+            double x_val = bdi_basis_colloc_point(s->basis[0], idx_grid_rho, s->nb_coeffs[0]);
             double d2q = bdi_qfunc_eval(s->qfunc, x_val, z_val, 2) +
                          bdi_qfunc_eval(s->qfunc, x_val, z_val, 6);
             int idx_grid = idx_grid_z * s->nb_colloc_points[0] + idx_grid_rho;