text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

#ifndef setuparea_h
#define setuparea_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include <stdbool.h>
#include <string.h>
#include "qfunctions/area/areacube.h"
#include "qfunctions/area/areasphere.h"

#if PETSC_VERSION_LT(3,14,0)
#  define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i)
#  define DMPlexRestoreClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexRestoreClosureIndices(a,b,c,d,f,g,i)
#endif

#if PETSC_VERSION_LT(3,14,0)
#  define DMPlexCreateSphereMesh(a,b,c,d,e) DMPlexCreateSphereMesh(a,b,c,e)
#endif

// -----------------------------------------------------------------------------
// PETSc Operator Structs
// -----------------------------------------------------------------------------

// Data for PETSc
typedef struct User_ *User;
struct User_ {
  MPI_Comm comm;
  DM dm;
  Vec Xloc, Yloc, diag;
  CeedVector xceed, yceed;
  CeedOperator op;
  Ceed ceed;
};

// -----------------------------------------------------------------------------
// libCEED Data Struct
// -----------------------------------------------------------------------------

// libCEED data struct for level
typedef struct CeedData_ *CeedData;
struct CeedData_ {
  Ceed ceed;
  CeedBasis basisx, basisu;
  CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi;
  CeedQFunction qf_apply;
  CeedOperator op_apply, op_restrict, op_interp;
  CeedVector qdata, uceed, vceed;
};

// -----------------------------------------------------------------------------
// Problem Option Data
// -----------------------------------------------------------------------------

// Problem options
typedef enum {
  CUBE = 0, SPHERE = 1
} problemType;
static const char *const problemTypes[] = {"cube", "sphere",
                                           "problemType", "AREA", NULL
                                          };

// Problem specific data
typedef struct {
  CeedInt ncompu, ncompx, qdatasize, qextra, topodim;
  CeedQFunctionUser setupgeo, apply;
  const char *setupgeofname, *applyfname;
  CeedEvalMode inmode, outmode;
  CeedQuadMode qmode;
} problemData;

static problemData problemOptions[6] = {
  [CUBE] = {
    .ncompx = 3,
    .ncompu = 1,
    .topodim = 2,
    .qdatasize = 1,
    .qextra = 1,
    .setupgeo = SetupMassGeoCube,
    .apply = Mass,
    .setupgeofname = SetupMassGeoCube_loc,
    .applyfname = Mass_loc,
    .inmode = CEED_EVAL_INTERP,
    .outmode = CEED_EVAL_INTERP,
    .qmode = CEED_GAUSS
  },
  [SPHERE] = {
    .ncompx = 3,
    .ncompu = 1,
    .topodim = 2,
    .qdatasize = 1,
    .qextra = 1,
    .setupgeo = SetupMassGeoSphere,
    .apply = Mass,
    .setupgeofname = SetupMassGeoSphere_loc,
    .applyfname = Mass_loc,
    .inmode = CEED_EVAL_INTERP,
    .outmode = CEED_EVAL_INTERP,
    .qmode = CEED_GAUSS
  }
};

// -----------------------------------------------------------------------------
// PETSc sphere auxiliary functions
// -----------------------------------------------------------------------------

// Utility function taken from petsc/src/dm/impls/plex/examples/tutorials/ex7.c
static PetscErrorCode ProjectToUnitSphere(DM dm) {
  Vec            coordinates;
  PetscScalar   *coords;
  PetscInt       Nv, v, dim, d;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr);
  ierr = VecGetLocalSize(coordinates, &Nv); CHKERRQ(ierr);
  ierr = VecGetBlockSize(coordinates, &dim); CHKERRQ(ierr);
  Nv  /= dim;
  ierr = VecGetArray(coordinates, &coords); CHKERRQ(ierr);
  for (v = 0; v < Nv; ++v) {
    PetscReal r = 0.0;

    for (d = 0; d < dim; ++d) r += PetscSqr(PetscRealPart(coords[v*dim+d]));
    r = PetscSqrtReal(r);
    for (d = 0; d < dim; ++d) coords[v*dim+d] /= r;
  }
  ierr = VecRestoreArray(coordinates, &coords); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// PETSc Finite Element space setup
// -----------------------------------------------------------------------------

// Create FE by degree
static int PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc,
                                 PetscBool isSimplex, const char prefix[],
                                 PetscInt order, PetscFE *fem) {
  PetscQuadrature q, fq;
  DM              K;
  PetscSpace      P;
  PetscDualSpace  Q;
  PetscInt        quadPointsPerEdge;
  PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
  PetscErrorCode  ierr;

  PetscFunctionBeginUser;
  /* Create space */
  ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix); CHKERRQ(ierr);
  ierr = PetscSpacePolynomialSetTensor(P, tensor); CHKERRQ(ierr);
  ierr = PetscSpaceSetFromOptions(P); CHKERRQ(ierr);
  ierr = PetscSpaceSetNumComponents(P, Nc); CHKERRQ(ierr);
  ierr = PetscSpaceSetNumVariables(P, dim); CHKERRQ(ierr);
  ierr = PetscSpaceSetDegree(P, order, order); CHKERRQ(ierr);
  ierr = PetscSpaceSetUp(P); CHKERRQ(ierr);
  ierr = PetscSpacePolynomialGetTensor(P, &tensor); CHKERRQ(ierr);
  /* Create dual space */
  ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
  CHKERRQ(ierr);
  ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); CHKERRQ(ierr);
  ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetDM(Q, K); CHKERRQ(ierr);
  ierr = DMDestroy(&K); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetNumComponents(Q, Nc); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetOrder(Q, order); CHKERRQ(ierr);
  ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetFromOptions(Q); CHKERRQ(ierr);
  ierr = PetscDualSpaceSetUp(Q); CHKERRQ(ierr);
  /* Create element */
  ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem); CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); CHKERRQ(ierr);
  ierr = PetscFESetFromOptions(*fem); CHKERRQ(ierr);
  ierr = PetscFESetBasisSpace(*fem, P); CHKERRQ(ierr);
  ierr = PetscFESetDualSpace(*fem, Q); CHKERRQ(ierr);
  ierr = PetscFESetNumComponents(*fem, Nc); CHKERRQ(ierr);
  ierr = PetscFESetUp(*fem); CHKERRQ(ierr);
  ierr = PetscSpaceDestroy(&P); CHKERRQ(ierr);
  ierr = PetscDualSpaceDestroy(&Q); CHKERRQ(ierr);
  /* Create quadrature */
  quadPointsPerEdge = PetscMax(order + 1,1);
  if (isSimplex) {
    ierr = PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0,
                                          &q); CHKERRQ(ierr);
    ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0,
                                          &fq); CHKERRQ(ierr);
  } else {
    ierr = PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0,
                                        &q); CHKERRQ(ierr);
    ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0,
                                        &fq); CHKERRQ(ierr);
  }
  ierr = PetscFESetQuadrature(*fem, q); CHKERRQ(ierr);
  ierr = PetscFESetFaceQuadrature(*fem, fq); CHKERRQ(ierr);
  ierr = PetscQuadratureDestroy(&q); CHKERRQ(ierr);
  ierr = PetscQuadratureDestroy(&fq); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// PETSc Setup for Level
// -----------------------------------------------------------------------------

// This function sets up a DM for a given degree
static int SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompu,
                           PetscInt dim) {
  PetscInt ierr;
  PetscFE fe;

  PetscFunctionBeginUser;

  // Setup FE
  ierr = PetscFECreateByDegree(dm, dim, ncompu, PETSC_FALSE, NULL, degree, &fe);
  CHKERRQ(ierr);
  ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr);

  // Setup DM
  ierr = DMCreateDS(dm); CHKERRQ(ierr);
  ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);
  CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// -----------------------------------------------------------------------------
// libCEED Setup for Level
// -----------------------------------------------------------------------------

