Initial commit.

1 files changed, 139 insertions, 0 deletions

diff --git a/src/pssolve.h b/src/pssolve.h
new file mode 100644
index 0000000..e6a4c1a
--- /dev/null
+++ b/src/pssolve.h
@@ -0,0 +1,139 @@
+/*
+ * Pseudospectral 2nd order 2D linear PDE solver
+ * Copyright (C) 2016 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 MD_PSSOLVE_H
+#define MD_PSSOLVE_H
+
+#include "common.h"
+
+#if HAVE_OPENCL
+#include <cl.h>
+#else
+typedef void* cl_context;
+typedef void* cl_command_queue;
+#endif
+
+#include <stdint.h>
+
+#include "basis.h"
+#include "threadpool.h"
+
+enum PSSolveDiffOrder {
+    PSSOLVE_DIFF_ORDER_00,
+    PSSOLVE_DIFF_ORDER_10,
+    PSSOLVE_DIFF_ORDER_01,
+    PSSOLVE_DIFF_ORDER_11,
+    PSSOLVE_DIFF_ORDER_20,
+    PSSOLVE_DIFF_ORDER_02,
+    PSSOLVE_DIFF_ORDER_NB,
+};
+
+typedef struct PSSolvePriv PSSolvePriv;
+
+typedef struct PSSolveContext {
+    /**
+     * Solver private data, not to be touched by the caller.
+     */
+    PSSolvePriv *priv;
+
+    /**
+     * Number of equations/unknown functions in the set.
+     * Set by md_pssolve_context_alloc().
+     */
+    unsigned int nb_equations;
+
+    /**
+     * The basis sets.
+     *
+     * basis[i][j] is the basis set used for i-th variable in j-th direction.
+     *
+     * The array is allocated by md_pssolve_context_alloc(), must be filled by
+     * by the caller before md_pssolve_context_init().
+     */
+    const MDBasisSetContext *(*basis)[2];
+
+    /**
+     * Order of the solver.
+     *
+     * solve_order[i][j] is the order of the solver (i.e. the number of the
+     * basis functions used) for i-th variable in j-th direction.
+     *
+     * Allocated by md_pssolve_context_alloc(), must be filled by the caller
+     * before md_pssolve_context_init().
+     */
+    unsigned int (*solve_order)[2];
+
+    /**
+     * Locations of the collocation points. The equation coefficients passed to
+     * md_pssolve_solve() should be evaluated at those grid positions.
+     *
+     * colloc_grid[i][j] is an array of length solve_order[i][j] and contains
+     * the collocation points for the i-th variable in the j-th direction.
+     *
+     * Set by the solver after md_pssolve_context_init().
+     */
+    double *(*colloc_grid)[2];
+
+    /**
+     * The thread pool used for multithreaded execution. May be set by the
+     * caller before md_pssolve_context_init(), otherwise a single thread will
+     * be used.
+     */
+    ThreadPoolContext *tp;
+
+    cl_context       ocl_ctx;
+    cl_command_queue ocl_queue;
+
+    uint64_t lu_solves_count;
+    uint64_t lu_solves_time;
+
+    uint64_t cg_solve_count;
+    uint64_t cg_iter_count;
+    uint64_t cg_time_total;
+
+    uint64_t construct_matrix_count;
+    uint64_t construct_matrix_time;
+} PSSolveContext;
+
+/**
+ * Allocate a new solver.
+ */
+int md_pssolve_context_alloc(PSSolveContext **ctx, unsigned int nb_equations);
+
+/**
+ * Initialize the solver for use after all the context options have been set.
+ */
+int md_pssolve_context_init(PSSolveContext *ctx);
+
+/**
+ * Free the solver and all its internal state.
+ */
+void md_pssolve_context_free(PSSolveContext **ctx);
+
+/**
+ * Solve a second order linear PDE in 2D with a pseudospectral method.
+ *
+ * @param eq_coeffs the equation coefficients.
+ * @param rhs the right-hand side of the equation at the collocation points.
+ * @param coeffs the spectral coefficients of the solution will be written here.
+ */
+int md_pssolve_solve(PSSolveContext *ctx,
+                     const double *(**eq_coeffs)[PSSOLVE_DIFF_ORDER_NB],
+                     const double *rhs, double *coeffs);
+
+#endif /* MD_PSSOLVE_H */