/*
 * 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;
}