// Destroy libCEED operator objects
static PetscErrorCode CeedDataDestroy(CeedData data) {
  PetscInt ierr;

  CeedVectorDestroy(&data->qdata);
  CeedVectorDestroy(&data->uceed);
  CeedVectorDestroy(&data->vceed);
  CeedBasisDestroy(&data->basisx);
  CeedBasisDestroy(&data->basisu);
  CeedElemRestrictionDestroy(&data->Erestrictu);
  CeedElemRestrictionDestroy(&data->Erestrictx);
  CeedElemRestrictionDestroy(&data->Erestrictqdi);
  CeedQFunctionDestroy(&data->qf_apply);
  CeedOperatorDestroy(&data->op_apply);
  ierr = PetscFree(data); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Auxiliary function to define CEED restrictions from DMPlex data
static int CreateRestrictionPlex(Ceed ceed, CeedInt P, CeedInt ncomp,
                                 CeedElemRestriction *Erestrict, DM dm) {
  PetscInt ierr;
  PetscInt c, cStart, cEnd, nelem, nnodes, *erestrict, eoffset;
  PetscSection section;
  Vec Uloc;

  PetscFunctionBegin;

  // Get Nelem
  ierr = DMGetSection(dm, &section); CHKERRQ(ierr);
  ierr = DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd); CHKERRQ(ierr);
  nelem = cEnd - cStart;

  // Get indices
  ierr = PetscMalloc1(nelem*P*P, &erestrict); CHKERRQ(ierr);
  for (c=cStart, eoffset = 0; c<cEnd; c++) {
    PetscInt numindices, *indices, i;
    ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE,
                                   &numindices, &indices, NULL, NULL);
    CHKERRQ(ierr);
    for (i=0; i<numindices; i+=ncomp) {
      for (PetscInt j=0; j<ncomp; j++) {
        if (indices[i+j] != indices[i] + (PetscInt)(copysign(j, indices[i])))
          SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
                   "Cell %D closure indices not interlaced", c);
      }
      // NO BC on closed surfaces
      PetscInt loc = indices[i];
      erestrict[eoffset++] = loc;
    }
    ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE,
                                       &numindices, &indices, NULL, NULL);
    CHKERRQ(ierr);
  }

  // Setup CEED restriction
  ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr);
  ierr = VecGetLocalSize(Uloc, &nnodes); CHKERRQ(ierr);

  ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr);
  CeedElemRestrictionCreate(ceed, nelem, P*P, ncomp, 1, nnodes,
                            CEED_MEM_HOST, CEED_COPY_VALUES, erestrict,
                            Erestrict);
  ierr = PetscFree(erestrict); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Set up libCEED for a given degree
static int SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
                                CeedInt topodim, CeedInt qextra,
                                PetscInt ncompx, PetscInt ncompu,
                                PetscInt xlsize, problemType problemChoice,
                                CeedData data) {
  int ierr;
  DM dmcoord;
  Vec coords;
  const PetscScalar *coordArray;
  CeedBasis basisx, basisu;
  CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi;
  CeedQFunction qf_setupgeo, qf_apply;
  CeedOperator op_setupgeo, op_apply;
  CeedVector xcoord, qdata, uceed, vceed;
  CeedInt P, Q, cStart, cEnd, nelem,
          qdatasize = problemOptions[problemChoice].qdatasize;

  // CEED bases
  P = degree + 1;
  Q = P + qextra;
  CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompu, P, Q,
                                  problemOptions[problemChoice].qmode, &basisu);
  CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompx, 2, Q,
                                  problemOptions[problemChoice].qmode, &basisx);

  // CEED restrictions
  ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr);
  ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL);
  CHKERRQ(ierr);

  ierr = CreateRestrictionPlex(ceed, 2, ncompx, &Erestrictx, dmcoord);
  CHKERRQ(ierr);
  ierr = CreateRestrictionPlex(ceed, P, ncompu, &Erestrictu, dm); CHKERRQ(ierr);

  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
  nelem = cEnd - cStart;

  CeedElemRestrictionCreateStrided(ceed, nelem, Q*Q, qdatasize,
                                   qdatasize*nelem*Q*Q,
                                   CEED_STRIDES_BACKEND, &Erestrictqdi);

  // Element coordinates
  ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr);
  ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr);

  CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);
  CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES,
                     (PetscScalar *)coordArray);
  ierr = VecRestoreArrayRead(coords, &coordArray);

  // Create the vectors that will be needed in setup and apply
  CeedInt nqpts;
  CeedBasisGetNumQuadraturePoints(basisu, &nqpts);
  CeedVectorCreate(ceed, qdatasize*nelem*nqpts, &qdata);
  CeedVectorCreate(ceed, xlsize, &uceed);
  CeedVectorCreate(ceed, xlsize, &vceed);

  // Create the Q-function that builds the operator (i.e. computes its
  // quadrature data) and set its context data
  CeedQFunctionCreateInterior(ceed, 1, problemOptions[problemChoice].setupgeo,
                              problemOptions[problemChoice].setupgeofname,
                              &qf_setupgeo);
  CeedQFunctionAddInput(qf_setupgeo, "x", ncompx, CEED_EVAL_INTERP);
  CeedQFunctionAddInput(qf_setupgeo, "dx", ncompx*topodim, CEED_EVAL_GRAD);
  CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT);
  CeedQFunctionAddOutput(qf_setupgeo, "qdata", qdatasize, CEED_EVAL_NONE);

  // Set up the mass operator
  CeedQFunctionCreateInterior(ceed, 1, problemOptions[problemChoice].apply,
                              problemOptions[problemChoice].applyfname,
                              &qf_apply);
  CeedQFunctionAddInput(qf_apply, "u", ncompu,
                        problemOptions[problemChoice].inmode);
  CeedQFunctionAddInput(qf_apply, "qdata", qdatasize, CEED_EVAL_NONE);
  CeedQFunctionAddOutput(qf_apply, "v", ncompu,
                         problemOptions[problemChoice].outmode);

  // Create the operator that builds the quadrature data for the operator
  CeedOperatorCreate(ceed, qf_setupgeo, NULL, NULL, &op_setupgeo);
  CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx,
                       CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx,
                       CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx,
                       CEED_VECTOR_NONE);
  CeedOperatorSetField(op_setupgeo, "qdata", Erestrictqdi,
                       CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

  // Create the mass or diff operator
  CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply);
  CeedOperatorSetField(op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_apply, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED,
                       qdata);
  CeedOperatorSetField(op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE);

  // Setup qdata
  CeedOperatorApply(op_setupgeo, xcoord, qdata, CEED_REQUEST_IMMEDIATE);

  // Cleanup
  CeedQFunctionDestroy(&qf_setupgeo);
  CeedOperatorDestroy(&op_setupgeo);
  CeedVectorDestroy(&xcoord);

  // Save libCEED data
  data->basisx = basisx;
  data->basisu = basisu;
  data->Erestrictx = Erestrictx;
  data->Erestrictu = Erestrictu;
  data->Erestrictqdi = Erestrictqdi;
  data->qf_apply = qf_apply;
  data->op_apply = op_apply;
  data->qdata = qdata;
  data->uceed = uceed;
  data->vceed = vceed;

  PetscFunctionReturn(0);
}

#endif // setuparea_h


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
#ifdef COSMIC_RAYS
#include "cosmic_rays.h"
#endif

void reconstruct_timebins(void)
{
  int i, n, prev, bin;
  long long glob_sum1, glob_sum2;

  for(bin = 0; bin < TIMEBINS; bin++)
    {
      TimeBinCount[bin] = 0;
      TimeBinCountSph[bin] = 0;
      FirstInTimeBin[bin] = -1;
      LastInTimeBin[bin] = -1;
#ifdef SFR
      TimeBinSfr[bin] = 0;
#endif
#ifdef BLACK_HOLES
      TimeBin_BH_mass[bin] = 0;
      TimeBin_BH_dynamicalmass[bin] = 0;
      TimeBin_BH_Mdot[bin] = 0;
      TimeBin_BH_Medd[bin] = 0;
#endif
    }

  for(i = 0; i < NumPart; i++)
    {
      bin = P[i].TimeBin;

      if(TimeBinCount[bin] > 0)
	{
	  PrevInTimeBin[i] = LastInTimeBin[bin];
	  NextInTimeBin[i] = -1;
	  NextInTimeBin[LastInTimeBin[bin]] = i;
	  LastInTimeBin[bin] = i;
	}
      else
	{
	  FirstInTimeBin[bin] = LastInTimeBin[bin] = i;
	  PrevInTimeBin[i] = NextInTimeBin[i] = -1;
	}
      TimeBinCount[bin]++;
      if(P[i].Type == 0)
	TimeBinCountSph[bin]++;

#ifdef SFR
      if(P[i].Type == 0)
	TimeBinSfr[bin] += SphP[i].Sfr;
#endif
#if BLACK_HOLES
      if(P[i].Type == 5)
	{
	  TimeBin_BH_mass[bin] += P[i].BH_Mass;
	  TimeBin_BH_dynamicalmass[bin] += P[i].Mass;
	  TimeBin_BH_Mdot[bin] += P[i].BH_Mdot;
	  TimeBin_BH_Medd[bin] += P[i].BH_Mdot / P[i].BH_Mass;
	}
#endif
    }

  FirstActiveParticle = -1;

  for(n = 0, prev = -1; n < TIMEBINS; n++)
    {
      if(TimeBinActive[n])
	for(i = FirstInTimeBin[n]; i >= 0; i = NextInTimeBin[i])
	  {
	    if(prev == -1)
	      FirstActiveParticle = i;

	    if(prev >= 0)
	      NextActiveParticle[prev] = i;

	    prev = i;
	  }
    }

  if(prev >= 0)
    NextActiveParticle[prev] = -1;

  sumup_large_ints(1, &NumForceUpdate, &glob_sum1);

  for(i = FirstActiveParticle, NumForceUpdate = 0; i >= 0; i = NextActiveParticle[i])
    {
      NumForceUpdate++;
      if(i >= NumPart)
	{
	  printf("Bummer i=%d\n", i);
	  endrun(12);
	}

    }

  sumup_large_ints(1, &NumForceUpdate, &glob_sum2);

  if(ThisTask == 0)
    {
      printf("sum1=%d%9d sum2=%d%9d\n",
	     (int) (glob_sum1 / 1000000000), (int) (glob_sum1 % 1000000000),
	     (int) (glob_sum2 / 1000000000), (int) (glob_sum2 % 1000000000));
    }

  if(glob_sum1 != glob_sum2 && All.NumCurrentTiStep > 0)
    endrun(121);
}

