Initial commit.

The following code is present:
* the basis API
* the BiCGSTAB solver
* the pseudospectral linear system solver
* helper APIs:
        - threadpool
        - logging
        - cpuid

1 files changed, 139 insertions, 0 deletions

diff --git a/tests/pssolve.c b/tests/pssolve.c
new file mode 100644
index 0000000..94afadd
--- /dev/null
+++ b/tests/pssolve.c
@@ -0,0 +1,139 @@
+/*
+ * 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/>.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "basis.h"
+#include "common.h"
+#include "log.h"
+#include "pssolve.h"
+
+#define N_X 8
+#define N_Z 4
+#define FUNC_X 4
+#define FUNC_Z 2
+#define TOL 5e-15
+static int scalar0(void)
+{
+    BasisSetContext *basis = NULL;
+    PSSolveContext    *ctx = NULL;
+
+    double eq_coeffs00[PSSOLVE_DIFF_ORDER_NB][N_X * N_Z];
+    const double *eq_coeffs0[PSSOLVE_DIFF_ORDER_NB];
+    const double *(*eq_coeffs)[PSSOLVE_DIFF_ORDER_NB];
+    double rhs[N_X * N_Z];
+    double coeffs[N_X * N_Z];
+    int ret;
+
+    ret = tdi_basis_init(&basis, BASIS_FAMILY_SB_EVEN, 1.0);
+    if (ret < 0)
+        goto finish;
+
+    ret = tdi_pssolve_context_alloc(&ctx, 1);
+    if (ret < 0)
+        return ret;
+
+    ctx->basis[0][0] = basis;
+    ctx->basis[0][1] = basis;
+
+    ctx->solve_order[0][0] = N_X;
+    ctx->solve_order[0][1] = N_Z;
+
+    ctx->logger.log = tdi_log_default_callback;
+
+    ret = tdi_pssolve_context_init(ctx);
+    if (ret < 0)
+        goto finish;
+
+    for (int j = 0; j < N_Z; j++) {
+        const double z = ctx->colloc_grid[0][1][j];
+
+        for (int i = 0; i < N_X; i++) {
+            const double x = ctx->colloc_grid[0][0][i];
+
+            const int idx = j * N_X + i;
+
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_20][idx] = 1.0;
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_02][idx] = 1.0;
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_10][idx] = 0.0;
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_01][idx] = 0.0;
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_11][idx] = 0.0;
+            eq_coeffs00[PSSOLVE_DIFF_ORDER_00][idx] = 0.0;
+
+            rhs[idx] = tdi_basis_eval(basis, BS_EVAL_TYPE_DIFF2, x, FUNC_X) *
+                       tdi_basis_eval(basis, BS_EVAL_TYPE_VALUE, z, FUNC_Z) +
+                       tdi_basis_eval(basis, BS_EVAL_TYPE_VALUE, x, FUNC_X) *
+                       tdi_basis_eval(basis, BS_EVAL_TYPE_DIFF2, z, FUNC_Z);
+        }
+    }
+    for (int i = 0; i < PSSOLVE_DIFF_ORDER_NB; i++)
+        eq_coeffs0[i] = eq_coeffs00[i];
+    eq_coeffs = &eq_coeffs0;
+
+    ret = tdi_pssolve_solve(ctx, &eq_coeffs, rhs, coeffs);
+    if (ret < 0)
+        goto finish;
+
+    for (int j = 0; j < N_Z; j++)
+        for (int i = 0; i < N_X; i++) {
+            const int idx = j * N_X + i;
+            const double val = (i == FUNC_X && j == FUNC_Z) ? 1.0 : 0.0;
+
+            if (fabs(coeffs[idx] - val) > TOL) {
+                fprintf(stderr, "unexpected value %g at %d/%d\n",
+                        coeffs[idx], i, j);
+                ret = -1;
+                goto finish;
+            }
+        }
+
+finish:
+    tdi_pssolve_context_free(&ctx);
+    tdi_basis_free(&basis);
+
+    return ret;
+}
+#undef N
+#undef FUNC_X
+#undef FUNC_Z
+#undef TOL
+
+static const struct {
+    const char *name;
+    int (*test)(void);
+} tests[] = {
+    { "scalar0", scalar0 },
+};
+
+int main(void)
+{
+    int ret = 0;
+
+    for (int i = 0; i < ARRAY_ELEMS(tests); i++) {
+        fprintf(stderr, "executing test '%s'\n", tests[i].name);
+        ret = tests[i].test();
+        if (ret < 0) {
+            fprintf(stderr, "test '%s' failed\n", tests[i].name);
+            return -1;
+        }
+        fprintf(stderr, "test '%s' succeeded\n", tests[i].name);
+    }
+
+    return 0;
+}