/*
 * Copyright 2019 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 MG2D_RESIDUAL_CALC_H
#define MG2D_RESIDUAL_CALC_H

#include <stddef.h>

#include <threadpool.h>

typedef struct ResidualCalcInternal ResidualCalcInternal;

/**
 * Context for computing the residual, allocated by mg2di_residual_calc_alloc(),
 * freed by mg2di_residual_calc_free()
 */
typedef struct ResidualCalcContext {
    ResidualCalcInternal *priv;

    /**
     * The thread pool, must be set by the caller before
     * mg2di_residual_calc_init().
     */
    TPContext *tp;

    /**
     * FD stencil size, must be set by the caller before
     * mg2di_residual_calc_init().
     */
    size_t fd_stencil;

    /**
     * Flags indicating supported CPU features,  must be set by the caller
     * before mg2di_residual_calc_init().
     */
    int cpuflags;
} ResidualCalcContext;

/**
 * Allocate and retur a new ResidualCalcContext.
 */
ResidualCalcContext *mg2di_residual_calc_alloc(void);
/**
 * Free a ResidualCalcContext and write NULL into the supplied pointer.
 */
void mg2di_residual_calc_free(ResidualCalcContext **ctx);

/**
 * Reinitialize the context for updated parameters. Must be called at least once
 * before mg2di_residual_calc().
 */
int mg2di_residual_calc_init(ResidualCalcContext *ctx);

/**
 * Calculate the residual.
 *
 * @param diff_coeffs Array containing all PDE coefficients.
 * @param diff_coeffs_stride Distance between adjacent lines for a given PDE
 *                           coefficient.
 * @param diff_coeffs_offset Distance between blocks containing different PDE
 *   coefficients. I.e. diff_coeffs[i * diff_coeffs_offset + j * diff_coeffs_stride + k]
 *   is the value of the i-th PDE coefficient (MG2D_DIFF_COEFF_*) at the
 *   (j, k)th gridpoint.
 *
 * @param reflect Flags indicating if dst should be extended by reflection.
 *      If the 'MG2D_BOUNDARY_x'th bit is set in reflect, then dst should be
 *      reflected reflect_dist points beyond its bounds in the corresponding
 *      direction.
 * @param reflect_dist number of points to fill by reflection
 */
int mg2di_residual_calc(ResidualCalcContext *ctx, size_t size[2],
                        double      *residual_max,
                        double               *dst, ptrdiff_t dst_stride,
                        const double           *u, ptrdiff_t u_stride,
                        const double         *rhs, ptrdiff_t rhs_stride,
                        const double *diff_coeffs, ptrdiff_t diff_coeffs_stride,
                                                   ptrdiff_t diff_coeffs_offset,
                        double u_mult, double res_mult,
                        int reflect, size_t reflect_dist);

#endif // MG2D_RESIDUAL_CALC_H