void drift_particle(int i, int time1)
{
  int j, time0, dt_step;
  double dt_drift, dt_gravkick, dt_hydrokick, dt_entr;

#ifdef DISTORTIONTENSORPS
  int j1, j2;
  /*
   Setup all numbers that are going to be written to caustic log file. Some of them will not be set, depending
   on Makefile options. Therefore we assign zero to each value in the beginning.
  */
  MyDouble s_1=0.0, s_2=0.0, s_3=0.0;
  MyDouble smear_x=0.0, smear_y=0.0, smear_z=0.0;
  MyDouble CurrentTurnaroundRadius=0.0;
  MyDouble caustic_counter = 0.0;
  MyDouble summed_annihilation = 0.0;
  MyDouble D_xP = 0.0;
  MyDouble init_density = 0.0;
  MyDouble rho_normed_cutoff = 0.0;
  MyDouble second_deriv_fac = 0.0;
  MyDouble analytic_add_to_annihilation = 0.0;

  /* things needed for caustics detection */
  MyDouble determinant =  1.0;
  MyDouble product_matrix[3][3];
  MyDouble current_normed_stream_density;

  MyDouble last_normed_stream_density;

#ifdef ANNIHILATION_RADIATION
  /*1: if annihilation radiation was integrated analytically through the caustic, 0: otherwise */
  int caustic_annihilation_flag = 0;
  MyDouble rho_0, kappa_caustic;
#endif

#if SIM_ADAPTIVE_SOFT
  CurrentTurnaroundRadius = All.CurrentTurnaroundRadius;
#endif
#endif /* DISTORTIONTENSORPS */

#ifdef COSMIC_RAYS
  int CRpop;
#endif

#ifdef SOFTEREQS
  double a3inv;

  if(All.ComovingIntegrationOn)
    a3inv = 1 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1;
#endif

  time0 = P[i].Ti_current;

  if(time1 < time0)
    {
      printf("i=%d time0=%d time1=%d\n", i, time0, time1);
      endrun(12);
    }

  if(time1 == time0)
    return;

  if(All.ComovingIntegrationOn)
    {
      dt_drift = get_drift_factor(time0, time1);
      dt_gravkick = get_gravkick_factor(time0, time1);
      dt_hydrokick = get_hydrokick_factor(time0, time1);
    }
  else
    {
      dt_drift = dt_gravkick = dt_hydrokick = (time1 - time0) * All.Timebase_interval;
    }

  for(j = 0; j < 3; j++)
    P[i].Pos[j] += P[i].Vel[j] * dt_drift;

/* START PHASE-SPACE ANALYSIS ------------------------------------------------------ */
#ifdef DISTORTIONTENSORPS

  /* do the DRIFT for distortion (this updates D_xq and D_xp; D_vq and D_vp are updated during kick) */
  for(j1 = 0; j1 < 3; j1++)
    for(j2 = 0; j2 < 6; j2++)
      P[i].distortion_tensorps[j1][j2] += P[i].distortion_tensorps[j1 + 3][j2] * dt_drift;

#ifdef REINIT_AT_TURNAROUND
  /* do we have to reinit distortion tensor? Only particles that are initially outside the turnaround radius are analyzed */
  if ((P[i].Pos[0]*P[i].Vel[0]+P[i].Pos[1]*P[i].Vel[1]+P[i].Pos[2]*P[i].Vel[2] < 0.0) && (P[i].turnaround_flag == 0))
    {
      /* yes, particle turned around */
      /* re-init distortion tensor */
      for(j1 = 0; j1 < 6; j1++)
	for(j2 = 0; j2 < 6; j2++)
	  {
	    if((j1 == j2))
	      {
		P[i].distortion_tensorps[j1][j2] = 1.0;
	      }
	    else
	      {
		P[i].distortion_tensorps[j1][j2] = 0.0;
	      }
	  }
      /* mark particle so that we know that it already turned around */
      P[i].turnaround_flag = 1;

      /* set new initial sheet orientation since we start distortion from this point (analytic solution) */
      for (j1=0; j1<=2; j1++)
        for (j2=0; j2<=2; j2++)
          P[i].V_matrix[j1][j2] = P[i].Pos[j1]*P[i].Pos[j2] / (P[i].Pos[0]*P[i].Pos[0]+P[i].Pos[1]*P[i].Pos[1]+P[i].Pos[2]*P[i].Pos[2]) *
                                  (3.0/4.0*M_PI)*(3.0/4.0*M_PI) / (3.0 + 1.0/All.SIM_epsilon) * 1.0 / All.Time;
      /*
       The correct stream init density at turn around is given by rho_crit * 9.0 * M_PI * M_PI / (16.0*(3.0*epsilon+1.0)).
       This is exactly the density at the turnaround radius for a given density perturbation index epsilon.
       (see for example, Sikivie et al 1997)
      */
      P[i].init_density = 9.0 * M_PI * M_PI / (16.0*(3.0*All.SIM_epsilon+1.0)) * 1.0 / (6.0 * M_PI * All.G * All.Time * All.Time);

      /* reset last_stream_determinant (needed for caustic identification) */
      P[i].last_stream_determinant = 1.0;
#ifdef ANNIHILATION_RADIATION
      /* set annihilation radiation to zero */
      P[i].annihilation = 0.0;
#endif
    }
#endif

  /* NOW CAUSTIC AND ANNIHILATION HANDLING */
#ifdef REINIT_AT_TURNAROUND
 /* only analyse particle initially outside the turnaround radius, those initially inside have been set to P[i].turnaround_flag=-1 */
 if (P[i].turnaround_flag >= 0)
 {
#endif
  /* CAUSTIC IDENTIFICATION PART */
  /* projection phase-space distrtion matrix (see Vogelsberger et al, 2008 Eq. (22)) -> needed to check for caustic in drift */
  product_matrix[0][0] = P[i].distortion_tensorps[0][0] +
                         P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][0] +
                         P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][0] +
                         P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][0];
  product_matrix[0][1] = P[i].distortion_tensorps[0][1] +
                         P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][1] +
                         P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][1] +
                         P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][1];
  product_matrix[0][2] = P[i].distortion_tensorps[0][2] +
                         P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][2] +
                         P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][2] +
                         P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][2];
  product_matrix[1][0] = P[i].distortion_tensorps[1][0] +
                         P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][0] +
                         P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][0] +
                         P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][0];
  product_matrix[1][1] = P[i].distortion_tensorps[1][1] +
                         P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][1] +
                         P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][1] +
                         P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][1];
  product_matrix[1][2] = P[i].distortion_tensorps[1][2] +
                         P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][2] +
                         P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][2] +
                         P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][2];
  product_matrix[2][0] = P[i].distortion_tensorps[2][0] +
                         P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][0] +
                         P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][0] +
                         P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][0];
  product_matrix[2][1] = P[i].distortion_tensorps[2][1] +
                         P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][1] +
                         P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][1] +
                         P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][1];
  product_matrix[2][2] = P[i].distortion_tensorps[2][2] +
                         P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][2] +
                         P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][2] +
                         P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][2];

  /* this determinant will change sign when we pass through a caustic -> criterion for caustics */
  determinant = ((product_matrix[0][0]) * (product_matrix[1][1]) * (product_matrix[2][2]) +
		 (product_matrix[0][1]) * (product_matrix[1][2]) * (product_matrix[2][0]) +
		 (product_matrix[0][2]) * (product_matrix[1][0]) * (product_matrix[2][1]) -
		 (product_matrix[0][2]) * (product_matrix[1][1]) * (product_matrix[2][0]) -
		 (product_matrix[0][0]) * (product_matrix[1][2]) * (product_matrix[2][1]) -
		 (product_matrix[0][1]) * (product_matrix[1][0]) * (product_matrix[2][2]));

  /*
    Current and last NORMED stream density, linear order result of the last and current timestep.
  */
  current_normed_stream_density = 1.0/fabs(determinant);
  last_normed_stream_density    = 1.0/fabs(P[i].last_stream_determinant);

