text:
/* matrix/gsl_matrix_uchar.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_UCHAR_H__
#define __GSL_MATRIX_UCHAR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uchar.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned char * data;
  gsl_block_uchar * block;
  int owner;
} gsl_matrix_uchar;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_view;

typedef _gsl_matrix_uchar_view gsl_matrix_uchar_view;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_const_view;

typedef const _gsl_matrix_uchar_const_view gsl_matrix_uchar_const_view;

/* Allocation */

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_block (gsl_block_uchar * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_matrix (gsl_matrix_uchar * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_uchar *
gsl_vector_uchar_alloc_row_from_matrix (gsl_matrix_uchar * m,
                                        const size_t i);

GSL_FUN gsl_vector_uchar *
gsl_vector_uchar_alloc_col_from_matrix (gsl_matrix_uchar * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_uchar_free (gsl_matrix_uchar * m);

/* Views */

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_submatrix (gsl_matrix_uchar * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_row (gsl_matrix_uchar * m, const size_t i);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_column (gsl_matrix_uchar * m, const size_t j);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_diagonal (gsl_matrix_uchar * m);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_superdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subrow (gsl_matrix_uchar * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uchar_view
gsl_matrix_uchar_subcolumn (gsl_matrix_uchar * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_array (unsigned char * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_array_with_tda (unsigned char * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector (gsl_vector_uchar * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector_with_tda (gsl_vector_uchar * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_submatrix (const gsl_matrix_uchar * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_row (const gsl_matrix_uchar * m,
                            const size_t i);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_column (const gsl_matrix_uchar * m,
                               const size_t j);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_diagonal (const gsl_matrix_uchar * m);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subdiagonal (const gsl_matrix_uchar * m,
                                    const size_t k);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_superdiagonal (const gsl_matrix_uchar * m,
                                      const size_t k);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subrow (const gsl_matrix_uchar * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subcolumn (const gsl_matrix_uchar * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array (const unsigned char * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array_with_tda (const unsigned char * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector (const gsl_vector_uchar * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector_with_tda (const gsl_vector_uchar * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_uchar_set_zero (gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_set_identity (gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_set_all (gsl_matrix_uchar * m, unsigned char x);

GSL_FUN int gsl_matrix_uchar_fread (FILE * stream, gsl_matrix_uchar * m) ;
GSL_FUN int gsl_matrix_uchar_fwrite (FILE * stream, const gsl_matrix_uchar * m) ;
GSL_FUN int gsl_matrix_uchar_fscanf (FILE * stream, gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_fprintf (FILE * stream, const gsl_matrix_uchar * m, const char * format);

GSL_FUN int gsl_matrix_uchar_memcpy(gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);
GSL_FUN int gsl_matrix_uchar_swap(gsl_matrix_uchar * m1, gsl_matrix_uchar * m2);

GSL_FUN int gsl_matrix_uchar_swap_rows(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_swap_columns(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_swap_rowcol(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uchar_transpose (gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_transpose_memcpy (gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);

GSL_FUN unsigned char gsl_matrix_uchar_max (const gsl_matrix_uchar * m);
GSL_FUN unsigned char gsl_matrix_uchar_min (const gsl_matrix_uchar * m);
GSL_FUN void gsl_matrix_uchar_minmax (const gsl_matrix_uchar * m, unsigned char * min_out, unsigned char * max_out);

GSL_FUN void gsl_matrix_uchar_max_index (const gsl_matrix_uchar * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_uchar_min_index (const gsl_matrix_uchar * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_uchar_minmax_index (const gsl_matrix_uchar * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_uchar_isnull (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_ispos (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_isneg (const gsl_matrix_uchar * m);
GSL_FUN int gsl_matrix_uchar_isnonneg (const gsl_matrix_uchar * m);

GSL_FUN int gsl_matrix_uchar_add (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_sub (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_mul_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_div_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_FUN int gsl_matrix_uchar_scale (gsl_matrix_uchar * a, const double x);
GSL_FUN int gsl_matrix_uchar_add_constant (gsl_matrix_uchar * a, const double x);
GSL_FUN int gsl_matrix_uchar_add_diagonal (gsl_matrix_uchar * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_uchar_get_row(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t i);
GSL_FUN int gsl_matrix_uchar_get_col(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t j);
GSL_FUN int gsl_matrix_uchar_set_row(gsl_matrix_uchar * m, const size_t i, const gsl_vector_uchar * v);
GSL_FUN int gsl_matrix_uchar_set_col(gsl_matrix_uchar * m, const size_t j, const gsl_vector_uchar * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned char   gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x);
GSL_FUN INLINE_DECL unsigned char * gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned char * gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned char
gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned char *
gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (unsigned char *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned char *
gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const unsigned char *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UCHAR_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
static char help[] =
"Solves a 1D advection plus diffusion problem using FD discretization\n"
"and a structured-grid (DMDA).  Option prefix -b1_.  Equation is\n"
"  - eps u'' + (a(x) u)' = 0\n"
"with a(x)=1, on domain [-1,1], and with Dirichlet boundary conditions\n"
"u(-1) = 1, u(1) = 0.  Default eps=0.01.  The diffusion discretized by\n"
"centered, as usual, but advection is by first-order upwinding, centered,\n"
"or van Leer limiter scheme.  An analytic Jacobian is implemented, except\n"
"for the van Leer limiter.  The limiters in the residual and Jacobian\n"
"evaluations are separately controllable.\n\n";

#include <petsc.h>

typedef enum {NONE, CENTERED, VANLEER} LimiterType;
static const char *LimiterTypes[] = {"none","centered","vanleer",
                                     "LimiterType", "", NULL};

static PetscReal centered(PetscReal theta) {
    return 0.5;
}

static PetscReal vanleer(PetscReal theta) {
    const PetscReal abstheta = PetscAbsReal(theta);
    return 0.5 * (theta + abstheta) / (1.0 + abstheta);
}

static void* limiterptr[] = {NULL, &centered, &vanleer};

typedef struct {
    PetscReal   eps;          // amount of diffusion; require: eps > 0
    PetscReal   (*limiter_fcn)(PetscReal),
                (*jac_limiter_fcn)(PetscReal);
} AdCtx;

static PetscReal u_exact(PetscReal x, AdCtx *usr) {
    return (1.0 - exp((x-1) / usr->eps)) / (1.0 - exp(- 2.0 / usr->eps));
}

static PetscReal wind_a(PetscReal x) {
    return 1.0;
}

extern PetscErrorCode FormUExact(DMDALocalInfo*, AdCtx*, Vec);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, PetscReal*,PetscReal*, AdCtx*);
extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*, PetscReal*, Mat, Mat, AdCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da, da_after;
    SNES           snes;
    Vec            u_initial, u, u_exact;
    PetscReal      hx, err2, errinf;
    DMDALocalInfo  info;
    LimiterType    limiter = NONE, jac_limiter;
    PetscBool      snesfdset, snesfdcolorset;
    AdCtx          user;

    PetscInitialize(&argc,&argv,(char*)0,help);

    user.eps = 0.01;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"b1_",
               "both1d (1D advection-diffusion solver) options",""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-eps","positive diffusion coefficient",
               "both1d.c",user.eps,&(user.eps),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-limiter","flux-limiter type",
               "both1d.c",LimiterTypes,
               (PetscEnum)limiter,(PetscEnum*)&limiter,NULL); CHKERRQ(ierr);
    jac_limiter = limiter;
    ierr = PetscOptionsEnum("-jac_limiter","flux-limiter type used in Jacobian evaluation",
               "both1d.c",LimiterTypes,
               (PetscEnum)jac_limiter,(PetscEnum*)&jac_limiter,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    if (user.eps <= 0.0) {
        SETERRQ1(PETSC_COMM_SELF,2,"eps=%.3f invalid ... eps > 0 required",user.eps);
    }
    user.limiter_fcn = limiterptr[limiter];
    ierr = PetscOptionsHasName(NULL,NULL,"-snes_fd",&snesfdset); CHKERRQ(ierr);
    ierr = PetscOptionsHasName(NULL,NULL,"-snes_fd_color",&snesfdcolorset); CHKERRQ(ierr);
    if (snesfdset || snesfdcolorset) {
        user.jac_limiter_fcn = NULL;
        jac_limiter = 4;   // corresponds to empty string
    } else
        user.jac_limiter_fcn = limiterptr[jac_limiter];

    ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,
                 3, // default to hx=1 grid
                 1,2, // d.o.f., stencil width
                 NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMDASetUniformCoordinates(da,-1.0,1.0,-1.0,1.0,-1.0,1.0); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
            (DMDASNESFunction)FormFunctionLocal,&user);CHKERRQ(ierr);
    ierr = DMDASNESSetJacobianLocal(da,
            (DMDASNESJacobian)FormJacobianLocal,&user); CHKERRQ(ierr);
    ierr = SNESSetApplicationContext(snes,&user); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

    ierr = DMGetGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = VecSet(u_initial,0.0); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);

    ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
         "done on %d point grid (eps = %g, limiter = %s, jac_limiter = %s)\n",
         info.mx,user.eps,LimiterTypes[limiter],LimiterTypes[jac_limiter]); CHKERRQ(ierr);

    ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
    ierr = FormUExact(&info,&user,u_exact); CHKERRQ(ierr);
    ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);    // u <- u + (-1.0) u_exact
    ierr = VecNorm(u,NORM_INFINITY,&errinf); CHKERRQ(ierr);
    ierr = VecNorm(u,NORM_2,&err2); CHKERRQ(ierr);
    hx = 2.0 / (info.mx - 1);
    err2 *= PetscSqrtReal(hx);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
         "numerical error:  |u-uexact|_inf = %.4e,  |u-uexact|_2 = %.4e\n",
         errinf,err2); CHKERRQ(ierr);

    VecDestroy(&u_exact);  SNESDestroy(&snes);
    return PetscFinalize();
}

PetscErrorCode FormUExact(DMDALocalInfo *info, AdCtx *usr, Vec uex) {
    PetscErrorCode  ierr;
    PetscInt   i;
    PetscReal  hx, x, *auex;

    hx = 2.0 / (info->mx - 1);
    ierr = DMDAVecGetArray(info->da, uex, &auex);CHKERRQ(ierr);
    for (i=info->xs; i<info->xs+info->xm; i++) {
        x = -1.0 + i * hx;
        auex[i] = u_exact(x,usr);
    }
    ierr = DMDAVecRestoreArray(info->da, uex, &auex);CHKERRQ(ierr);
    return 0;
}

/* compute residuals:
     F_i = (- eps u'' + (a(x) u)') * hx   at interior points
     F_i = c (u - (b.c.))          at boundary points         */
PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal *au,
                                 PetscReal *aF, AdCtx *usr) {
    const PetscReal  eps = usr->eps,
                     hx = 2.0 / (info->mx - 1),
                     halfx = hx / 2.0,
                     hx2 = hx * hx,
                     scdiag = (2.0 * eps) / hx + 1.0;
    PetscReal        x, uE, uW, uxx, a, u_up, flux, u_dn, u_far, theta;
    PetscInt         i;

    // for each owned cell, non-advective part of residual at cell center
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if (i == 0) {
            aF[i] = scdiag * (au[i] - 1.0);
        } else if (i == info->mx-1) {
            aF[i] = scdiag * (au[i] - 0.0);
        } else {
            uW = (i == 1)          ? 1.0 : au[i-1];
            uE = (i == info->mx-2) ? 0.0 : au[i+1];
            uxx = (uW - 2.0 * au[i] + uE) / hx2;
            aF[i] = - eps * uxx * hx;
        }
    }
    // for each E face of an owned cell, compute flux at the face center
    //     and then add that to the correct residual
    // note -1 start to get W faces of owned cells living on ownership
    //     boundaries
    for (i=info->xs-1; i<info->xs+info->xm; i++) {
        // if cell center is outside [-1,1], or on x=1 boundary, then no need
        // to compute a flux
        if (i < 0 || i == info->mx-1)
            continue;
        // traverse E cell face center points x_{i+1/2} for flux contributions
        // get pth component of wind and first-order upwind flux
        x = -1.0 + i * hx;
        a = wind_a(x + halfx);
        if (a >= 0.0) {
            if (i == 0) {
                u_up = 1.0;
            } else {
                u_up = au[i];
            }
        } else {
            if (i+1 == info->mx-1) {
                u_up = 0.0;
            } else {
                u_up = au[i+1];
            }
        }
        flux = a * u_up;
        if (usr->limiter_fcn != NULL) {
            // flux correction from high-order formula with psi(theta)
            if (a >= 0)
                u_dn = (i+1 < info->mx-1) ? au[i+1] : 0.0;
            else
                u_dn = au[i];
            if (u_dn != u_up) {
                if (a >= 0)
                   u_far = (i-1 > 0) ? au[i-1] : 1.0;
                else
                   u_far = (i+2 < info->mx-1) ? au[i+2] : 0.0;
                theta = (u_up - u_far) / (u_dn - u_up);
                flux += a * (*(usr->limiter_fcn))(theta) * (u_dn - u_up);
            }
        }
        // update non-boundary and owned F_i on both sides of computed flux
        // note aF[] does not have stencil width
        if (i > 0 && i >= info->xs)
            aF[i] += flux;  // flux out of i at E
        if (i+1 < info->mx-1 && i+1 < info->xs + info->xm)
            aF[i+1] -= flux;  // flux into i+1 at W
    }
    return 0;
}

PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscReal *u,
                                 Mat J, Mat P, AdCtx *usr) {
    PetscErrorCode ierr;
    const PetscReal eps = usr->eps,
                    hx = 2.0 / (info->mx - 1),
                    halfx = hx / 2.0,
                    scdiag = (2.0 * eps) / hx + 1.0;
    PetscInt        i, col[3];
    PetscReal       x, aE, aW, v[3];

    if (usr->jac_limiter_fcn == &vanleer) {
        SETERRQ(PETSC_COMM_SELF,1,"Jacobian for vanleer limiter is not implemented");
    }

    ierr = MatZeroEntries(P); CHKERRQ(ierr);

    for (i=info->xs; i<info->xs+info->xm; i++) {
        if (i == 0 || i == info->mx-1) {
            v[0] = scdiag;
            ierr = MatSetValues(P,1,&i,1,&i,v,ADD_VALUES); CHKERRQ(ierr);
        } else {
            // diffusive part
            col[0] = i;
            v[0] = (2.0 * eps) / hx;
            col[1] = i-1;
            v[1] = (i-1 > 0) ? - eps / hx : 0.0;
            col[2] = i+1;
            v[2] = (i+1 < info->mx-1) ? - eps / hx : 0.0;
            ierr = MatSetValues(P,1,&i,3,col,v,ADD_VALUES); CHKERRQ(ierr);
            // advective part: from each adjacent face
            x = -1.0 + i * hx;
            aE = wind_a(x + halfx);
            aW = wind_a(x - halfx);
            if (aE >= 0.0) {
                col[0] = i;
                v[0] = aE;
            } else {
                col[0] = i+1;
                v[0] = (i+1 < info->mx-1) ? aE : 0.0;  // check if i+1 is boundary
            }
            if (aW >= 0.0) {
                col[1] = i-1;
                v[1] = (i-1 > 0) ? - aW : 0.0;  // check if i-1 is boundary
            } else {
                col[1] = i;
                v[1] = - aW;
            }
            ierr = MatSetValues(P,1,&i,2,col,v,ADD_VALUES); CHKERRQ(ierr);
            if (usr->jac_limiter_fcn == &centered) {
                col[0] = i+1;
                col[1] = i;
                if (aE >= 0.0) {
                    v[0] = (i+1 < info->mx-1) ? aE/2.0 : 0.0;  // check if i+1 is boundary
                    v[1] = - aE/2.0;
                } else {
                    v[0] = (i+1 < info->mx-1) ? - aE/2.0 : 0.0;  // check if i+1 is boundary
                    v[1] = aE/2.0;
                }
                ierr = MatSetValues(P,1,&i,2,col,v,ADD_VALUES); CHKERRQ(ierr);
                col[0] = i;
                col[1] = i-1;
                if (aW >= 0.0) {
                    v[0] = - aW/2.0;
                    v[1] = (i-1 > 0) ? aW/2.0 : 0.0;  // check if i-1 is boundary
                } else {
                    v[0] = aW/2.0;
                    v[1] = (i-1 > 0) ? - aW/2.0 : 0.0;  // check if i-1 is boundary
                }
                ierr = MatSetValues(P,1,&i,2,col,v,ADD_VALUES); CHKERRQ(ierr);
            }
        }
    }
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

ext: c


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef type_3f1a28e2_0da6_44b6_84db_29542f5a65c0_h
#define type_3f1a28e2_0da6_44b6_84db_29542f5a65c0_h

#include <utility>
#include <gslib/basetype.h>
#include <gslib/math.h>

__gslib_begin__

template<class _ty, class _protopt>
struct point_t:
    public _protopt
{
    typedef _ty type;
    typedef _protopt proto;
    typedef point_t<_ty, _protopt> myref;

public:
    point_t() { this->x = 0; this->y = 0; }
    point_t(const proto& p): proto(p) {}
    point_t(type a, type b) { this->x = a; this->y = b; }
    void offset(type u, type v) { this->x += u; this->y += v; }
    void offset(const proto& p) { this->x += p.x; this->y += p.y; }
    void set_point(type a, type b) { this->x = a; this->y = b; }
    bool operator == (const myref& that) const { return this->x == that.x && this->y == that.y; }
    bool operator != (const myref& that) const { return this->x != that.x || this->y != that.y; }
};

struct vec2i { int x, y; };
typedef point_t<int, vec2i> point;
typedef point_t<float, vec2> pointf;

template<class _ty, class _ptcls>
struct rect_t
{
    typedef _ty type;
    typedef rect_t<_ty, _ptcls> myref;
    typedef _ptcls point;

public:
    type    left, top, right, bottom;

public:
    rect_t()
    {
        left = 0;
        top = 0;
        right = 0;
        bottom = 0;
    }
    rect_t(type l, type t, type w, type h) { set_rect(l, t, w, h); }
    type width() const { return right - left; }
    type height() const { return bottom - top; }
    void set_rect(type l, type t, type w, type h)
    {
        left = l;
        top = t;
        right = l + w;
        bottom = t + h;
    }
    void set_ltrb(type l, type t, type r, type b)
    {
        left = l;
        top = t;
        right = r;
        bottom = b;
    }
    void set_by_pts(const point& p1, const point& p2)
    {
        left = gs_min(p1.x, p2.x);
        top = gs_min(p1.y, p2.y);
        right = gs_max(p1.x, p2.x);
        bottom = gs_max(p1.y, p2.y);
    }
    bool in_rect(const point& pt) const { return pt.x >= left && pt.x < right && pt.y >= top && pt.y < bottom; }
    void offset(type x, type y) { left += x; right += x; top += y; bottom += y; }
    void deflate(type u, type v);
    void move_to(const point& pt) { move_to(pt.x, pt.y); }
    void move_to(type x, type y);
    bool operator == (const myref& that) const { return left == that.left && right == that.right && top == that.top && bottom == that.bottom; }
    bool operator != (const myref& that) const { return left != that.left || right != that.right || top != that.top || bottom != that.bottom; }
    type area() const { return width() * height(); }
    point center() const
    {
        point c;
        c.x = (left + right) / 2;
        c.y = (top + bottom) / 2;
        return c;
    }
    point top_left() const { return point(left, top); }
    point top_right() const { return point(right, top); }
    point bottom_left() const { return point(left, bottom); }
    point bottom_right() const { return point(right, bottom); }
};

typedef rect_t<int, point> rect;
typedef rect_t<float, pointf> rectf;

inline rectf to_rectf(const rect& rc)
{
    return std::move(rectf((float)rc.left, (float)rc.top, (float)rc.width(), (float)rc.height()));
}

inline rect to_aligned_rect(const rectf& rc)
{
    return std::move(rect(round(rc.left), round(rc.top), round(rc.width()), round(rc.height())));
}

gs_export extern bool intersect_rect(rect& rc, const rect& rc1, const rect& rc2);
gs_export extern bool is_rect_intersected(const rect& rc1, const rect& rc2);
gs_export extern void union_rect(rect& rc, const rect& rc1, const rect& rc2);
gs_export extern bool substract_rect(rect& rc, const rect& rc1, const rect& rc2);
gs_export extern bool intersect_rect(rectf& rc, const rectf& rc1, const rectf& rc2);
gs_export extern bool is_rect_intersected(const rectf& rc1, const rectf& rc2);
gs_export extern bool is_rect_contained(const rect& rc1, const rect& rc2);
gs_export extern bool is_rect_contained(const rectf& rc1, const rectf& rc2);
gs_export extern void union_rect(rectf& rc, const rectf& rc1, const rectf& rc2);
gs_export extern bool substract_rect(rectf& rc, const rectf& rc1, const rectf& rc2);
gs_export extern bool is_line_rect_overlapped(const point& p1, const point& p2, const rect& rc);
gs_export extern bool is_line_rect_overlapped(const pointf& p1, const pointf& p2, const rectf& rc);

__gslib_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <string.h>
#include <math.h>
#include <mpi.h>

#include <gsl/gsl_rng.h>

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>
#include "pmpfft.h"

/* The following functions fill the gaussian field*/
static void
pmic_fill_gaussian_gadget(PM * pm, FastPMFloat * delta_k, int seed);
static void
pmic_fill_gaussian_fast(PM * pm, FastPMFloat * delta_k, int seed);
static void
pmic_fill_gaussian_slow(PM * pm, FastPMFloat * delta_k, int seed);

void
fastpm_ic_fill_gaussiank(PM * pm, FastPMFloat * delta_k, int seed, enum FastPMFillDeltaKScheme scheme)
{

    /* clear the memory to avoid any modes that we forget to set. */
    memset(delta_k, 0, pm_allocsize(pm) * sizeof(delta_k[0]));

    switch(scheme) {
        case FASTPM_DELTAK_GADGET:
            pmic_fill_gaussian_gadget(pm, delta_k, seed);
            break;
        case FASTPM_DELTAK_FAST:
            pmic_fill_gaussian_fast(pm, delta_k, seed);
            break;
        case FASTPM_DELTAK_SLOW:
            pmic_fill_gaussian_slow(pm, delta_k, seed);
            break;
        default:
            pmic_fill_gaussian_gadget(pm, delta_k, seed);
            break;
    }
}

struct PofK {
    fastpm_fkfunc func;
    void * data;
    double Volume;
} ;

static double _powerspec_to_transfer(double k, struct PofK * pk)
{
    double f = sqrt(pk->func(k, pk->data));
    f *= sqrt(1.0 / pk->Volume);
    return f;
}

void
fastpm_ic_induce_correlation(PM * pm, FastPMFloat * delta_k, fastpm_fkfunc pkfunc, void * data)
{
    struct PofK pk;
    pk.func = pkfunc;
    pk.data = data;
    pk.Volume = pm->Volume;

    fastpm_apply_any_transfer(pm, delta_k, delta_k, (fastpm_fkfunc) _powerspec_to_transfer, &pk);
}

void
fastpm_ic_remove_variance(PM * pm, FastPMFloat * delta_k)
{
#pragma omp parallel
    {
        PMKIter kiter;

        for(pm_kiter_init(pm, &kiter);
            !pm_kiter_stop(&kiter);
            pm_kiter_next(&kiter)) {
            double k2 = 0;
            int d;
            for(d = 0; d < 3; d++) {
                k2 += kiter.kk[d][kiter.iabs[d]];
            }

            /* https://en.wikipedia.org/wiki/Atan2 */
            double a = delta_k[kiter.ind];
            double b = delta_k[kiter.ind + 1];

            if(a == 0 && b == 0)   {
                delta_k[kiter.ind + 0] = 0;
                delta_k[kiter.ind + 1] = 0;
            } else {
                double phase = atan2(b, a);
                delta_k[kiter.ind + 0] = cos(phase);
                delta_k[kiter.ind + 1] = sin(phase);
            }

        }

    }
}

static inline void
SETSEED(PM * pm, unsigned int * table[2][2], int i, int j, gsl_rng * rng)
{
    unsigned int seed = 0x7fffffff * gsl_rng_uniform(rng);

    int ii[2] = {i, (pm->Nmesh[0] - i) % pm->Nmesh[0]};
    int jj[2] = {j, (pm->Nmesh[1] - j) % pm->Nmesh[1]};
    int d1, d2;
    for(d1 = 0; d1 < 2; d1++) {
        ii[d1] -= pm->ORegion.start[0];
        jj[d1] -= pm->ORegion.start[1];
    }
    for(d1 = 0; d1 < 2; d1++)
    for(d2 = 0; d2 < 2; d2++) {
        if( ii[d1] >= 0 &&
            ii[d1] < pm->ORegion.size[0] &&
            jj[d2] >= 0 &&
            jj[d2] < pm->ORegion.size[1]
        ) {
            table[d1][d2][ii[d1] * pm->ORegion.size[1] + jj[d2]] = seed;
        }
    }
}
static inline unsigned int
GETSEED(PM * pm, unsigned int * table[2][2], int i, int j, int d1, int d2)
{
    i -= pm->ORegion.start[0];
    j -= pm->ORegion.start[1];
    if(i < 0) abort();
    if(j < 0) abort();
    if(i >= pm->ORegion.size[0]) abort();
    if(j >= pm->ORegion.size[1]) abort();
    return table[d1][d2][i * pm->ORegion.size[1] + j];
}

static void
SAMPLE(gsl_rng * rng, double * ampl, double * phase)
{
    *phase = gsl_rng_uniform(rng) * 2 * M_PI;
    *ampl = 0;
    do *ampl = gsl_rng_uniform(rng); while(*ampl == 0);
}

static void
pmic_fill_gaussian_gadget(PM * pm, FastPMFloat * delta_k, int seed)
{
    /* Fill delta_k with gadget scheme */
    int d;
    int i, j, k;

    memset(delta_k, 0, sizeof(delta_k[0]) * pm->allocsize);

    gsl_rng * rng = gsl_rng_alloc(gsl_rng_ranlxd1);
    gsl_rng_set(rng, seed);

    unsigned int * seedtable[2][2];
    for(i = 0; i < 2; i ++)
    for(j = 0; j < 2; j ++) {
            seedtable[i][j] = calloc(pm->ORegion.size[0] * pm->ORegion.size[1], sizeof(int));
    }

    for(i = 0; i < pm->Nmesh[0] / 2; i++) {
        for(j = 0; j < i; j++) SETSEED(pm, seedtable, i, j, rng);
        for(j = 0; j < i + 1; j++) SETSEED(pm, seedtable, j, i, rng);
        for(j = 0; j < i; j++) SETSEED(pm, seedtable, pm->Nmesh[0] - 1 - i, j, rng);
        for(j = 0; j < i + 1; j++) SETSEED(pm, seedtable, pm->Nmesh[1] - 1 - j, i, rng);
        for(j = 0; j < i; j++) SETSEED(pm, seedtable, i, pm->Nmesh[1] - 1 - j, rng);
        for(j = 0; j < i + 1; j++) SETSEED(pm, seedtable, j, pm->Nmesh[0] - 1 - i, rng);
        for(j = 0; j < i; j++) SETSEED(pm, seedtable, pm->Nmesh[0] - 1 - i, pm->Nmesh[1] - 1 - j, rng);
        for(j = 0; j < i + 1; j++) SETSEED(pm, seedtable, pm->Nmesh[1] - 1 - j, pm->Nmesh[0] - 1 - i, rng);
    }
    gsl_rng_free(rng);

    ptrdiff_t irel[3];
    for(i = pm->ORegion.start[0];
        i < pm->ORegion.start[0] + pm->ORegion.size[0];
        i ++) {

        gsl_rng * lower_rng = gsl_rng_alloc(gsl_rng_ranlxd1);
        gsl_rng * this_rng = gsl_rng_alloc(gsl_rng_ranlxd1);

        int ci = pm->Nmesh[0] - i;
        if(ci >= pm->Nmesh[0]) ci -= pm->Nmesh[0];

        for(j = pm->ORegion.start[1];
            j < pm->ORegion.start[1] + pm->ORegion.size[1];
            j ++) {
            /* always pull the gaussian from the lower quadrant plane for k = 0
             * plane*/
            /* always pull the whitenoise from the lower quadrant plane for k = 0
             * plane and k == Nmesh / 2 plane*/
            int d1 = 0, d2 = 0;
            int cj = pm->Nmesh[1] - j;
            if(cj >= pm->Nmesh[1]) cj -= pm->Nmesh[1];

            /* d1, d2 points to the conjugate quandrant */
            if( (ci == i && cj < j)
             || (ci < i && cj != j)
             || (ci < i && cj == j)) {
                d1 = 1;
                d2 = 1;
            }

            unsigned int seed_conj, seed_this;
            /* the lower quadrant generator */
            seed_conj = GETSEED(pm, seedtable, i, j, d1, d2);
            gsl_rng_set(lower_rng, seed_conj);

            seed_this = GETSEED(pm, seedtable, i, j, 0, 0);
            gsl_rng_set(this_rng, seed_this);

            for(k = 0; k <= pm->Nmesh[2] / 2; k ++) {
                int use_conj = (d1 != 0 || d2 != 0) && (k == 0 || k == pm->Nmesh[2] / 2);

                double ampl, phase;
                if(use_conj) {
                    /* on k = 0 and Nmesh/2 plane, we use the lower quadrant generator,
                     * then hermit transform the result if it is nessessary */
                    SAMPLE(this_rng, &ampl, &phase);
                    SAMPLE(lower_rng, &ampl, &phase);
                } else {
                    SAMPLE(lower_rng, &ampl, &phase);
                    SAMPLE(this_rng, &ampl, &phase);
                }

                ptrdiff_t iabs[3] = {i, j, k};
                ptrdiff_t ip = 0;
                for(d = 0; d < 3; d ++) {
                    irel[d] = iabs[d] - pm->ORegion.start[d];
                    ip += pm->ORegion.strides[d] * irel[d];
                }

                if(irel[2] < 0) continue;
                if(irel[2] >= pm->ORegion.size[2]) continue;

                /* we want two numbers that are of std ~ 1/sqrt(2) */
                ampl = sqrt(- log(ampl));

                (delta_k + 2 * ip)[0] = ampl * cos(phase);
                (delta_k + 2 * ip)[1] = ampl * sin(phase);

                if(use_conj) {
                    (delta_k + 2 * ip)[1] *= -1;
                }

                if((pm->Nmesh[0] - iabs[0]) % pm->Nmesh[0] == iabs[0] &&
                   (pm->Nmesh[1] - iabs[1]) % pm->Nmesh[1] == iabs[1] &&
                   (pm->Nmesh[2] - iabs[2]) % pm->Nmesh[2] == iabs[2]) {
                    /* The mode is self conjuguate, thus imaginary mode must be zero */
                    (delta_k + 2 * ip)[1] = 0;
                    (delta_k + 2 * ip)[0] = ampl * cos(phase);
                }

                if(iabs[0] == 0 && iabs[1] == 0 && iabs[2] == 0) {
                    /* the mean is zero */
                    (delta_k + 2 * ip)[0] = 0;
                    (delta_k + 2 * ip)[1] = 0;
                }
            }
        }
        gsl_rng_free(lower_rng);
        gsl_rng_free(this_rng);
    }
    for(i = 0; i < 2; i ++)
    for(j = 0; j < 2; j ++) {
        free(seedtable[i][j]);
    }
/*
    char * fn[1000];
    sprintf(fn, "canvas.dump.f4.%d", pm->ThisTask);
    fwrite(pm->canvas, sizeof(pm->canvas[0]), pm->ORegion.total * 2, fopen(fn, "w"));
*/
}

static void
pmic_fill_gaussian_fast(PM * pm, FastPMFloat * delta_k, int seed)
{
    ptrdiff_t ind;
    int d;

    gsl_rng* random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);

    /* set uncorrelated seeds */
    gsl_rng_set(random_generator, seed);
    for(d = 0; d < pm->ThisTask * 8; d++) {
        seed = 0x7fffffff * gsl_rng_uniform(random_generator);
    }

    gsl_rng_set(random_generator, seed);

    FastPMFloat * g_x = pm_alloc(pm);

    for(ind = 0; ind < pm->IRegion.total; ind += 2) {
        double phase = gsl_rng_uniform(random_generator) * 2 * M_PI;
        double ampl;
        do
            ampl = gsl_rng_uniform(random_generator);
        while(ampl == 0.0);

        /* we need two gaussians of std=1.0 in real space (see footnote 1) */
        ampl = sqrt(-2.0 * log(ampl));
        /* r2c will reduce the variance, so we compensate here. */
        ampl *= sqrt(pm_norm(pm));

        g_x[ind] = ampl * sin(phase);
        g_x[ind + 1] = ampl * cos(phase);
    }
    pm_r2c(pm, g_x, delta_k);
    pm_free(pm, g_x);
}

static void
pmic_fill_gaussian_slow(PM * pm, FastPMFloat * delta_k, int seed)
{
    ptrdiff_t i[3] = {0};
    int d;
    gsl_rng* random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);

    gsl_rng_set(random_generator, seed);

    FastPMFloat * g_x = pm_alloc(pm);

    for(i[0] = 0; i[0] < pm->Nmesh[0]; i[0]++)
    for(i[1] = 0; i[1] < pm->Nmesh[1]; i[1]++)
    for(i[2] = 0; i[2] < pm->Nmesh[2]; i[2]++) {
        double phase = gsl_rng_uniform(random_generator) * 2 * M_PI;
        double ampl;
        do
            ampl = gsl_rng_uniform(random_generator);
        while(ampl == 0.0);
        ptrdiff_t ii[3];
        ptrdiff_t ind = 0;
        for(d = 0; d < 3; d ++) {
            if(i[d] < pm->IRegion.start[d]) goto next;
            if(i[d] >= pm->IRegion.start[d] + pm->IRegion.size[d]) goto next;
            ii[d] = i[d] - pm->IRegion.start[d];
            ind += ii[d] * pm->IRegion.strides[d];
        }
        /* we need two gaussians of std=1.0 in real space */
        ampl = sqrt(-2.0 * log(ampl));

        /* r2c will reduce the variance, so we compensate here. */
        ampl *= sqrt(pm_norm(pm));

        g_x[ind] = ampl * sin(phase);
        next:
        continue;
    }
    pm_r2c(pm, g_x, delta_k);
    pm_free(pm, g_x);
    gsl_rng_free(random_generator);
}

/* Footnotes */

/* 1):
 * We want delta(k) = delta_real + I delta_imag, where delta_real and
 * delta_imag are Gaussian random variables with variance given by
 * power spectrum, \sigma^2=P(k). We can obtain this equivalently as
 *
 *   delta(k) = A exp(i phase),
 *
 * where the phase is random (i.e. sampled from a uniform distribution)
 * and the amplitude A follows a Rayleigh distribution (see
 * https://en.wikipedia.org/wiki/Rayleigh_distribution). To sample from
 * Rayleigh distribution, use inverse transform sampling
 * (see https://en.wikipedia.org/wiki/Inverse_transform_sampling), i.e.
 * start from uniform random variable in [0,1] and then apply inverse of CDF
 * of Rayleigh distribution. From F(A)=CDF(A)=1-e^{-A^2/(2\sigma^2)} we get
 * A = \sigma \sqrt{-2 ln(1-CDF)}. So if x is uniform random number in [0,1], then
 * A = \sigma \sqrt(-2 ln(x)) follows Rayleigh distribution as desired.
 * Here we used x instead of 1-x because this does not make a difference for a
 * uniform random number in [0,1]. In the code below, we start with \sigma=1 and
 * multiply by sqrt(P(k)) later.
 */

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include "gsl/gsl_roots.h"
#include "psrsalsa.h"
int minimize_1D_double_refine_borders(int findroot, double (*funk)(double, void *), void *params, int gridsearch, int investigateLocalMinima, double *x_lower, double *x_upper, int debug_verbose)
{
  double x, y, ymin[3], ymin_new[3], dy[3], ymax, x_lower_new, x_upper_new, sign_old, gridsearch_margin, imin_new[3], i_originalgrid, imin1, imin2;
  int imin[3], unset[3], i, minima;
  unset[0] = unset[1] = unset[2] = 1;
  for(i = 0; i < gridsearch; i++) {
    x = *x_lower + i*(*x_upper - *x_lower)/(double)(gridsearch-1);
    y = funk(x, params);
    if(debug_verbose) {
      fflush(stdout); fprintf(stderr, "x=%f y=%f\n", x, y);
    }
    if(findroot) {
      if(i == 0) {
 sign_old = y;
 imin[0] = 0;
      }else {
 if(y < 0 && sign_old > 0) {
   imin[0] = i;
 }else if(y > 0 && sign_old < 0) {
   imin[0] = i;
 }
 if(imin[0] != 0)
   break;
      }
    }else {
      if(unset[0]) {
 imin[0] = i;
 ymin[0] = y;
 unset[0] = 0;
      }else if(unset[1]) {
 imin[1] = i;
 ymin[1] = y;
 unset[1] = 0;
      }else if(unset[2]) {
 imin[2] = i;
 ymin[2] = y;
 unset[2] = 0;
      }else {
 for(minima = 0; minima < 3; minima++) {
   dy[minima] = ymin[minima] - y;
 }
 if(dy[0] > dy[1] && dy[0] > dy[2]) {
   if(dy[0] > 0) {
     imin[0] = i;
     ymin[0] = y;
   }
 }else if(dy[1] > dy[0] && dy[1] > dy[2]) {
   if(dy[1] > 0) {
     imin[1] = i;
     ymin[1] = y;
   }
 }else {
   if(dy[2] > 0) {
     imin[2] = i;
     ymin[2] = y;
   }
 }
      }
      if(i == 0 || y > ymax) {
 ymax = y;
      }
    }
  }
  if(findroot == 0) {
    imin1 = imin[0];
    imin2 = imin[0];
    if(imin[1] < imin1)
      imin1 = imin[1];
    if(imin[2] < imin1)
      imin1 = imin[2];
    if(imin[1] > imin2)
      imin2 = imin[1];
    if(imin[2] > imin2)
      imin2 = imin[2];
    if(imin2-imin1 != 2 && debug_verbose) {
      fflush(stdout); fprintf(stderr, "Found more than one possible minima (range = %.2f .. %.2f)\n", imin1, imin2);
    }
    if(imin2-imin1 > 0.3*gridsearch) {
      if(investigateLocalMinima == 0) {
 imin1 = imin[0];
 imin2 = imin[0];
 if(ymin[1] < ymin[0]) {
   imin1 = imin[1];
   imin2 = imin[1];
 }else if(ymin[2] < ymin[0] && ymin[2] < ymin[1]) {
   imin1 = imin[2];
   imin2 = imin[2];
 }
      }else {
 if(debug_verbose) {
   fflush(stdout); fprintf(stderr, "Refining gridsearch around local minima\n");
 }
 unset[0] = unset[1] = unset[2] = 1;
 for(minima = 0; minima < 3; minima++) {
   for(x = *x_lower + (imin[minima]-1.5)*(*x_upper - *x_lower)/(double)(gridsearch-1); x < *x_lower + (imin[minima]+1.5)*(*x_upper - *x_lower)/(double)(gridsearch-1); x += 0.09*(*x_upper - *x_lower)/(double)(gridsearch-1)) {
     y = funk(x, params);
     if(debug_verbose) {
       fflush(stdout); fprintf(stderr, "x=%f y=%f\n", x, y);
     }
     i_originalgrid = (x - *x_lower)*(double)(gridsearch-1)/(*x_upper - *x_lower);
     if(unset[0] || y < ymin_new[0]) {
       imin_new[0] = i_originalgrid;
       ymin_new[0] = y;
       unset[0] = 0;
     }else if(unset[1] || y < ymin_new[1]) {
       imin_new[1] = i_originalgrid;
       ymin_new[1] = y;
       unset[1] = 0;
     }else if(unset[2] || y < ymin_new[2]) {
       imin_new[2] = i_originalgrid;
       ymin_new[2] = y;
       unset[2] = 0;
     }
   }
 }
 imin1 = imin_new[0];
 imin2 = imin_new[0];
 if(imin_new[1] < imin1)
   imin1 = imin_new[1];
 if(imin_new[2] < imin1)
   imin1 = imin_new[2];
 if(imin_new[1] > imin2)
   imin2 = imin_new[1];
 if(imin_new[2] > imin2)
   imin2 = imin_new[2];
 if(imin2-imin1 != 2 && debug_verbose) {
   fflush(stdout); fprintf(stderr, "Found more than one possible minima (range = %.2f .. %.2f)\n", imin1, imin2);
 }
 if(imin2-imin1 > 0.3*gridsearch) {
   fflush(stdout); fprintf(stderr, "minimize_1D_double_refine_borders: Refining gridsearch around local minima failed\n");
   exit(0);
 }
      }
    }
  }
  gridsearch_margin = 0.1*(double)gridsearch;
  if(gridsearch_margin < 1)
    gridsearch_margin = 1;
  if(findroot) {
    if(imin[0] == 0) {
      return 2;
    }
    x_lower_new = *x_lower + (imin[0]-gridsearch_margin)*(*x_upper - *x_lower)/(double)(gridsearch-1);
    x_upper_new = *x_lower + (imin[0]+gridsearch_margin-1)*(*x_upper - *x_lower)/(double)(gridsearch-1);
    *x_lower = x_lower_new;
    *x_upper = x_upper_new;
  }else {
    x_lower_new = *x_lower + (imin1-gridsearch_margin)*(*x_upper - *x_lower)/(double)(gridsearch-1);
    x_upper_new = *x_lower + (imin2+gridsearch_margin)*(*x_upper - *x_lower)/(double)(gridsearch-1);
    *x_lower = x_lower_new;
    *x_upper = x_upper_new;
  }
  return 0;
}
int minimize_1D_double(int findroot, double (*funk)(double, void *), void *params, double x_lower, double x_upper, int gridsearch, int investigateLocalMinima, int nested, double *x_minimum, int max_iter, double epsabs, double epsrel, int verbose, int debug_verbose)
{
  const gsl_min_fminimizer_type *T_minimizer;
  gsl_min_fminimizer *s_minimizer;
  const gsl_root_fsolver_type *T_root;
  gsl_root_fsolver *s_root;
  int status;
  int iter, i, ret, nest;
  gsl_function F;
  F.function = funk;
  F.params = params;
  iter = 0;
#if GSL_VERSION_NUMBER < 102
  printerror(0, "ERROR minimize_1D_double: Not supported for GSL < 1.2");
  exit(0);
#endif
  if(verbose) {
    for(i = 0; i < verbose - 1; i++)
      printf(" ");
    if(findroot == 0)
      printf("Minimising function between %f and %f\n", x_lower, x_upper);
    else
      printf("Finding root of function between %f and %f\n", x_lower, x_upper);
  }
  if(gridsearch > 1) {
    if(verbose) {
      for(i = 0; i < verbose - 1; i++)
 printf(" ");
      printf("  Refining boundaries with a %d point grid search\n", gridsearch);
    }
    for(nest = 0; nest < 1+nested; nest++) {
      ret = minimize_1D_double_refine_borders(findroot, funk, params, gridsearch, investigateLocalMinima, &x_lower, &x_upper, debug_verbose);
      if(ret != 0) {
 if(ret == 1) {
   if(verbose) {
     for(i = 0; i < verbose - 1; i++)
       printf(" ");
     printf("  Maximum itterations (%d) exceeded while refining borders\n", max_iter);
   }
 }
 return ret;
      }
      if(verbose) {
 for(i = 0; i < verbose - 1; i++)
   printf(" ");
 printf("  Refined boundaries are %f and %f\n", x_lower, x_upper);
      }
    }
  }
  *x_minimum = 0.5*(x_upper+x_lower);
  if(findroot == 0) {
    T_minimizer = gsl_min_fminimizer_brent;
    s_minimizer = gsl_min_fminimizer_alloc(T_minimizer);
    fflush(stdout);
    gsl_min_fminimizer_set(s_minimizer, &F, *x_minimum, x_lower, x_upper);
  }else {
    T_root = gsl_root_fsolver_brent;
    s_root = gsl_root_fsolver_alloc(T_root);
    if(gridsearch == 0) {
      double val1, val2;
      val1 = funk(x_lower, params);
      val2 = funk(x_upper, params);
      if((val1 > 0 && val2 > 0) || (val1 < 0 && val2 < 0)) {
 return 3;
      }
    }
    gsl_root_fsolver_set(s_root, &F, x_lower, x_upper);
  }
  do {
    iter++;
    if(findroot == 0) {
      status = gsl_min_fminimizer_iterate(s_minimizer);
#if GSL_VERSION_NUMBER >= 102
      *x_minimum = gsl_min_fminimizer_x_minimum (s_minimizer);
#endif
      x_lower = gsl_min_fminimizer_x_lower(s_minimizer);
      x_upper = gsl_min_fminimizer_x_upper(s_minimizer);
    }else {
      status = gsl_root_fsolver_iterate(s_root);
      *x_minimum = gsl_root_fsolver_root(s_root);
      x_lower = gsl_root_fsolver_x_lower(s_root);
      x_upper = gsl_root_fsolver_x_upper(s_root);
    }
    status = gsl_min_test_interval (x_lower, x_upper, epsabs, epsrel);
  }while(status == GSL_CONTINUE && iter < max_iter);
  if(findroot)
    gsl_root_fsolver_free(s_root);
  else
    gsl_min_fminimizer_free(s_minimizer);
  if(status == GSL_SUCCESS) {
    if(verbose) {
      for(i = 0; i < verbose - 1; i++)
 printf(" ");
      if(findroot)
 printf("  Root found at %f with error %f\n", *x_minimum, fabs(x_upper - x_lower));
      else
 printf("  Minimum found at %f with error %f\n", *x_minimum, fabs(x_upper - x_lower));
    }
    return 0;
  }else if(iter == max_iter) {
    if(verbose) {
      for(i = 0; i < verbose - 1; i++)
 printf(" ");
      printf("  Maximum itterations (%d) exceeded. Value confined to [%f, %f]\n", max_iter, x_lower, x_upper);
    }
    return 1;
  }else {
    if(verbose) {
      for(i = 0; i < verbose - 1; i++)
 printf(" ");
      printf("  Unknown error during minimization\n");
    }
    return 5;
  }
  return status;
}
double (*internal_funk_provided_by_user)(double *, void *);
int internal_paramnr;
double *internal_xminimum;
double internal_desired_chi2;
double internal_find_1D_error_funk(double x, void *params)
{
  double chi2;
  internal_xminimum[internal_paramnr] = x;
  chi2 = internal_funk_provided_by_user(internal_xminimum, params);
  chi2 -= internal_desired_chi2;
  return chi2;
}
int find_1D_error(double (*funk)(double *, void *), double *xminimum, int paramnr, int nrparameters, double dx, double dxmax, void *params, double sigma, double chi2min, int max_itr, double epsabs, double epsrel, double *errorbar, int verbose)
{
  int ittr, i, ret, debug_verbose;
  double x_lower, x_upper, diff, xval, *xminimum_fiddle;
  debug_verbose = 0;
  dx = fabs(dx);
  if(verbose) {
    for(i = 0; i < verbose - 1; i++)
      printf(" ");
    printf("Finding error at value %f with stepsize %f\n", xminimum[paramnr], dx);
  }
  xminimum_fiddle = malloc(nrparameters*sizeof(double));
  if(xminimum_fiddle == NULL) {
    fflush(stdout);
    fprintf(stderr, "ERROR find_1D_error: Memory allocation error\n");
    return 4;
  }
  memcpy(xminimum_fiddle, xminimum, nrparameters*sizeof(double));
  internal_paramnr = paramnr;
  internal_xminimum = xminimum_fiddle;
  internal_funk_provided_by_user = funk;
  internal_desired_chi2 = chi2min*(1+fabs(sigma));
  if(verbose) {
    for(i = 0; i < verbose - 1; i++)
      printf(" ");
    printf("  Mimimum chi2 = %f, so %f sigma point corresponds to %f\n", chi2min, fabs(sigma), internal_desired_chi2);
  }
  x_lower = xminimum[paramnr];
  x_upper = xminimum[paramnr];
  ittr = 0;
  diff = internal_find_1D_error_funk(x_upper, params);
  if(diff > 0) {
    fflush(stdout);
    fprintf(stderr, "ERROR find_1D_error: Function called with initial parameters outside specified sigma limit: chi2 = %f higher than sigma border (%f).\n", diff, internal_desired_chi2);
    exit(0);
  }
  do {
    if(sigma >= 0) {
      x_upper += dx;
      if(dxmax >= 0) {
 if(fabs(x_upper - xminimum[paramnr]) > dxmax) {
   *errorbar = dxmax;
   free(xminimum_fiddle);
   return 0;
 }
      }
      diff = internal_find_1D_error_funk(x_upper, params);
    }else {
      x_lower -= dx;
      if(dxmax >= 0) {
 if(fabs(x_lower - xminimum[paramnr]) > dxmax) {
   *errorbar = dxmax;
   free(xminimum_fiddle);
   return 0;
 }
      }
      diff = internal_find_1D_error_funk(x_lower, params);
    }
    ittr++;
    if(ittr == max_itr) {
      free(xminimum_fiddle);
      return 1;
    }
  }while(diff < 0);
  if(verbose) {
    for(i = 0; i < verbose - 1; i++)
      printf(" ");
    printf("  Found brackets: [%f %f]\n", x_lower, x_upper);
    verbose += 2;
  }
  ret = minimize_1D_double(1, internal_find_1D_error_funk, params, x_lower, x_upper, 0, 0, 0, &xval, max_itr, epsabs, epsrel, verbose, debug_verbose);
  *errorbar = fabs(xval - xminimum[paramnr]);
  if(verbose) {
    verbose -= 2;
    for(i = 0; i < verbose - 1; i++)
      printf(" ");
    if(ret == 0)
      printf("  Found errorbar: %f\n", *errorbar);
    else if(ret == 1)
      printf("  Maximum nr of itterations reached, did not converge fully\n");
    else if(ret == 2)
      printf("  Did not find root in specified range\n");
    else if(ret == 3)
      printf("  Lower and upper limit do not bracket a root\n");
    else {
      ret = 5;
      printf("  Other unspecified error\n");
    }
  }
  free(xminimum_fiddle);
  return ret;
}

max_stars_repo_name: David-McKenna/psrsalsa

max_issues_repo_name: David-McKenna/psrsalsa

max_forks_repo_name: David-McKenna/psrsalsa


text:
/*
 * SPDX-License-Identifier:  BSD-3-Clause
 *
 * example_08-StructOfArrays-CellLinkedList-OuterOmp.c :
 *      Example of SPH Density Calculation using
 *      fast neighbor search the main density loop via
 *      Cell Linked List method, Struct of Arrays (SoA)
 *      data layout, OpenMP parallelization at the
 *      cell level, SIMD directives in the kernel
 *      and in the inner-most loop.
 *
 * (C) Copyright 2021 José Hugo Elsas
 * Author: José Hugo Elsas <jhelsas@gmail.com>
 *
 * Command Line Options:
 *   -runs  <int>   : Set the number of repetitions (runs) for
 *                      calculating the density. The value of
 *                      the density is based on the last
 *                      iteration.
 *                    Default value: 1
 *   -run_seed <int>: Flag to set an alternative seed use for
 *                      for the PRNG. Instead of feeding seed
 *                      to the PRNG directly, it feeds
 *                      seed + iteration, as to generate different
 *                      configurations for each iteration.
 *                    Default value: 0 - (possible 0/1)
 *   -seed     <int>: Set the seed to use for the SPH particles
 *                      uniform position generation in the box
 *                    Default value: 123123123
 *
 *   -N        <int>: Set the number of SPH particles to be used
 *                    Default value: 1e5 = 100,000
 *   -h      <float>: Set the value of the smoothing kernel
 *                      parameter h, which corresponds to half
 *                      of the support of the kernel.
 *                    Default value: 0.05
 *
 *   -Nx       <int>: Set the number of Cells in the X direction
 *                    Default value: 10
 *   -Ny       <int>: Set the number of Cells in the Y direction
 *                    Default value: 10
 *   -Nz       <int>: Set the number of Cells in the Z direction
 *                    Default value: 10
 *
 *   -Xmin   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Y direction for
 *                    the Cell Linked List box
 *                      Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *
 *   -Xmax   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Ymax   <float>: Set the lower bound in the Y direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Zmax   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 */

#include <math.h>
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <stdbool.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_utils.h"

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#define COMPUTE_BLOCKS 5

int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr, linkedListBox *box, SPHparticle *lsph, double *times);

int compute_density_3d_outerOmp(int N, double h, SPHparticle *lsph, linkedListBox *box);

int compute_density_3d_chunk_noomp(int64_t node_begin, int64_t node_end,
                                   int64_t nb_begin, int64_t nb_end,double h,
                                   double* restrict x, double* restrict y,
                                   double* restrict z, double* restrict nu,
                                   double* restrict rho);

int count_box_pairs(linkedListBox *box);

int setup_box_pairs(linkedListBox *box,
                    int64_t *node_begin,int64_t *node_end,
                    int64_t *nb_begin,int64_t *nb_end);

double w_bspline_3d_constant(double h);

#pragma omp declare simd
double w_bspline_3d_simd(double q);

int main(int argc, char **argv){
  bool run_seed = false;       // By default the behavior is is to use the same seed
  int runs = 1,err;            // it only runs once
  long int seed = 123123123;   // The default seed is 123123123
  int64_t N = 100000;          // The default number of particles is N = 1e5 = 100,000
  double h=0.05;               // The default kernel smoothing length is h = 0.05
  linkedListBox *box;          // Uninitialized Box containing the cells for the cell linked list method
  SPHparticle *lsph;           // Uninitialized array of SPH particles

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox)); // Create a box representing the entire 3d domain

  // allow for command line customization of the run
  arg_parse(argc,argv,&N,&h,&seed,&runs,&run_seed,box);  // Parse the command line options
                                                         // line arguments and override default values

  err = SPHparticle_SoA_malloc(N,&lsph);
  if(err)
    fprintf(stderr,"error in SPHparticle_SoA_malloc\n");

  void *swap_arr = malloc(N*sizeof(double));
  double times[runs*COMPUTE_BLOCKS];

  for(int run=0;run<runs;run+=1)
    main_loop(run,run_seed,N,h,seed,swap_arr,box,lsph,times);

  bool is_cll = true;
  const char *prefix = "ex08,cll,SoA,outer,simd";
  print_time_stats(prefix,is_cll,N,h,seed,runs,lsph,box,times);
  print_sph_particles_density(prefix,is_cll,N,h,seed,runs,lsph,box);

  SPHparticleSOA_safe_free(N,&lsph);
  safe_free_box(box);
  free(swap_arr);

  return 0;
}

/*
 *  Function main_loop:
 *    Runs the main loop of the program, including the particle array generation,
 *    density calculation and the timings annotations.
 *
 *    Arguments:
 *       run <int>            : index (or value) or the present iteration
 *       run_seed <bool>      : boolean defining whether to use run index for seed or not
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       seed <long int>      : seed for GSL PRNG generator to generate particle positions
 *       box  <linkedListBox> : Box of linked list cells, encapsulating the 3d domain
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *       times <double>       : Array to store the computation timings to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 *       times <double>       : Times is updated by reference
 */
int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr, linkedListBox *box, SPHparticle *lsph, double *times)
{
  int err;

  if(run_seed)
    err = gen_unif_rdn_pos_box(N,seed+run,box,lsph);
  else
    err = gen_unif_rdn_pos_box(N,seed,box,lsph);

  if(err)
    fprintf(stderr,"error in gen_unif_rdn_pos\n");

  // ------------------------------------------------------ //

  double t0,t1,t2,t3,t4,t5;

  t0 = omp_get_wtime();

  err = compute_hash_MC3D(N,lsph,box);                      // Compute Morton Z 3D hash based on the
  if(err)                                                   // cell index for each of the X, Y and Z
    fprintf(stderr,"error in compute_hash_MC3D\n");         // directions, in which a given particle reside

  t1 = omp_get_wtime();

  qsort(lsph->hash,N,2*sizeof(int64_t),compare_int64_t);    // Sort the Particle Hash Hashes, getting the shuffled
                                                            // index necessary to re-shuffle the remaining arrays
  t2 = omp_get_wtime();

  err = reorder_lsph_SoA(N,lsph,swap_arr);                  // Reorder all arrays according to the sorted hash,
  if(err)                                                   // As to have a quick way to retrieve a cell
    fprintf(stderr,"error in reorder_lsph_SoA\n");          // given its hash.

  t3 = omp_get_wtime();

  err = setup_interval_hashtables(N,lsph,box);              // Annotate the begining and end of each cell
  if(err)                                                   // on the cell linked list method for fast
    fprintf(stderr,"error in setup_interval_hashtables\n"); // neighbor search

  t4 = omp_get_wtime();

  err = compute_density_3d_outerOmp(N,h,lsph,box);          // Compute the density of the particles based
  if(err)                                                   // on the cell linked list method for fast
    fprintf(stderr,"error in compute_density\n");           // neighbor search

  // ------------------------------------------------------ //

  t5 = omp_get_wtime();

  times[COMPUTE_BLOCKS*run+0] = t1-t0;                      // Time for compute morton Z 3d hash
  times[COMPUTE_BLOCKS*run+1] = t2-t1;                      // Time for sorting the particles' hashes
  times[COMPUTE_BLOCKS*run+2] = t3-t2;                      // Time for reordering all other arrays accordingly
  times[COMPUTE_BLOCKS*run+3] = t4-t3;                      // Time for setting up the interval hash tables
  times[COMPUTE_BLOCKS*run+4] = t5-t4;                      // Time for computing the SPH particle densities

  return 0;
}

/*
 *  Function compute_density_3d_outerOmp:
 *    Computes the SPH density from the particles using cell linked list with
 *    vectorization at the compute_density_3d_chunk level, but the parallelization
 *    done at the level of the outer-most loop of the compute_density_3d_cll_outerOmp
 *    function, not at the chunk level.
 *
 *    Arguments:
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 */
int compute_density_3d_outerOmp(int N, double h, SPHparticle *lsph, linkedListBox *box){
  memset(lsph->rho,(int)0,N*sizeof(double));                                             // Pre-initialize the density to zero

  #pragma omp parallel for                                                                 // Execute the iteration in parallel
  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){ // Iterate over each receiver cell begin index
    int64_t node_hash=-1,node_begin=0, node_end=0;                                         // Start initializing the node indexes on the array
    int64_t nb_begin= 0, nb_end = 0;                                                       // initialize the neighbor indexes
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];                    // prepare a list of potential neighbor hashes

    if (kh_exist(box->hbegin, kbegin)){                                                    // verify if that given iterator actually exists
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));                  // Then get the end of the receiver cell iterator

      node_hash = kh_key(box->hbegin, kbegin);                                             // Then get the hash corresponding to it
      node_begin = kh_value(box->hbegin, kbegin);                                          // Get the receiver cell begin index in the array
      node_end   = kh_value(box->hend, kend);                                              // Get the receiver cell end index in the array

      neighbour_hash_3d(node_hash,nblist,box->width,box);                                  // then find the hashes of its neighbors
      for(int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){              // and the iterate over them
        if(nblist[j]>=0){                                                                  // if a given neighbor actually has particles
          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );            // then get the contributing cell begin index
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );            // and get the contributing cell end index

          compute_density_3d_chunk_noomp(node_begin,node_end,nb_begin,nb_end,h,            // and compute the density contribution from
                                         lsph->x,lsph->y,lsph->z,lsph->nu,lsph->rho);      // the contributing cell to the receiver cell
        }
      }
    }
  }

  return 0;
}

/*
 *  Function compute_density_3d_chunk_noomp:
 *    Computes the SPH density contribution for a pair of cells, from nb_ indexes
 *    to the node_ indexes. No parallelization is performed  with vectorization
 *    performed in the inner loop.
 *
 *    Arguments:
 *       node_begin <int>     : Begin index of the receiver cell
 *       node_end   <int>     : End   index of the receiver cell
 *       nb_begin <int>       : Begin index of the sender (neighbor) cell
 *       nb_end   <int>       : End   index of the sender (neighbor) cell
 *       h       <double>     : Smoothing Length for the Smoothing Kernel w_bspline
 *       x       <double*>    : Array of particles' X positions
 *       y       <double*>    : Array of particles' Y positions
 *       z       <double*>    : Array of particles' Z positions
 *       nu      <double*>    : Array of particles' density weights (i.e. masses)
 *    Returns:
 *       0                    : error code returned
 *       rho       <double*>  : Array of particles' densities
 */
int compute_density_3d_chunk_noomp(int64_t node_begin, int64_t node_end,
                                   int64_t nb_begin, int64_t nb_end,double h,
                                   double* restrict x, double* restrict y,
                                   double* restrict z, double* restrict nu,
                                   double* restrict rho){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);

  for(int64_t ii=node_begin;ii<node_end;ii+=1){           // Iterate over the ii index of the chunk
    double xii = x[ii];                                   // Load the X component of the ii particle position
    double yii = y[ii];                                   // Load the Y component of the ii particle position
    double zii = z[ii];                                   // Load the Z component of the ii particle position
    double rhoii = 0.0;                                   // Initialize the chunk contribution to density

    #pragma omp simd                                      // Hint at the compiler to vectorize
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){             // Iterate over the each other particle in jj loop
      double q = 0.;                                      // Initialize the distance

      double xij = xii-x[jj];                             // Load and subtract jj particle's X position component
      double yij = yii-y[jj];                             // Load and subtract jj particle's Y position component
      double zij = zii-z[jj];                             // Load and subtract jj particle's Z position component

      q += xij*xij;                                       // Add the jj contribution to the ii distance in X
      q += yij*yij;                                       // Add the jj contribution to the ii distance in X
      q += zij*zij;                                       // Add the jj contribution to the ii distance in X

      q = sqrt(q)*inv_h;                                  // Sqrt to compute the distance

      rhoii += nu[jj]*w_bspline_3d_simd(q);               // Add up the contribution from the jj particle
    }                                                     // to the intermediary density and then
    rho[ii] += rhoii*kernel_constant;                     // add the intermediary density to the full density
  }

  return 0;
}

/*
 *  Function w_bspline_3d_constant:
 *    Returns the 3d normalization constant for the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *    Returns:
 *       3d bspline normalization density <double>
 */
double w_bspline_3d_constant(double h){
  return 3./(2.*M_PI*h*h*h);  // 3d normalization value for the b-spline kernel
}

/*
 *  Function w_bspline_3d_simd:
 *    Returns the un-normalized value of the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       q <double>           : Distance between particles normalized by the smoothing length h
 *    Returns:
 *       wq <double>          : Unnormalized value of the kernel
 *
 *    Observation:
 *       Why not else if(q<2.)?
 *       Because if you use "else if", the compiler refuses to vectorize,
 *       This results in a large slowdown, as of 2.5x slower for example_04
 */
#pragma omp declare simd
double w_bspline_3d_simd(double q){
  double wq=0;
  double wq1 = (0.6666666666666666 - q*q + 0.5*q*q*q);             // The first polynomial of the spline
  double wq2 = 0.16666666666666666*(2.-q)*(2.-q)*(2.-q);           // The second polynomial of the spline

  if(q<2.)                                                         // If the distance is below 2
    wq = wq2;                                                      // Use the 2nd polynomial for the spline

  if(q<1.)                                                         // If the distance is below 1
    wq = wq1;                                                      // Use the 1st polynomial for the spline

  return wq;                                                       // return which ever value corresponds to the distance
}

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Begin constants */
#define VERBOSE 0 //when set to 1, running smc_sampler(...) will generate A LOT of output.
#define M 200
#define MH_STEPS 10
#define N_MIXTURE_COMPONENTS 4
#define SIGMA 0.55
#define HLINE printf("\n--------------------------------------------------------------------------------\n")

const double MIXTURE_WEIGHTS[N_MIXTURE_COMPONENTS] = {0.25,0.25,0.25,0.25};
static double* xNew; //an array used for swapping

/* begin function-signatures */
int test(void);
void allocateParticle_swap(void);
double dnorm(const double x,const double mu);
double log_likelihood(const double * y, const int N_y, const double* mu);
_Bool ESS_is_largeEnough(const double * W,const int N_W);
double pi_updateStep(const double* y,const int N_y, const double* x, const int n);
void smc_sampler(const double* yObs,const int N_yObs ,const int n_particles , double** particles, double* weights);
void randParticle(double * x);
void resample(double ** particles, double * weights,const int n_particles, const double w0, gsl_rng * r, gsl_ran_discrete_t * wSampler);
void propagate_particle(double * x,const int n, const double* yObs,const int N_yObs,gsl_rng * r);
void smc_sampler_for_R( double* yObs, int* N_yObs, int* N_particles, double* X_vec, double* W, double* time);
void freeParticle_swap(void);

int main(){
  allocateParticle_swap();
  return test();
  freeParticle_swap();
}

/* a wrapper function for calling smc_sampler from R */
void smc_sampler_for_R(
  double* yObs,
  int* N_yObs,
  int* N_particles,
  double* X_vec,
  double* W,
  double* t_SMC )
  {

  // allocate memmory for swapping
  allocateParticle_swap();

  // Generate an (N_particles x N_MIXTURE_COMPONENTS) array used to
  // internally represent particles
  double** X;
  X = (double**) malloc(sizeof(double *) * (*N_particles));
  size_t i = 0;
  for(i=0 ; i < *N_particles ; ++i){
    X[i] = (double*) malloc(sizeof(double)*N_MIXTURE_COMPONENTS);
  }

  //start timer
  time_t t0 = time(NULL);

  //run sampler
  smc_sampler(yObs,*N_yObs,*N_particles,X,W);

  //stop timer
  time_t t1 = time(NULL);

  //return runtime in secconds
  *t_SMC = difftime(t1,t0);

  //write to output vector
  size_t j;
  for( i = 0 ; i < (*N_particles) ; ++i){
    for( j = 0 ; j < N_MIXTURE_COMPONENTS ;  ++j){
      X_vec[N_MIXTURE_COMPONENTS * i + j] = X[i][j];
    }
  }

  //free up memory;
  for(i=0 ; i < *N_particles ; ++i){
    free(X[i]);
  }
  free(X);

  freeParticle_swap();

  /* //for debugging.
  int i;
  for(i=0 ; i < *N_particles ; ++i){
    X_vec[i] = (double) i;
  }
  */
}

/* Main function for running SMC sampler
 yObs is an a vector of N_yObs observations

 n_particles is the desired number of particles to be simulated

 particles should be a double[n_particles][N_MIXTURE_COMPONENTS] array.
 It will be used for storing results and doing computaations.

 weights should be a double[n_particles] array.
 It will be used for storing results and doing computations.
*/
void smc_sampler(
  const double* yObs,
  const int N_yObs,
  const int n_particles,
  double** particles,
  double* weights)
  {

  time_t t0 = time(NULL);
  time_t t1;
  size_t i; //summation indices

  gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937);
  gsl_ran_discrete_t * wSampler;

  if(VERBOSE) printf("Initializinng...");
  /* Initialize  weights and particles*/
  double w0 = 1 / (double) n_particles; // default weight
  for(i=0; i < n_particles ; ++i){
    weights[i] = w0;
    randParticle(particles[i]);
  }
  if(VERBOSE) printf("DONE!\n");

  if(VERBOSE) time(&t1);
  if(VERBOSE) printf("Elapsed time = %.0f sec\n",difftime(t1,t0));

  if(VERBOSE) printf("Running sampling steps from 1 to 200...\n");

  /* Iterate from n=1 to n=M */
  int n = 1;
  unsigned int resampleCounter = 0;
  while(1){

    // for diagnosing running time and output
    if((n-1)%10 == 0){
      if(VERBOSE) time(&t1);
      if(VERBOSE) printf("Running steps %i...%i; \t",n,n+9);
      if(VERBOSE) printf("elapsed time so far = %.0f sec\n",difftime(t1,t0));

      if(VERBOSE){
        printf("Printing the first 10 particles to check validity of output:\n");
        size_t j;
        for(i=0;i<10;++i){
          printf("  x[%i]_(%i) = ",(int) i,n);
          for(j=0; j<N_MIXTURE_COMPONENTS ; ++j){
            if(particles[i][j] >= 0) printf(" ");
            if(fabs(particles[i][j]) < 10) printf(" ");
            printf("%.3f\t",particles[i][j]);
          }
          printf(" \tW[%i](%i) = %f",(int) i , n, weights[i]);
          printf("\n");
        }
      }

      double sumW = 0;
      for(i = 0 ; i < n_particles ; ++i){
        sumW += weights[i];
      }
      if(VERBOSE) printf("sum of weights = %f",sumW);
      if(VERBOSE) HLINE;

     }

    /* Resample and reset weights, if nessecary */
    if(!(ESS_is_largeEnough(weights,n_particles))){
      if(VERBOSE) ++resampleCounter;
      if(VERBOSE) HLINE;
      if(VERBOSE) printf("resampling in step n = %i\n",n);
      if(VERBOSE) HLINE;
      //Preprocessing for sampling from a discrete distribution
      wSampler = gsl_ran_discrete_preproc (n_particles,weights);

      //Resample (AND reset) weights
      resample(particles,weights,n_particles,w0,r,wSampler);
    }

    /* Update particle weights */
    double sumW = 0;
    for(i = 0 ; i<n_particles ; ++i){
      weights[i] *= pi_updateStep(yObs, N_yObs, particles[i] ,n);
      sumW += weights[i];
    }
    double Z = (double) 1 / sumW;
    for(i = 0 ; i<n_particles ; ++i){
      weights[i] *= Z;
    }

    /* break out when n == M */
    if(n==M) break;

    /* Propagate particles */
    for(i=0; i<n_particles ; ++i){
      propagate_particle(particles[i], n, yObs, N_yObs,r);
    }

    /*increment counter*/
    ++n;
  }
  if(VERBOSE) printf("total number of resampling-steps = %i\n",resampleCounter);
  gsl_rng_free(r);
  // gsl_ran_discrete_free(wSampler);
}

/**********************
* Auxiliary Functions *
**********************/

void allocateParticle_swap(void){
  //modify in case of parralelization
  xNew = (double*) malloc(sizeof(double)*N_MIXTURE_COMPONENTS);
}

void freeParticle_swap(void){
  //modify in case of parralelization
  free(xNew);
}

double dnorm(const double x, const double mu){
  return  gsl_sf_erf_Z((x - mu)/SIGMA)/SIGMA;
}

/* returns log( p( y | mu , sigma , weights ) )
 y should be a vector of doubles corresponding to observations
 N_y corresponds to the length of y
 mu coresponds to a vector of modes of length N_MIXTURE_COMPONENTS */
double log_likelihood(const double* y, const int N_y,const double* mu){
  int i; int j; // initialize summation-indices

  /* Verify that all modes are in [-10,10] */
  for(i=0; i<N_MIXTURE_COMPONENTS; ++i){
    //INFINITY is a macro defined in math.h
    if( mu[i] > 10.0 || mu[i] < -10.0 ){
      if(VERBOSE) printf("WOAH! particle out of range\n");
      return -INFINITY;
    }
  }
  /* Compute the log-likelihood iteratively */
  double sum; // inner summation (over mu)
  double Sum = 0; //outer summation (over y)
  for(i = 0;i < N_y; ++i){
    sum = 0;
    for(j = 0 ; j < N_MIXTURE_COMPONENTS ; ++j){
      sum += MIXTURE_WEIGHTS[j] * dnorm(y[i], mu[j]);
    }
    Sum += log(sum);
  }
  return Sum;
}

/* Generates a ranodm particle ie. 4 ind. samples uniformly from -10 to 10 */
/* Results are stored in x */
void randParticle(double* x){
  unsigned char i;
  for(i=0 ; i<N_MIXTURE_COMPONENTS; ++i){
    x[i] = (double) -10 + (double) 20 * (double) rand() / (double) RAND_MAX;
  }
}

_Bool ESS_is_largeEnough(const double * W,const int N_W){
  double sumSq = 0;
  int i;
  for(i = 0 ; i < N_W ; ++i){
    sumSq += W[i] * W[i];
  }
  return (_Bool) 1 >= sumSq * N_W / 2.0;
}

void resample(
  double ** particles,
  double * weights,
  const int n_particles,
  const double w0,
  gsl_rng * r,
  gsl_ran_discrete_t * wSampler)
  {
  size_t i; size_t j;
  for(i=0; i< n_particles; ++i){
    //iterate over particles

    size_t i_new = gsl_ran_discrete(r,wSampler);
    // copy the particle in question
    for(j = 0 ; j < N_MIXTURE_COMPONENTS ; ++j){
      particles[i][j] = particles[i_new][j];
    }
  }
  //reset particle weights
  for(i=0; i < n_particles ; ++i){
    weights[i] = w0;
  }
}

/*
A n auxiliary function for computing pi_n(x_(n-1)) / pi_(n-1)(x_(n-1))
\eqn{\propto} p( x_(n-1) | y ) ^ ((2n -1)/M^2). Used when updating weights.
@param y Observed data
@param x Our proposed value for Mu (a vector of modes)
*/
double pi_updateStep(const double* y,const int N_y, const double* x,const int n){
  return exp( ( ((float) 2*n -1)/((float) (M*M)) ) * log_likelihood(y,N_y,x));
}

void propagate_particle(double * x,const int n, const double* yObs,const int N_yObs,gsl_rng * r){

  //NOTE: xNew MUST have been allocated before calling this method!

  size_t i = 0; size_t j;
  for(i = 0 ; i < MH_STEPS ; ++i){

    // Generate new observation
    for(j = 0 ; j < N_MIXTURE_COMPONENTS ; ++j){
      xNew[j] = x[j] + gsl_ran_gaussian(r,SIGMA);
    }

/*
    // REMOVE WHEN COMPILING FOR SPEED!
    if(VERBOSE){
      // printf("Acceptance! alpha = %0.3f\n",alpha);
      printf("  x\t= ");
      for(j=0; j<N_MIXTURE_COMPONENTS ; ++j){
        if(x[j] >= 0) printf(" ");
        if(fabs(x[j]) < 10) printf(" ");
        printf("%.3f\t",x[j]);
      }
      printf("\n  x_new\t=");
      for(j=0; j<N_MIXTURE_COMPONENTS ; ++j){
        if(xNew[j] >= 0) printf(" ");
        if(fabs(xNew[j]) < 10) printf(" ");
        printf("%.3f\t",xNew[j]);
      }
      printf("\n");
    }
*/

    // double logL_new = log_likelihood(yObs, N_yObs, xNew);
    // double logL_old = log_likelihood(yObs, N_yObs, x);
    double logL_Diff = log_likelihood(yObs, N_yObs, xNew) - log_likelihood(yObs, N_yObs, x);
    double alpha = 1;
    _Bool autoAccept = 1;

    // Iff logL_Diff < 0 true, the sign of the argument fo exp() will be <1.
    if(logL_Diff < 0){
      double exponent = (double) n / (double) M;
      exponent *= exponent; // square the exponent
      alpha = exp(exponent * logL_Diff );
      autoAccept = 0;
      if(VERBOSE) printf("alpha = %.4f",alpha);
    }

    if( autoAccept || ((double) rand() / (double) RAND_MAX) < alpha){
      if( VERBOSE && !autoAccept ) printf("\t new point accepted anyway\n");
      for(j = 0 ; j < N_MIXTURE_COMPONENTS ; ++j){
        x[j] = xNew[j];
      }
    } else {
      if(VERBOSE) printf("\n");
    }
  }
  if(VERBOSE) printf("\n");
}

/* A generic function to test during development
Returns 0 if test succcessfull*/
int test(void){
/*
  double x = 1; double y;
  double mu = 0; double sigma = 5;
  y = dnorm(x,mu,sigma);
  printf("Hello World!\nNormal(1;0,5) = %.10f\n",y);
*/

/*
  double y[2] = {0,2};
  double mu[4] = {-3,0,3,6};
  double x = log_likelihood(y,2,mu);
  printf("Hello World!\nx = %.10f\n",x);
*/

/*
  double mu[4] = {-3,0,3,6};
  if( mu[2] < 10 ) printf("Hello World!\n");
*/

/*
  int N_W = 6;
  double W1[6] = {1,0,0,0,0,0};
  double W2[6] = {0.2,0.2,0.2,0.2,0.2,0};

  if(!(ESS_is_largeEnough(W1,N_W))) printf("ESS1 is too small!\n");
  if(ESS_is_largeEnough(W2,N_W)) printf("ESS2 on the other hand is groovy!\n");
*/

/*
  double y[2] = {0,2};
  double mu[4] = {-3,0,3,6};
  double x = pi_updateStep(y,2,mu,180);
  printf("Hello World!\nx = %.10f\n",x);
*/

/*
  int i; int j;
  double particles[10][4];
  for(i = 0 ; i<10; ++i){
    randParticle(particles[i]);
  }

  for(i=0 ; i<10 ; ++i){
    printf("rand particles = (");
    for(j=0 ; j< 4 ; ++j){
      printf(" %.2f ",particles[i][j]);
    }
    printf( ")\n");
  }
*/
/*
  int N = 1000;

  double* W;
  W = (double*) malloc(sizeof(double)*N);

  double** X;
  X = (double**) malloc(sizeof(double *) * N);
  size_t i = 0;
  for(i=0 ; i < N ; ++i){
    X[i] = (double*) malloc(sizeof(double*)*N_MIXTURE_COMPONENTS);
  }

  // double W;
  // double* X[N][N_MIXTURE_COMPONENTS];

  // data pasted in from R
  // generated using sampleMM(100) from mixtureModel_SMC.R
  int N_yObs = 100;
  double yObs[100] =
   {0.270469630282116 , 0.0948544266135598 , 4.82353957890236 ,
    2.71176569013062 , -2.10140287048951 , -2.4639397242585 ,
    -2.47096320311181 , 5.59163596897931 , -3.50775981783882 ,
    3.26972797997964 , 3.32256384988976 , 5.73616849148524 ,
    5.74816296903159 , -2.97620992270215 , -0.211030444261816 ,
    5.56902016608288 , -3.69397287391297 , 5.71718272922272 ,
    -0.375056085230512 , 0.290727561775271 , 3.53116678441137 ,
    -2.92396872511173 , 5.44623183851942 , 5.60306594746839 ,
    5.56855197639855 , 2.28105552916777 , 2.71638272346431 ,
    3.7321825583144 , 6.10978497180474 , 6.7089442115652 ,
    3.14182326021059 , -2.72341774114777 , -0.0919451805310823 ,
    -3.05334875104901 , 6.00697089830724 , -2.81827687988051 ,
    6.03028599825704 , 3.77081092101111 , -0.17554472859167 ,
    5.57845673340545 , 6.25654424410659 , 6.40046681340197 ,
    -2.39322730700275 , 0.937990298556008 , -3.61930611127869 ,
    2.99419558349795 , -2.95427533963002 , 2.53999917081115 ,
    -4.35170347149833 , 2.56389630336686 , 2.78219272857881 ,
    -2.53796084626858 , 2.75783249644654 , 0.265653681397089 ,
    -4.08489661583957 , 0.141285949971235 , -3.09744117607388 ,
    -3.29683116823359 , -0.56348650880506 , 2.13298104427159 ,
    5.81313774243827 , 2.90131666078861 , -0.870942485382193 ,
    -0.653182465687133 , 3.24697479240955 , 3.46323186816872 ,
    3.70853665171654 , -3.79582094987125 , 6.55169236217858 ,
    -3.63836564940765 , 0.132152004366257 , 5.24827155186988 ,
    2.53806422373366 , 5.6550753862474 , 2.42737960470922 ,
    5.88968129179486 , 6.83133504027115 , -0.473386824297909 ,
    -0.883286661476473 , 6.94044430889499 , 2.73913195305667 ,
    3.02781182874733 , 6.53392928375413 , 6.53623960926176 ,
    6.44586265929949 , 0.737757409838858 , -2.81327223680057 ,
    0.534431297323472 , -3.28504638650391 , -0.0839400339092673 ,
    -3.30271081947717 , 5.97247453668681 , -3.55468774001656 ,
    -0.211152666301356 , 6.12181217531004 , 2.56643856048767 ,
    0.200557127297898 , -3.97176332345347 , 0.0176492920069539 ,
    2.87381823065611};

  printf("Running SMC-sampler with %i particles, based on %i observations.\n",N,N_yObs);
  time_t t1 = time(NULL);
  smc_sampler(yObs,N_yObs,N,X,W);
  time_t t2 = time(NULL);
  printf("DONE!\nTotal elapsed time = %.0f\n",difftime(t2,t1));

  HLINE;

  printf("Printing the first 10 particles to check validity of output:\n");
  size_t j;
  for(i=0;i<10;++i){
    printf("  x[%i] = ",(int) i);
    for(j=0; j<N_MIXTURE_COMPONENTS ; ++j){
      if(X[i][j] >= 0) printf(" ");
      if(fabs(X[i][j]) < 10) printf(" ");
      printf("%.3f\t",X[i][j]);
    }
    printf("\n");
  }

  free(W);
  // free(X)
  for(i=0 ; i<N; ++i){
    free(X[i]);
  }
  free(X);
*/

/*
  int N_yObs = 100;

  //DATA from R, so that the S-method and the .c method can be directly compated
  double yObs[100] = {6.42614409154078,0.260646235868524,6.20077626783781,-3.54124388162526,3.16810352365003,-3.02631343546788,6.22743478871312,6.8131742684717,6.24776924281844,3.65204896173014,6.01102502511306,-2.78204404408537,6.54496067301763,2.86661732139087,6.41368106980899,-2.54549938986914,6.08430821908108,2.65819528014716,-4.10369954335742,-3.70836228135375,3.55973960494948,-3.99408019214758,-3.06751044524848,-0.560058224387682,0.625592793169344,2.58052593365431,0.689049384542289,0.494930560119956,0.390397430673082,2.87213201595585,6.08140477180053,-3.63073724259346,1.76368358083185,0.0978387885170348,6.49576360025151,6.65592881975105,-2.31116539097971,6.69822824187681,0.0197128170022586,-3.6584094557839,6.37219590110012,5.98608504754753,5.76477302651333,-2.49473884182791,3.68861876383596,6.71753596416937,-2.86065603329243,2.36234293889529,-2.13229209410826,-3.1821797074668,3.67640450999392,2.35200323180728,-2.90422204445312,3.85237110998024,5.73406945312462,3.61994874049274,-4.0561902744938,5.66696946993138,5.96245114627653,6.59816990256299,0.231879318383607,2.34472949459441,2.72937582033292,-0.444334698629388,2.22482121394123,3.15359838488669,5.79310254628412,2.95627482456672,0.0754282164405023,-3.27128946837834,-2.77302236660845,7.25710145940199,6.2324724323204,5.9296066592542,5.60999425487975,-0.286465744070193,0.587073138580418,-2.67791031284446,6.98129381109518,3.36137498170339,-3.16490919686681,-3.62156309636938,5.41982668735859,-3.07991914749654,-2.89255426727392,5.61010574600101,-3.10612256078391,3.48991878214866,2.58896632171936,0.703862729104843,4.72001741709978,-0.244802119787893,6.45395925612269,0.253461730596016,2.69213945052578,-2.90338536462303,3.28208239426602,-2.57167632631676,0.920364746339414,0.993277872559861};

  int N_X = 5;
  double X[5][4] = {
    {0,3,-3,6},
    {-3,0,3,6},
    {0,0,0,0},
    {9,9,9,9},
    {0,6,3,15}
  };

  size_t i;
  for(i = 0 ; i<N_X ; ++i){
    double LL = log_likelihood(yObs,N_yObs,X[i]);
    printf("log_L(X[%i]) = %f\n",(int) i,LL);
  }
*/
  return 0;
}