#ifdef ANNIHILATION_RADIATION
  /* avarage stream density */
  rho_0 = (current_normed_stream_density + last_normed_stream_density) / 2;

  /* extract phase-space information for annihilation radiation calculation -> cutoff density */
  P[i].rho_normed_cutoff_last    = P[i].rho_normed_cutoff_current;
  P[i].rho_normed_cutoff_current = analyse_phase_space(i, &s_1, &s_2, &s_3, &smear_x, &smear_y, &smear_z, &D_xP, &second_deriv_fac);
#endif

  /* CAUSTIC FOUND? -> was there a caustics between the last and the current timestep and does the particle type fit? */
  if((determinant * P[i].last_stream_determinant < 0.0) && (((1 << P[i].Type) & (CAUSTIC_FINDER))))
    {

#ifdef ANNIHILATION_RADIATION
      /* get normed_cutoff density */
      rho_normed_cutoff = P[i].rho_normed_cutoff_current;

      /*
       If rho_0 = (current_normed_stream_density + last_normed_stream_density) / 2 (so the average of current and last
       normed stream density) is below current normed cutoff density, we integrate analytically through the caustic
       for that time step.

       rho
       ^
       |         *  <--- rho_normed_cutoff
       |        ***
       |       *****
       |      *******   <-- area under curve ~ diverges logarithmicly
       |     *********
       |    ***********
       |   *************  <--- rho_0
       |---|-----|-----|----> t
         t_last t_c   t_current

       Note that the stream density is assumed to be symmetric around the caustic point: t_caustic = (t_last + t_current) / 2

       So we can center everything around t=0:

       rho
       ^
       |         *  <--- rho_normed_cutoff
       |        ***
       |       *****
       |      *******   <-- area under curve ~ diverges logarithmicly
       |     *********
       |    ***********
       |   *************  <--- rho_0
       |---|-----|-----|----> t
    -dt_drift/2  0   +dt_drift/2

       t_caustic = 0  and  dt_drift = t_current - t_last

       The fact that the numerical normed density near the caustic is below the cutoff, means that we cannot
       numerically resolve the very high caustic density within the time step.

       The analytic integration exploits the logarithmic divergence:

       Ansatz:
       rho(t) = kappa / |t|   --->   rho_0 = kappa / |dt_drift/2|   --->  kappa = dt_drift/2 * rho_0

       Integrated annihilation rate:

       2 \int_{-dt_drift/2}^{+dt_drift/2} dt rho(t) =
       2 - 2 kappa \int_{-dt_drift/2}^{0} dt 1/t  =  [ substitute drho / dt = + kappa/t^2 for t <0]
       2 2 kappa \int_{rho_0}^{rho_c} drho 1/rho =
       2 2 kappa ln(rho_c/rho_0)

       Therefore we need to add 2 kappa ln(\rho_c/\rho_0) when passing through an unresolved caustic,
       where rho_c is the calculated maximum density in the caustic.

       Note:
       We assume that the maximum density for the current and last time step do not differ much, so
       that we can compare rho_0 always to the current maximum density.

      */
#ifdef REINIT_AT_TURNAROUND
      /*
       We only integrate annihilation for particles that turned around already (otherwise wrong initial density)
       This is a security check, since each particle going through a caustic in general already turned around before.
      */
      if (P[i].turnaround_flag == 1)
      {
#endif
      /* if this condition is true we have to integrate analyitcally through the caustic */
      if (rho_0 < P[i].rho_normed_cutoff_current)
       {
        /* we do have an analytic integration through the caustic -> increase the counter and mark the flag */
        P[i].analytic_caustics          += 1.0;
        caustic_annihilation_flag        = 1;
        /* calculate analytic contribution */
        kappa_caustic                    = dt_drift / 2.0 * rho_0;
        analytic_add_to_annihilation     = 2.0 * kappa_caustic * log(fabs(P[i].rho_normed_cutoff_current/rho_0));

        /* add contribution to annihilation radiation */
        P[i].annihilation += analytic_add_to_annihilation;
       }
      summed_annihilation  = P[i].annihilation;
#ifdef REINIT_AT_TURNAROUND
      } /* end if (P[i].turnaround_flag == 1) */
#endif

#endif /* ANNIHILATION_RADIATION */

      /* increase caustic counter by one */
      P[i].caustic_counter += 1.0;
      caustic_counter       = P[i].caustic_counter;

#ifdef OUTPUT_LAST_CAUSTIC
      P[i].lc_Time    = All.Time;
      P[i].lc_Pos[0]  = P[i].Pos[0];
      P[i].lc_Pos[1]  = P[i].Pos[1];
      P[i].lc_Pos[2]  = P[i].Pos[2];
      P[i].lc_Vel[0]  = P[i].Vel[0];
      P[i].lc_Vel[1]  = P[i].Vel[1];
      P[i].lc_Vel[2]  = P[i].Vel[2];
#ifdef ANNIHILATION_RADIATION
      P[i].lc_rho_normed_cutoff = P[i].rho_normed_cutoff_current;
#else
      P[i].lc_rho_normed_cutoff = 0.0;
#endif
#endif

      init_density = P[i].init_density;

      /* write data to caustic log file -> allows caustics to be tracked on time-stepping frequency */
      fprintf(FdCaustics, "%g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
	      (MyOutputFloat) All.Time,
              (MyOutputFloat) determinant,
               P[i].ID,
	      (MyOutputFloat) P[i].Pos[0], (MyOutputFloat) P[i].Pos[1], (MyOutputFloat) P[i].Pos[2],
	      (MyOutputFloat) P[i].Vel[0], (MyOutputFloat) P[i].Vel[1], (MyOutputFloat) P[i].Vel[2],
              caustic_counter,
              s_1, s_2, s_3, smear_x, smear_y, smear_z,
              summed_annihilation,
              CurrentTurnaroundRadius,
              D_xP,
	      init_density / (1.0/(6.0*M_PI*All.G*All.Time*All.Time)),
              rho_normed_cutoff,
              second_deriv_fac,
              analytic_add_to_annihilation);
      fflush(FdCaustics);
    }

#ifdef ANNIHILATION_RADIATION
  /* is the numerical stream density above the normed_cutoff -> cut it */
  if (rho_0 > P[i].rho_normed_cutoff_current)
   {
    rho_0 = P[i].rho_normed_cutoff_current;
   }

  /* save cutted averaged stream density in physical units -> multiply with initial stream density */
  P[i].stream_density = P[i].init_density * rho_0;
#ifdef REINIT_AT_TURNAROUND
  /*
    We only integrate annihilation for particles that turned around already (otherwise wrong initial density)
    This is a security check, since each particle going through a caustic in general already turned around before.
  */
  if (P[i].turnaround_flag == 1)
   {
#endif
    /* integrate normed annihilation radiation -> multiplication with initial stream density is done in io.c to reduce round-off */
    if (caustic_annihilation_flag == 0)
      P[i].annihilation += rho_0 * dt_drift;
#ifdef REINIT_AT_TURNAROUND
   } /* end if (P[i].turnaround_flag == 1) */
#endif

#ifdef NO_CENTER_ANNIHILATION
  /* ignore the center for annihilation */
  if (sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]) < All.SofteningTable[P[i].Type])
   {
    /* set everything to zero in center */
    P[i].annihilation   = 0.0;
    P[i].stream_density = 0.0;
   }
#endif
#endif /*ANNIHILATION_RADIATION */

#ifdef REINIT_AT_TURNAROUND
 } /* end if (P[i].turnaround_flag >= 0) */
#endif

 /* update determinant, so we can identify the next caustic along the orbit */
 P[i].last_stream_determinant = determinant;

#endif /* DISTORTIONTENSORPS */
/* END PHASE-SPACE ANALYSIS ------------------------------------------------------ */

#ifndef HPM
  if(P[i].Type == 0)
    {
#ifdef PMGRID
      for(j = 0; j < 3; j++)
	SphP[i].VelPred[j] +=
	  (P[i].g.GravAccel[j] + P[i].GravPM[j]) * dt_gravkick + SphP[i].a.HydroAccel[j] * dt_hydrokick;
#else
      for(j = 0; j < 3; j++)
	SphP[i].VelPred[j] += P[i].g.GravAccel[j] * dt_gravkick + SphP[i].a.HydroAccel[j] * dt_hydrokick;
#endif

      SphP[i].d.Density *= exp(-SphP[i].v.DivVel * dt_drift);
      PPP[i].Hsml *= exp(0.333333333333 * SphP[i].v.DivVel * dt_drift);

      if(PPP[i].Hsml < All.MinGasHsml)
	PPP[i].Hsml = All.MinGasHsml;

#ifndef WAKEUP
      dt_step = (P[i].TimeBin ? (1 << P[i].TimeBin) : 0);
#else
      dt_step = P[i].dt_step;
#endif
      dt_entr = (time1 - (P[i].Ti_begstep + dt_step / 2)) * All.Timebase_interval;

#ifndef EOS_DEGENERATE
#ifndef MHM
#ifndef SOFTEREQS
#ifndef VORONOI_MESHRELAX

      SphP[i].Pressure = (SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density, GAMMA);

#endif
#else
      if(SphP[i].d.Density * a3inv >= All.PhysDensThresh)
	SphP[i].Pressure =
	  All.FactorForSofterEQS * (SphP[i].Entropy +
				    SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density,
									 GAMMA) + (1 -
										   All.
										   FactorForSofterEQS) *
	  GAMMA_MINUS1 * SphP[i].d.Density * All.InitGasU;
      else
	SphP[i].Pressure = (SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density, GAMMA);
#endif
#else
      /* Here we use an isothermal equation of state */
      SphP[i].Pressure = GAMMA_MINUS1 * SphP[i].d.Density * All.InitGasU;
      SphP[i].Entropy = SphP[i].Pressure / pow(SphP[i].d.Density, GAMMA);
#endif
#else
      /* call tabulated eos with physical units */
#ifdef EOS_ENERGY
      eos_calc_egiven2(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc,
		       SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr * All.UnitTime_in_s, &SphP[i].temp,
		       &SphP[i].Pressure);
      SphP[i].Pressure /= All.UnitPressure_in_cgs;
#else
      eos_calc_sgiven(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc,
		      SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr * All.UnitTime_in_s, &SphP[i].temp,
		      &SphP[i].Pressure, &SphP[i].u);
      SphP[i].Pressure /= All.UnitPressure_in_cgs;
#endif
#endif

#ifdef COSMIC_RAYS
#if defined( CR_UPDATE_PARANOIA )
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	CR_Particle_Update(SphP + i, CRpop);
#endif
#ifndef CR_NOPRESSURE
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	SphP[i].Pressure += CR_Comoving_Pressure(SphP + i, CRpop);
#endif
#endif

#if defined(MAGNETIC) && !defined(EULERPOTENTIALS)
      for(j = 0; j < 3; j++)
	SphP[i].BPred[j] += SphP[i].DtB[j] * dt_entr;
#ifdef DIVBCLEANING_DEDNER
      SphP[i].PhiPred += SphP[i].DtPhi * dt_entr;
#endif
#endif

    }
#endif /* end of HPM */

  P[i].Ti_current = time1;
}

void move_particles(int time1)
{
  int i;

  if(ThisTask == 0)
    printf("MOVE\n");

  for(i = 0; i < NumPart; i++)
    drift_particle(i, time1);
}

/*! This function makes sure that all particle coordinates (Pos) are
 *  periodically mapped onto the interval [0, BoxSize].  After this function
 *  has been called, a new domain decomposition should be done, which will
 *  also force a new tree construction.
 */
#ifdef PERIODIC
void do_box_wrapping(void)
{
  int i, j;
  double boxsize[3];

  for(j = 0; j < 3; j++)
    boxsize[j] = All.BoxSize;

#ifdef LONG_X
  boxsize[0] *= LONG_X;
#endif
#ifdef LONG_Y
  boxsize[1] *= LONG_Y;
#endif
#ifdef LONG_Z
  boxsize[2] *= LONG_Z;
#endif

  for(i = 0; i < NumPart; i++)
    for(j = 0; j < 3; j++)
      {
	while(P[i].Pos[j] < 0)
	  P[i].Pos[j] += boxsize[j];

	while(P[i].Pos[j] >= boxsize[j])
	  P[i].Pos[j] -= boxsize[j];
      }
}
#endif

/*

#ifdef XXLINFO
#ifdef MAGNETIC
  double MeanB_part = 0, MeanB_sum;

#ifdef TRACEDIVB
  double MaxDivB_part = 0, MaxDivB_all;
  double dmax1, dmax2;
#endif
#endif
#ifdef TIME_DEP_ART_VISC
  double MeanAlpha_part = 0, MeanAlpha_sum;
#endif
#endif

#ifdef XXLINFO
      if(Flag_FullStep == 1)
        {
#ifdef MAGNETIC
          MeanB_part += sqrt(SphP[i].BPred[0] * SphP[i].BPred[0] +
                             SphP[i].BPred[1] * SphP[i].BPred[1] + SphP[i].BPred[2] * SphP[i].BPred[2]);
#ifdef TRACEDIVB
          MaxDivB_part = DMAX(MaxDivB, fabs(SphP[i].divB));
#endif
#endif
#ifdef TIME_DEP_ART_VISC
          MeanAlpha_part += SphP[i].alpha;
#endif
        }
#endif

#ifdef XXLINFO
  if(Flag_FullStep == 1)
    {
#ifdef MAGNETIC
      MPI_Reduce(&MeanB_part, &MeanB_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      if(ThisTask == 0)
	MeanB = MeanB_sum / All.TotN_gas;
#ifdef TRACEDIVB
      MPI_Reduce(&MaxDivB_part, &MaxDivB_all, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
      if(ThisTask == 0)
	MaxDivB = MaxDivB_all;
#endif
#endif
#ifdef TIME_DEP_ART_VISC
      MPI_Reduce(&MeanAlpha_part, &MeanAlpha_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      if(ThisTask == 0)
	MeanAlpha = MeanAlpha_sum / All.TotN_gas;
#endif
    }
#endif
*/

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
// rand.h
#ifndef RAND_H
#define RAND_H
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define URNG generator

class generator{
    const gsl_rng_type * T;
    public:
    gsl_rng * r;
    generator();
    ~generator();
    void seed(unsigned long int s);
};

#include "types.h"

class r_lognorm{
    double zeta;
    double sigma;
    public:
    r_lognorm(REAL mean,REAL sigma);
    ~r_lognorm();
    REAL operator()(URNG& g);
};

REAL lognorm_func(URNG& g,REAL mean,REAL sig);

class r_hypgeo{
    unsigned int n1;
    unsigned int n2;
    unsigned int t;
    public:
    r_hypgeo(unsigned int population1,unsigned int population2,unsigned int samples);
    ~r_hypgeo();
    unsigned int  operator()(URNG& g);
};

unsigned int hypgeo_func(URNG& g,unsigned int population1,unsigned int population2,unsigned int samples);

class r_uni{
    double a;
    double b;
    public:
    r_uni(REAL min,REAL max);
    ~r_uni();
    REAL operator()(URNG& g);
};

REAL unif_func(URNG& g,REAL min,REAL max);

REAL exp_func(URNG& g,REAL mu);

int ber_func(URNG& g ,REAL p);

#endif

ext: h


text:
/***********************************
V2 IMPLEMENTATION USING DATASETS
************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <cblas.h>
#include <string.h>
#include <mpi.h>

#include "../inc/knnring.h"
#include "../inc/readlib.h"
#include "../inc/vptring.h"

#define B 3

int main(int argc,char *argv[]){

  int p_rank;
  int p_size;
  int n;
  int d;
  int k = 10;
  double *X;

  //CHOOSE WHICH DATASET
  char file_path[200] = "../data/ColorHistogram.asc";

  //INITIALIZE MPI COMMUNICATIONS
  //MPI_Init(&argc,&argv);   doesn't need any argument

  MPI_Init(NULL,NULL);

  MPI_Comm_size(MPI_COMM_WORLD, &p_size);
  MPI_Comm_rank(MPI_COMM_WORLD, &p_rank);
  MPI_Status status;                      // STATUS COMMUNICATION
  int tag = 1;                            // COMMUNICATION TAG

  //time variables to measure time
  struct timespec ts_start;
  struct timespec ts_end;
  double time;

  if(p_rank == 0){

    //read all the data
    double *X_total;
    int n_total;
    X_total = read_data(file_path,&n_total,&d);

    //send n,d variables and the points coordinates to other processes.
    n = n_total/p_size;
    X = (double *)malloc(n*d*sizeof(double));
    for(int i=1;i<p_size;i++){
      for(int j=0;j<n*d;j++){
        X[j] = X_total[i*n*d+j];
      }

      //send to every PROCESS
      MPI_Send(&n,1, MPI_INT,    i, tag, MPI_COMM_WORLD);
      MPI_Send(&d,1, MPI_INT,    i, tag, MPI_COMM_WORLD);
      MPI_Send(X,n*d,MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
    }

    //keep the first n points for myself.
    for(int j=0;j<n*d;j++){
      X[j] = X_total[j];
    }
    free(X_total);

    clock_gettime(CLOCK_MONOTONIC, &ts_start);
    knnresult results = distrV2(X,n,d,k,B);
    clock_gettime(CLOCK_MONOTONIC, &ts_end);
    time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
    printf("TOTAL V2 time: %lf\n",time);

    //WRITE TO FILE
    char str[200];
    snprintf(str,sizeof(str),"results/dist_v2_%d_%d_%d.txt",n,d,k);
    write_to_file(str,results.nidx,results.ndist,n*k,results.k);

    for(int p=1;p<p_size;p++){
      MPI_Recv(results.nidx, n*k,    MPI_INT,p,tag,MPI_COMM_WORLD, &status);
      MPI_Recv(results.ndist,n*k, MPI_DOUBLE,p,tag,MPI_COMM_WORLD, &status);
      write_to_file(str,results.nidx,results.ndist,n*k,results.k);
    }
  }
  else{
    //receive number of points n and distances.
    MPI_Recv(&n,1,MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
    MPI_Recv(&d,1,MPI_INT, 0, tag, MPI_COMM_WORLD, &status);

    //receive coordinates
    X = (double *)malloc(n*d*sizeof(double));
    MPI_Recv(X,n*d,MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);

    clock_gettime(CLOCK_MONOTONIC, &ts_start);
    knnresult results = distrV2(X,n,d,k,B);
    clock_gettime(CLOCK_MONOTONIC, &ts_end);
    time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
    //printf("TOTAL V2 time: %lf from process %d\n",time,p_rank);

    //send back the results
    MPI_Send(results.nidx, n*k, MPI_INT   ,0,tag,MPI_COMM_WORLD);
    MPI_Send(results.ndist,n*k, MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
  }
  //execution is over
  MPI_Finalize();
  printf("DONE ....=>EXITING\n");
  return 0;
}


text:
/* -*- linux-c -*- */
/* fewbody.h

   Copyright (C) 2002-2004 John M. Fregeau

   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 2 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#ifndef _FEWBODY_H
#define _FEWBODY_H 1

#include <stdio.h>
#include <gsl/gsl_nan.h>
#include <gsl/gsl_rng.h>

/* version information */
#define FB_VERSION "0.26"
#define FB_NICK "Oblivion"
#define FB_DATE "Mon Sep  6 09:05:47 EDT 2010"

/* dimensionless constants */
#define FB_CONST_PI 3.141592653589793238462643

/* constants, in cgs units */
#define FB_CONST_MSUN 1.989e+33
#define FB_CONST_RSUN 6.9599e+10
#define FB_CONST_C 2.99792458e+10
#define FB_CONST_G 6.67259e-8
#define FB_CONST_AU 1.496e+13
#define FB_CONST_PARSEC 3.0857e+18
#define FB_CONST_YR 3.155693e+7

/* these usually shouldn't need to be changed */
#define FB_H 1.0e-2
#define FB_SSTOP GSL_POSINF
#define FB_AMIN GSL_POSINF
#define FB_RMIN GSL_POSINF
#define FB_ROOTSOLVER_MAX_ITER 100
#define FB_ROOTSOLVER_ABS_ACC 1.0e-11
#define FB_ROOTSOLVER_REL_ACC 1.0e-11
#define FB_MAX_STRING_LENGTH 2048
#define FB_MAX_LOGENTRY_LENGTH (32 * FB_MAX_STRING_LENGTH)

/* a struct containing the units used */
typedef struct{
	double v; /* velocity */
	double l; /* length */
	double t; /* time */
	double m; /* mass */
	double E; /* energy */
} fb_units_t;

/* the fundamental object */
typedef struct fb_obj{
	int ncoll; /* total number of stars collided together in this star */
	long *id; /* numeric id array */
	char idstring[FB_MAX_STRING_LENGTH]; /* string id */
	double m; /* mass */
	double R; /* radius */
	double Eint; /* internal energy (used to check energy conservation) */
	double Lint[3]; /* internal ang mom (used to check ang mom conservation) */
	double x[3]; /* position */
	double v[3]; /* velocity */
	int n; /* total number of stars in hierarchy */
	struct fb_obj *obj[2]; /* pointers to children */
	double a; /* semimajor axis */
	double e; /* eccentricity */
	double Lhat[3]; /* angular momentum vector */
	double Ahat[3]; /* Runge-Lenz vector */
	double t; /* time at which node was upsynced */
	double mean_anom; /* mean anomaly when node was upsynced */
} fb_obj_t;

/* parameters for the K-S integrator */
typedef struct{
	int nstar; /* number of actual stars */
	int kstar; /* nstar*(nstar-1)/2, number of separations */
	double *m; /* m[nstar] */
	double *M; /* M[kstar] */
	double **amat; /* amat[nstar][kstar] */
	double **Tmat; /* Tmat[kstar][kstar] */
	double Einit; /* initial energy used in integration scheme */
} fb_ks_params_t;

/* parameters for the non-regularized integrator */
typedef struct{
	int nstar; /* number of actual stars */
	double *m; /* m[nstar] */
} fb_nonks_params_t;

/* the hierarchy data structure */
typedef struct{
	int nstarinit; /* initial number of stars (may not equal nstar if there are collisions) */
	int nstar; /* number of stars */
	int nobj; /* number of binary trees */
	int *hi; /* hierarchical index array */
	int *narr; /* narr[i] = number of hierarchical objects with i elements */
	fb_obj_t *hier; /* memory location of hierarchy information */
	fb_obj_t **obj; /* array of pointers to top nodes of binary trees */
} fb_hier_t;

/* input parameters */
typedef struct{
	int ks; /* 0=no regularization, 1=K-S regularization */
	double tstop; /* stopping time, in units of t_dyn */
	int Dflag; /* 0=don't print to stdout, 1=print to stdout */
	double dt; /* time interval between printouts will always be greater than this value */
	double tcpustop; /* cpu stopping time, in units of seconds */
	double absacc; /* absolute accuracy of the integrator */
	double relacc; /* relative accuracy of the integrator */
	int ncount; /* number of integration steps between each call to fb_classify() */
	double tidaltol; /* tidal tolerance */
	char firstlogentry[FB_MAX_LOGENTRY_LENGTH]; /* first entry to put in printout log */
	double fexp; /* expansion factor for a merger product: R = f_exp (R_1+R_2) */
} fb_input_t;

/* return parameters */
typedef struct{
	long count; /* number of integration steps */
	int retval; /* return value; 1=success, 0=failure */
	long iclassify; /* number of times classify was called */
	double tcpu; /* cpu time taken */
	double DeltaE; /* change in energy */
	double DeltaEfrac; /* change in energy, as a fraction of initial energy */
	double DeltaL; /* change in ang. mom. */
	double DeltaLfrac; /* change in ang. mom., as a fraction of initial ang. mom. */
	double Rmin; /* minimum distance of close approach during interaction */
	int Rmin_i; /* index of star i participating in minimum close approach */
	int Rmin_j; /* index of star j participating in minimum close approach */
	int Nosc; /* number of oscillations of the quantity s^2 (McMillan & Hut 1996) (Nosc=Nmin-1, so resonance if Nosc>=1) */
} fb_ret_t;

/* fewbody.c */
fb_ret_t fewbody(fb_input_t input, fb_hier_t *hier, double *t);

/* fewbody_classify.c */
int fb_classify(fb_hier_t *hier, double t, double tidaltol);
int fb_is_stable(fb_obj_t *obj);
int fb_is_stable_binary(fb_obj_t *obj);
int fb_is_stable_triple(fb_obj_t *obj);
int fb_is_stable_quad(fb_obj_t *obj);
int fb_mardling(fb_obj_t *obj, int ib, int is);

/* fewbody_coll.c */
int fb_is_collision(double r, double R1, double R2);
int fb_collide(fb_hier_t *hier, double f_exp);
void fb_merge(fb_obj_t *obj1, fb_obj_t *obj2, int nstarinit, double f_exp);

/* fewbody_hier.c */
void fb_malloc_hier(fb_hier_t *hier);
void fb_init_hier(fb_hier_t *hier);
void fb_free_hier(fb_hier_t hier);
void fb_trickle(fb_hier_t *hier, double t);
void fb_elkcirt(fb_hier_t *hier, double t);
int fb_create_indices(int *hi, int nstar);
int fb_n_hier(fb_obj_t *obj);
char *fb_sprint_hier(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH]);
char *fb_sprint_hier_hr(fb_hier_t hier, char string[FB_MAX_STRING_LENGTH]);
void fb_upsync(fb_obj_t *obj, double t);
void fb_randorient(fb_obj_t *obj, gsl_rng *rng);
void fb_downsync(fb_obj_t *obj, double t);
void fb_objcpy(fb_obj_t *obj1, fb_obj_t *obj2);

/* fewbody_int.c */
void fb_malloc_ks_params(fb_ks_params_t *ks_params);
void fb_init_ks_params(fb_ks_params_t *ks_params, fb_hier_t hier);
void fb_free_ks_params(fb_ks_params_t ks_params);
void fb_malloc_nonks_params(fb_nonks_params_t *nonks_params);
void fb_init_nonks_params(fb_nonks_params_t *nonks_params, fb_hier_t hier);
void fb_free_nonks_params(fb_nonks_params_t nonks_params);

/* fewbody_io.c */
void fb_print_version(FILE *stream);
void fb_print_story(fb_obj_t *star, int nstar, double t, char *logentry);

/* fewbody_isolate.c */
int fb_collapse(fb_hier_t *hier, double t, double tidaltol);
int fb_expand(fb_hier_t *hier, double t, double tidaltol);

/* fewbody_ks.c */
double fb_ks_dot(double x[4], double y[4]);
double fb_ks_mod(double x[4]);
void fb_calc_Q(double q[4], double Q[4]);
void fb_calc_ksmat(double Q[4], double Qmat[4][4]);
void fb_calc_amat(double **a, int nstar, int kstar);
void fb_calc_Tmat(double **a, double *m, double **T, int nstar, int kstar);
int fb_ks_func(double s, const double *y, double *f, void *params);
double fb_ks_Einit(const double *y, fb_ks_params_t params);
void fb_euclidean_to_ks(fb_obj_t **star, double *y, int nstar, int kstar);
void fb_ks_to_euclidean(double *y, fb_obj_t **star, int nstar, int kstar);

/* fewbody_nonks.c */
int fb_nonks_func(double t, const double *y, double *f, void *params);
int fb_nonks_jac(double t, const double *y, double *dfdy, double *dfdt, void *params);
void fb_euclidean_to_nonks(fb_obj_t **star, double *y, int nstar);
void fb_nonks_to_euclidean(double *y, fb_obj_t **star, int nstar);

/* fewbody_scat.c */
void fb_init_scattering(fb_obj_t *obj0, fb_obj_t *obj1, double vinf, double b, double rtid);
void fb_normalize(fb_hier_t *hier, fb_units_t units);

/* fewbody_utils.c */
double *fb_malloc_vector(int n);
double **fb_malloc_matrix(int nr, int nc);
void fb_free_vector(double *v);
void fb_free_matrix(double **m);
double fb_sqr(double x);
double fb_cub(double x);
double fb_dot(double x[3], double y[3]);
double fb_mod(double x[3]);
int fb_cross(double x[3], double y[3], double z[3]);
int fb_angmom(fb_obj_t *star, int nstar, double L[3]);
void fb_angmomint(fb_obj_t *star, int nstar, double L[3]);
double fb_einttot(fb_obj_t *star, int nstar);
double fb_petot(fb_obj_t *star, int nstar);
double fb_ketot(fb_obj_t *star, int nstar);
double fb_outerpetot(fb_obj_t **obj, int nobj);
double fb_outerketot(fb_obj_t **obj, int nobj);
double fb_kepler(double e, double mean_anom);
double fb_keplerfunc(double mean_anom, void *params);
double fb_reltide(fb_obj_t *bin, fb_obj_t *single, double r);

/* fewbody_ui.c */
int fbui_new_hier(fb_hier_t *hier, int n);
int fbui_delete_hier(fb_hier_t *hier);
fb_obj_t *fbui_hierarchy_element(fb_hier_t *hier, int n, int m);
fb_obj_t *fbui_hierarchy_single(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_binary(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_triple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_quadruple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_quintuple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_sextuple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_septuple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_octuple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_nonuple(fb_hier_t *hier, int m);
fb_obj_t *fbui_hierarchy_decuple(fb_hier_t *hier, int m);
int fbui_make_pair(fb_obj_t *parentobj, fb_obj_t *obj1, fb_obj_t *obj2);
int fbui_initialize_single(fb_obj_t *single, long id, char idstring[FB_MAX_STRING_LENGTH]);
fb_obj_t *fbui_tree(fb_hier_t *hier, int n);
int fbui_obj_ncoll_set(fb_obj_t *obj, int ncoll);
int fbui_obj_ncoll_get(fb_obj_t *obj);
int fbui_obj_id_set(fb_obj_t *obj, int index, long id);
long fbui_obj_id_get(fb_obj_t *obj, int index);
int fbui_obj_idstring_set(fb_obj_t *obj, char idstring[FB_MAX_STRING_LENGTH]);
char *fbui_obj_idstring_get(fb_obj_t *obj);
int fbui_obj_mass_set(fb_obj_t *obj, double mass);
double fbui_obj_mass_get(fb_obj_t *obj);
int fbui_obj_radius_set(fb_obj_t *obj, double radius);
double fbui_obj_radius_get(fb_obj_t *obj);
int fbui_obj_Eint_set(fb_obj_t *obj, double Eint);
double fbui_obj_Eint_get(fb_obj_t *obj);
int fbui_obj_Lint_set(fb_obj_t *obj, double Lint[3]);
int fbui_obj_Linti_set(fb_obj_t *obj, double Lint0, double Lint1, double Lint2);
double *fbui_obj_Lint_get(fb_obj_t *obj);
int fbui_obj_x_set(fb_obj_t *obj, double x[3]);
int fbui_obj_xi_set(fb_obj_t *obj, double x0, double x1, double x2);
double *fbui_obj_x_get(fb_obj_t *obj);
int fbui_obj_v_set(fb_obj_t *obj, double v[3]);
int fbui_obj_vi_set(fb_obj_t *obj, double v0, double v1, double v2);
double *fbui_obj_v_get(fb_obj_t *obj);
int fbui_obj_n_set(fb_obj_t *obj, int n);
int fbui_obj_n_get(fb_obj_t *obj);
int fbui_obj_obj_set(fb_obj_t *obj, int i, fb_obj_t *objtopointto);
int fbui_obj_left_child_set(fb_obj_t *obj, fb_obj_t *objtopointto);
int fbui_obj_right_child_set(fb_obj_t *obj, fb_obj_t *objtopointto);
fb_obj_t *fbui_obj_obj_get(fb_obj_t *obj, int i);
fb_obj_t *fbui_obj_left_child_get(fb_obj_t *obj);
fb_obj_t *fbui_obj_right_child_get(fb_obj_t *obj);
int fbui_obj_a_set(fb_obj_t *obj, double a);
double fbui_obj_a_get(fb_obj_t *obj);
int fbui_obj_e_set(fb_obj_t *obj, double e);
double fbui_obj_e_get(fb_obj_t *obj);
int fbui_obj_Lhat_set(fb_obj_t *obj, double Lhat[3]);
int fbui_obj_Lhati_set(fb_obj_t *obj, double Lhat0, double Lhat1, double Lhat2);
double *fbui_obj_Lhat_get(fb_obj_t *obj);
int fbui_obj_Ahat_set(fb_obj_t *obj, double Ahat[3]);
int fbui_obj_Ahati_set(fb_obj_t *obj, double Ahat0, double Ahat1, double Ahat2);
double *fbui_obj_Ahat_get(fb_obj_t *obj);
int fbui_obj_t_set(fb_obj_t *obj, double t);
double fbui_obj_t_get(fb_obj_t *obj);
int fbui_obj_mean_anom_set(fb_obj_t *obj, double mean_anom);
double fbui_obj_mean_anom_get(fb_obj_t *obj);

/* macros */
/* The variadic macro syntax here conforms to the C99 standard, but for some
   reason won't compile on Mac OSX with gcc. */
/* #define fb_dprintf(...) if (fb_debug) fprintf(stderr, __VA_ARGS__) */
/* The variadic macro syntax here is the old gcc standard, and compiles on
   Mac OSX with gcc. */
#define fb_dprintf(args...) if (fb_debug) fprintf(stderr, args)
#define FB_MIN(a, b) ((a)<=(b)?(a):(b))
#define FB_MAX(a, b) ((a)>=(b)?(a):(b))
#define FB_DELTA(i, j) ((i)==(j)?1:0)
#define FB_KS_K(i, j, nstar) ((i)*(nstar)-((i)+1)*((i)+2)/2+(j))

/* there is just one global variable */
extern int fb_debug;

#endif /* fewbody.h */

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/* Copyright (C) 2010-2021 Barcelona Supercomputing Center and University of
 * Illinois at Urbana-Champaign
 * SPDX-License-Identifier: MIT
 */

/** \file
 * \brief Wrapper routines for GSL random number functions.
 */

/* clang-format off */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <stdio.h>
#include <time.h>

/** \brief Private internal-use variable to store the random number
 * generator.
 */
static gsl_rng *camp_rand_gsl_rng = NULL;

/** \brief Result code indicating successful completion.
 */
#define CAMP_RAND_GSL_SUCCESS      0
/** \brief Result code indicating initialization failure.
 */
#define CAMP_RAND_GSL_INIT_FAIL    1
/** \brief Result code indicating the generator was not initialized
 * when it should have been.
 */
#define CAMP_RAND_GSL_NOT_INIT     2
/** \brief Result code indicating the generator was already
 * initialized when an initialization was attempted.
 */
#define CAMP_RAND_GSL_ALREADY_INIT 3

/** \brief Initialize the random number generator with the given seed.
 *
 * This must be called before any other GSL random number functions
 * are called.
 *
 * \param seed The random seed to use.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 * \sa camp_rand_finalize_gsl() to cleanup the generator.
 */
int camp_srand_gsl(int seed)
{
        if (camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_ALREADY_INIT;
        }
        gsl_set_error_handler_off(); // turn off automatic error handling
        camp_rand_gsl_rng = gsl_rng_alloc(gsl_rng_mt19937);
        if (camp_rand_gsl_rng == NULL) {
                return CAMP_RAND_GSL_INIT_FAIL;
        }
        gsl_rng_set(camp_rand_gsl_rng, seed);
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Cleanup and deallocate the random number generator.
 *
 * This must be called after camp_srand_gsl().
 *
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_finalize_gsl()
{
        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        gsl_rng_free(camp_rand_gsl_rng);
        camp_rand_gsl_rng = NULL;
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Generate a uniform random number in \f$[0,1)\f$.
 *
 * \param harvest A pointer to the generated random number.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_gsl(double *harvest)
{
        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        *harvest = gsl_rng_uniform(camp_rand_gsl_rng);
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Generate a uniform random integer in \f$[1,n]\f$.
 *
 * \param n The upper limit of the random integer.
 * \param harvest A pointer to the generated random number.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_int_gsl(int n, int *harvest)
{
        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        *harvest = gsl_rng_uniform_int(camp_rand_gsl_rng, n) + 1;
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Generate a normally-distributed random number.
 *
 * \param mean The mean of the distribution.
 * \param stddev The standard deviation of the distribution.
 * \param harvest A pointer to the generated random number.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_normal_gsl(double mean, double stddev, double *harvest)
{
        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        *harvest = gsl_ran_gaussian(camp_rand_gsl_rng, stddev) + mean;
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Generate a Poisson-distributed random integer.
 *
 * \param mean The mean of the distribution.
 * \param harvest A pointer to the generated random number.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_poisson_gsl(double mean, int *harvest)
{
        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        *harvest = gsl_ran_poisson(camp_rand_gsl_rng, mean);
        return CAMP_RAND_GSL_SUCCESS;
}

/** \brief Generate a Binomial-distributed random integer.
 *
 * \param n The sample size for the distribution.
 * \param p The sample probability for the distribution.
 * \param harvest A pointer to the generated random number.
 * \return CAMP_RAND_GSL_SUCCESS on success, otherwise an error code.
 */
int camp_rand_binomial_gsl(int n, double p, int *harvest)
{
        unsigned int u;

        if (!camp_rand_gsl_rng) {
                return CAMP_RAND_GSL_NOT_INIT;
        }
        u = n;
        *harvest = gsl_ran_binomial(camp_rand_gsl_rng, p, u);
        return CAMP_RAND_GSL_SUCCESS;
}

/* clang-format on */

max_stars_repo_name: open-atmos/camp

max_issues_repo_name: open-atmos/camp

max_forks_repo_name: open-atmos/camp


text:
/*
   For more information, please see: http://software.sci.utah.edu

   The MIT License

   Copyright (c) 2015 Scientific Computing and Imaging Institute,
   University of Utah.

   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.
*/

///@brief
///@file  DenseMatrixMultiplication.h

#ifndef CORE_DATATYPES_DENSEMATRIXMULTIPLICATION_H
#define CORE_DATATYPES_DENSEMATRIXMULTIPLICATION_H

#if defined(HAVE_CBLAS)
#if defined(__APPLE__)
#include <vecLib/cblas.h>
#else
extern "C"{
#include <cblas.h>
}
#endif
#endif

namespace SCIRun {

template <typename T>
void
DenseMatrixGeneric<T>::mult(const ColumnMatrix& x, ColumnMatrix& b,
                            index_type beg, index_type end,
                            int spVec) const
{
  ASSERTEQ(x.nrows(), this->ncols_);
  ASSERTEQ(b.nrows(), this->nrows_);

  if (beg == -1) beg = 0;
  if (end == -1) end = this->nrows_;

#if defined(HAVE_CBLAS)
  double ALPHA = 1.0;
  double BETA = 0.0;
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, (end-beg),
    1, this->ncols_, ALPHA, dataptr_+(beg*this->ncols_), this->ncols_,
    x.get_data_pointer(), 1, BETA,
    b.get_data_pointer()+beg, 1);
#else

  double* xdata = x.get_data_pointer();
  double* bdata = b.get_data_pointer();

  size_type m8 = (this->ncols_)/8;
  size_type m = (this->ncols_)-m8*8;

  index_type i, j;
  if(!spVec)
  {
    for (i=beg; i<end; i++)
    {
      double sum=0;
      double* row=data[i];
      double* xd=xdata;

      for (j=0; j<m8; j++)
      {
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
        sum+=(*row)*(*xd);
        row++; xd++;
      }

      for (j=0; j<m; j++)
      {
        sum+=(*row)*(*xd);
        row++; xd++;
      }
      (*bdata)=sum;
      bdata++;
    }
  }
  else
  {
    for (i=beg; i<end; i++) b[i]=0;
    for (j=0; j<this->ncols_; j++)
    {
      if (x[j])
        for (i=beg; i<end; i++)
        {
          b[i]+=data[i][j]*x[j];
        }
    }
  }

#endif
}

template <typename T>
void
DenseMatrixGeneric<T>::multiply(ColumnMatrix& x, ColumnMatrix& b) const
{
  index_type i, j;

  double* xdata = x.get_data_pointer();
  double* bdata = b.get_data_pointer();

  size_type m8 = (this->ncols_)/8;
  size_type m = (this->ncols_)-m8*8;

  for (i=0; i<this->nrows_; i++)
  {
    double sum=0;
    double* row = data[i];
    double* xd = xdata;
    for (j=0; j<m8; j++)
    {
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
      sum+=(*row)*(*xd);
      row++; xd++;
    }

    for (j=0; j<m; j++)
    {
      sum+=(*row)*(*xd);
      row++; xd++;
    }

    *bdata=sum;
    bdata++;
  }
}

template <typename T>
void
DenseMatrixGeneric<T>::mult_transpose(const ColumnMatrix& x, ColumnMatrix& b,
                                      index_type beg, index_type end,
                                      int spVec) const
{
  // Compute At*x=b
  ASSERT(x.nrows() == this->nrows_);
  ASSERT(b.nrows() == this->ncols_);
  if (beg == -1) beg = 0;
  if (end == -1) end = this->ncols_;
  index_type i, j;
  if (!spVec)
  {
    for (i=beg; i<end; i++)
    {
      double sum=0;
      for (j=0; j<this->nrows_; j++)
      {
        sum+=data[j][i]*x[j];
      }
      b[i]=sum;
    }
  }
  else
  {
    for (i=beg; i<end; i++) b[i]=0;
    for (j=0; j<this->nrows_; j++)
      if (x[j])
      {
        double *row=data[j];
        for (i=beg; i<end; i++)
          b[i]+=row[i]*x[j];
      }
  }
}

template <typename T>
DenseMatrix*
DenseMatrixGeneric<T>::make_diagonal_from_column(const ColumnMatrix& column, size_type rows, size_type cols)
{
  DenseMatrix* result = zero_matrix(rows, cols);
  for (size_type i = 0; i < column.nrows(); ++i)
    (*result)[i][i] = column[i];

  return result;
}

} // End namespace SCIRun

#endif


text:
#pragma once

#include <gsl/gsl>
#include <codecvt>
#include <locale>
#include <string>

#ifdef _MSC_VER
#include <string_view>
#else
// llvm-libc++ included with Android NDK r15c hasn't yet promoted the C++1z library fundamentals
// technical specification to the final C++17 specification for string_view.
// Force promote the types we need up into the std namespace.
#include <experimental/string_view>
namespace std
{
    using string_view = experimental::string_view;
    using wstring_view = experimental::wstring_view;
}
#endif

namespace arcana
{
    inline std::string utf16_to_utf8(gsl::cwzstring<> input)
    {
        std::wstring_convert<std::codecvt_utf8_utf16<wchar_t>> converter;
        return converter.to_bytes(input);
    }

    inline std::wstring utf8_to_utf16(gsl::czstring<> input)
    {
        std::wstring_convert<std::codecvt_utf8<wchar_t>> converter;
        return converter.from_bytes(input);
    }

    inline std::string utf16_to_utf8(std::wstring_view input)
    {
        std::wstring_convert<std::codecvt_utf8_utf16<wchar_t>> converter;
        return converter.to_bytes(input.data(), input.data() + input.size());
    }

    inline std::wstring utf8_to_utf16(std::string_view input)
    {
        std::wstring_convert<std::codecvt_utf8_utf16<wchar_t>> converter;
        return converter.from_bytes(input.data(), input.data() + input.size());
    }

    struct string_compare
    {
        using is_transparent = std::true_type;

        bool operator()(gsl::cstring_span<> a, gsl::cstring_span<> b) const
        {
            return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());
        }

        bool operator()(gsl::czstring_span<> a, gsl::czstring_span<> b) const
        {
            return (*this)(a.as_string_span(), b.as_string_span());
        }

        bool operator()(gsl::czstring<> a, gsl::czstring<> b) const
        {
            return strcmp(a, b) < 0;
        }
    };
}