text:
#pragma once

#include <memory>
#include <gsl\gsl>
#include "RTTI.h"

namespace Library
{
	class Game;

	class AbstractContentTypeReader : public RTTI
	{
		RTTI_DECLARATIONS(AbstractContentTypeReader, RTTI)

	public:
		AbstractContentTypeReader(const AbstractContentTypeReader&) = default;
		AbstractContentTypeReader& operator=(const AbstractContentTypeReader&) = default;
		AbstractContentTypeReader(AbstractContentTypeReader&&) = default;
		AbstractContentTypeReader& operator=(AbstractContentTypeReader&&) = default;
		virtual ~AbstractContentTypeReader() = default;

		std::uint64_t TargetTypeId() const;
		virtual std::shared_ptr<RTTI> Read(const std::wstring& assetName) = 0;

	protected:
		AbstractContentTypeReader(Game& game, const std::uint64_t targetTypeId);

		gsl::not_null<Game*> mGame;
		const std::uint64_t mTargetTypeId;
	};

	template <typename T>
	class ContentTypeReader : public AbstractContentTypeReader
	{
	public:
		ContentTypeReader(const ContentTypeReader&) = default;
		ContentTypeReader& operator=(const ContentTypeReader&) = default;
		ContentTypeReader(ContentTypeReader&&) = default;
		ContentTypeReader& operator=(ContentTypeReader&&) = default;
		virtual ~ContentTypeReader() = default;

		virtual std::shared_ptr<RTTI> Read(const std::wstring& assetName) override;

	protected:
		ContentTypeReader(Game& game, const std::uint64_t targetTypeId);

		virtual std::shared_ptr<T> _Read(const std::wstring& assetName) = 0;
	};
}

#include "ContentTypeReader.inl"


text:
#include <cblas.h>

#include "tasks.h"
#include "utils.h"

void gemm_task_par_reconfigure(int nth)
{
    // empty
}

void gemm_task_par_finalize(void)
{
    // empty
}

void gemm_task_par(void *ptr, int nth, int me)
{
    struct gemm_task_arg *arg = (struct gemm_task_arg*) ptr;

    int m        = arg->m;
    int n        = arg->n;
    int k        = arg->k;
    double *A21  = arg->A21;
    double *A21T = arg->A21T;
    double *A22  = arg->A22;
    int ldA      = arg->ldA;

    // Cut A21T into panels of size m-by-blksz. Compute nominal block size.
    int blksz = iceil(n, nth);

    // Determine my share of A21T.
    int my_first_col = blksz * me;
    int my_num_cols  = min(blksz, n - my_first_col);

    // Compute A22 := A22 - A21 * A21T', using my share of A22 and A21T.
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
                m, my_num_cols, k,
                -1.0, A21,                      ldA,
                      A21T + my_first_col,      ldA,
                 1.0, A22 + my_first_col * ldA, ldA);
}

ext: c

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */

/* This file is part of phonopy. */

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

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

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

/* * Neither the name of the phonopy project 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 <lapacke.h>
#include <stdlib.h>
#include <phonoc_array.h>
#include <phonoc_utils.h>
#include <phonon4_h/frequency_shift.h>
#include <phonon4_h/real_to_reciprocal.h>

static void get_fc4_normal_for_frequency_shift_at_gp
(double *fc4_normal_real,
 const double *frequencies,
 const lapack_complex_double *eigenvectors,
 const int grid_point0,
 const int grid_point1,
 const int *grid_address,
 const int *mesh,
 const double *fc4,
 const Darray *shortest_vectors,
 const Iarray *multiplicity,
 const double *masses,
 const int *p2s_map,
 const int *s2p_map,
 const Iarray *band_indices,
 const double cutoff_frequency);
static lapack_complex_double fc4_sum(const int bi0,
				     const int bi1,
				     const lapack_complex_double *eigvecs0,
				     const lapack_complex_double *eigvecs1,
				     const lapack_complex_double *fc4_reciprocal,
				     const double *masses,
				     const int num_atom);
static int collect_undone_grid_points(int *undone,
				      char *phonon_done,
				      const int num_grid_points,
				      const int *grid_points);

void get_fc4_frequency_shifts(double *frequency_shifts,
			      const double *fc4_normal_real,
			      const double *frequencies,
			      const Iarray *grid_points1,
			      const Darray *temperatures,
			      const int *band_indicies,
			      const int num_band0,
			      const int num_band,
			      const double unit_conversion_factor)
{
  int i, j, k, l;
  double shift, num_phonon;

  for (i = 0; i < temperatures->dims[0]; i++) {
    for (j = 0; j < num_band0; j++) {
      shift = 0;
#pragma omp parallel for reduction(+:shift) private(k, l, num_phonon)
      for (k = 0; k < grid_points1->dims[0]; k++) {
	for (l = 0; l < num_band; l++) {
	  if (temperatures->data[i] > 0) {
	    num_phonon = 2 *
	      bose_einstein(frequencies[grid_points1->data[k] * num_band + l],
			    temperatures->data[i]) + 1;
	  } else {
	    num_phonon = 1;
	  }
	  shift += unit_conversion_factor * fc4_normal_real
	    [k * num_band0 * num_band + j * num_band + l] * num_phonon;
	}
      }
      frequency_shifts[i * num_band0 + j] = shift;
    }
  }
}

void
get_fc4_normal_for_frequency_shift(double *fc4_normal_real,
				   const double *frequencies,
				   const lapack_complex_double *eigenvectors,
				   const int grid_point0,
				   const Iarray *grid_points1,
				   const int *grid_address,
				   const int *mesh,
				   const double *fc4,
				   const Darray *shortest_vectors,
				   const Iarray *multiplicity,
				   const double *masses,
				   const int *p2s_map,
				   const int *s2p_map,
				   const Iarray *band_indicies,
				   const double cutoff_frequency)
{
  int i, num_atom, num_band, num_band0;

  num_atom = multiplicity->dims[1];
  num_band = num_atom * 3;
  num_band0 = band_indicies->dims[0];

#pragma omp parallel for private(i)
  for (i = 0; i < grid_points1->dims[0]; i++) {
    get_fc4_normal_for_frequency_shift_at_gp(fc4_normal_real +
					     i * num_band0 * num_band,
					     frequencies,
					     eigenvectors,
					     grid_point0,
					     grid_points1->data[i],
					     grid_address,
					     mesh,
					     fc4,
					     shortest_vectors,
					     multiplicity,
					     masses,
					     p2s_map,
					     s2p_map,
					     band_indicies,
					     cutoff_frequency);
  }
}

void set_phonons_for_frequency_shift(Darray *frequencies,
				     Carray *eigenvectors,
				     char *phonon_done,
				     const Iarray *grid_points,
				     const int *grid_address,
				     const int *mesh,
				     const Darray *fc2,
				     const Darray *svecs_fc2,
				     const Iarray *multi_fc2,
				     const double *masses_fc2,
				     const int *p2s_fc2,
				     const int *s2p_fc2,
				     const double unit_conversion_factor,
				     const double *born,
				     const double *dielectric,
				     const double *reciprocal_lattice,
				     const double *q_direction,
				     const double nac_factor,
				     const char uplo)
{
  int num_undone;
  int *undone;

  undone = (int*)malloc(sizeof(int) * frequencies->dims[0]);
  num_undone = collect_undone_grid_points(undone,
					  phonon_done,
					  grid_points->dims[0],
					  grid_points->data);

  get_undone_phonons(frequencies,
		     eigenvectors,
		     undone,
		     num_undone,
		     grid_address,
		     mesh,
		     fc2,
		     svecs_fc2,
		     multi_fc2,
		     masses_fc2,
		     p2s_fc2,
		     s2p_fc2,
		     unit_conversion_factor,
		     born,
		     dielectric,
		     reciprocal_lattice,
		     q_direction,
		     nac_factor,
		     uplo);

  free(undone);
}

void reciprocal_to_normal4(lapack_complex_double *fc4_normal,
			   const lapack_complex_double *fc4_reciprocal,
			   const double *freqs0,
			   const double *freqs1,
			   const lapack_complex_double *eigvecs0,
			   const lapack_complex_double *eigvecs1,
			   const double *masses,
			   const int *band_indices,
			   const int num_band0,
			   const int num_band,
			   const double cutoff_frequency)
{
  int i, j, bi, num_atom;
  lapack_complex_double fc4_sum_elem;

  num_atom = num_band / 3;

  for (i = 0; i < num_band0; i++) {
    bi = band_indices[i];
    if (freqs0[bi] > cutoff_frequency) {
      for (j = 0; j < num_band; j++) {
	if (freqs1[j] > cutoff_frequency) {
	  fc4_sum_elem = fc4_sum(bi,
				 j,
				 eigvecs0,
				 eigvecs1,
				 fc4_reciprocal,
				 masses,
				 num_atom);
	  fc4_normal[i * num_band + j] = lapack_make_complex_double
	    (lapack_complex_double_real(fc4_sum_elem) / freqs0[bi] / freqs1[j],
	     lapack_complex_double_imag(fc4_sum_elem) / freqs0[bi] / freqs1[j]);
	} else {
	  fc4_normal[i * num_band + j] = lapack_make_complex_double(0, 0);
	}
      }
    } else {
      for (j = 0; j < num_band; j++) {
	fc4_normal[i * num_band + j] = lapack_make_complex_double(0, 0);
      }
    }
  }
}

static void get_fc4_normal_for_frequency_shift_at_gp
(double *fc4_normal_real,
 const double *frequencies,
 const lapack_complex_double *eigenvectors,
 const int grid_point0,
 const int grid_point1,
 const int *grid_address,
 const int *mesh,
 const double *fc4,
 const Darray *shortest_vectors,
 const Iarray *multiplicity,
 const double *masses,
 const int *p2s_map,
 const int *s2p_map,
 const Iarray *band_indices,
 const double cutoff_frequency)
{
  int i, num_atom, num_band, num_band0;
  lapack_complex_double *fc4_reciprocal, *fc4_normal;
  double q[12];

  num_atom = multiplicity->dims[1];
  num_band = num_atom * 3;
  num_band0 = band_indices->dims[0];

  fc4_reciprocal = (lapack_complex_double*)
    malloc(sizeof(lapack_complex_double) *
	   num_atom * num_atom * num_atom * num_atom * 81);
  fc4_normal = (lapack_complex_double*)
    malloc(sizeof(lapack_complex_double) * num_band0 * num_band);

  for (i = 0; i < 3; i++) {
    q[i + 3] = (double)grid_address[grid_point0 * 3 + i] / mesh[i];
    q[i] = -q[i + 3];
    q[i + 6] = (double)grid_address[grid_point1 * 3 + i] / mesh[i];
    q[i + 9] = -q[i + 6];
  }

  real_to_reciprocal4(fc4_reciprocal,
		      q,
		      fc4,
		      shortest_vectors,
		      multiplicity,
		      p2s_map,
		      s2p_map);
  reciprocal_to_normal4(fc4_normal,
			fc4_reciprocal,
			frequencies + grid_point0 * num_band,
			frequencies + grid_point1 * num_band,
			eigenvectors + grid_point0 * num_band * num_band,
			eigenvectors + grid_point1 * num_band * num_band,
			masses,
			band_indices->data,
			num_band0,
			num_band,
			cutoff_frequency);

  for (i = 0; i < num_band0 * num_band; i++) {
    fc4_normal_real[i] = lapack_complex_double_real(fc4_normal[i]);
  }

  free(fc4_reciprocal);
  free(fc4_normal);
}

static lapack_complex_double fc4_sum(const int bi0,
				     const int bi1,
				     const lapack_complex_double *eigvecs0,
				     const lapack_complex_double *eigvecs1,
				     const lapack_complex_double *fc4_reciprocal,
				     const double *masses,
				     const int num_atom)
{
  int i, j, k, l, m, n, p, q;
  double sum_real, sum_imag, sum_real_cart, sum_imag_cart, mmm;
  lapack_complex_double eig_prod, eigvec0conj, eigvec1conj;

  sum_real = 0;
  sum_imag = 0;
  for (i = 0; i < num_atom; i++) {
    for (j = 0; j < num_atom; j++) {
      for (k = 0; k < num_atom; k++) {
	for (l = 0; l < num_atom; l++) {
	  sum_real_cart = 0;
	  sum_imag_cart = 0;
	  mmm = sqrt(masses[i] * masses[j] * masses[k] * masses[l]);
	  for (m = 0; m < 3; m++) {
	    eigvec0conj = lapack_make_complex_double
	      (lapack_complex_double_real
	       (eigvecs0[(i * 3 + m) * num_atom * 3 + bi0]),
	       -lapack_complex_double_imag
	       (eigvecs0[(i * 3 + m) * num_atom * 3 + bi0]));
	    for (n = 0; n < 3; n++) {
	    for (p = 0; p < 3; p++) {
	    for (q = 0; q < 3; q++) {
	      eigvec1conj = lapack_make_complex_double
		(lapack_complex_double_real
		 (eigvecs1[(l * 3 + q) * num_atom * 3 + bi1]),
		 -lapack_complex_double_imag
		 (eigvecs1[(l * 3 + q) * num_atom * 3 + bi1]));
	      eig_prod =
		phonoc_complex_prod(eigvec0conj,
		phonoc_complex_prod(eigvecs0[(j * 3 + n) * num_atom * 3 + bi0],
                phonoc_complex_prod(eigvecs1[(k * 3 + p) * num_atom * 3 + bi1],
                phonoc_complex_prod(eigvec1conj,
                fc4_reciprocal[i * num_atom * num_atom * num_atom * 81 +
			       j * num_atom * num_atom * 81 +
			       k * num_atom * 81 +
			       l * 81 +
			       m * 27 +
			       n * 9 +
			       p * 3 +
			       q]))));
	      sum_real_cart += lapack_complex_double_real(eig_prod);
	      sum_imag_cart += lapack_complex_double_imag(eig_prod);
	    }
	    }
	    }
	  }
	  sum_real += sum_real_cart / mmm;
	  sum_imag += sum_imag_cart / mmm;
	}
      }
    }
  }
  return lapack_make_complex_double(sum_real, sum_imag);
}

static int collect_undone_grid_points(int *undone,
				      char *phonon_done,
				      const int num_grid_points,
				      const int *grid_points)
{
  int i, gp, num_undone;

  num_undone = 0;

  for (i = 0; i < num_grid_points; i++) {
    gp = grid_points[i];
    if (phonon_done[gp] == 0) {
      undone[num_undone] = gp;
      num_undone++;
      phonon_done[gp] = 1;
    }
  }

  return num_undone;
}

max_stars_repo_name: atztogo/forcefit

max_issues_repo_name: atztogo/forcefit

max_forks_repo_name: atztogo/forcefit


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>

#ifdef SUBFIND

#include "fof.h"

#include "allvars.h"
#include "proto.h"
#include "domain.h"
#include "subfind.h"

static struct id_list
{
  MyIDType ID;
  int GrNr;
  int SubNr;
  float BindingEgy;

#ifdef SUBFIND_SAVE_PARTICLELISTS
  float Pos[3];
  float Vel[3];
  int Type;
#ifdef STELLARAGE
  float Mass;
  float StellarAge;
#endif
#endif
}
 *ID_list;

static int Nids;

void subfind(int num)
{
  double t0, t1, tstart, tend;
  int i, gr, nlocid, offset, limit, ncount;

#ifdef DENSITY_SPLIT_BY_TYPE
  struct unbind_data *d;
  int j, n, count[6], countall[6];
  double a3inv, dmax1, dmax2;
#endif

  if(ThisTask == 0)
    printf("\nWe now execute a parallel version of SUBFIND.\n");

  tstart = second();

#ifndef ONLY_PRODUCE_HSML_FILES
  if(!All.ComovingIntegrationOn)
    {
      if(ThisTask == 0)
	printf("works only for comoving integration.\n");
      endrun(0);
    }
#endif

#ifdef DENSITY_SPLIT_BY_TYPE
  a3inv = 1 / (All.Time * All.Time * All.Time);

  for(j = 0; j < 6; j++)
    count[j] = 0;

  /* let's count number of particles of selected species */
  for(i = 0; i < NumPart; i++)
    count[P[i].Type]++;

  MPI_Allreduce(count, countall, 6, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

  /* do first loop: basically just defining the hsml for different species */
  for(j = 0; j < 6; j++)
    {
      if((1 << j) & (DENSITY_SPLIT_BY_TYPE))
	{
#ifdef BLACK_HOLES
	  if(j == 5)
	    countall[j] = 0;	/* this will prevent that the black holes are treated separately */
#endif

	  force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);

	  if(countall[j] > All.DesNumNgb)
	    {
	      /* build index list of particles of selectes species */
	      d = (struct unbind_data *) mymalloc(count[j] * sizeof(struct unbind_data));
	      for(i = 0, n = 0; i < NumPart; i++)
		if(P[i].Type == j)
		  d[n++].index = i;

	      t0 = second();
	      if(ThisTask == 0)
		printf("Tree construction for species %d (%d).\n", j, countall[j]);

	      CPU_Step[CPU_FOF] += measure_time();

	      force_treebuild(count[j], d);

	      myfree(d);

	      t1 = second();
	      if(ThisTask == 0)
		printf("tree build for species %d took %g sec\n", j, timediff(t0, t1));
	    }
	  else
	    {
	      t0 = second();
	      if(ThisTask == 0)
		printf("Tree construction.\n");

	      CPU_Step[CPU_FOF] += measure_time();

	      force_treebuild(NumPart, NULL);

	      t1 = second();
	      if(ThisTask == 0)
		printf("tree build took %g sec\n", timediff(t0, t1));
	    }

	  /* let's determine the local densities */
	  t0 = second();
	  subfind_setup_smoothinglengths(j);
	  subfind_density(j);
	  t1 = second();
	  if(ThisTask == 0)
	    printf("density and smoothing length for species %d took %g sec\n", j, timediff(t0, t1));

	  force_treefree();

	  /* let's save density contribution of own species */
	  for(i = 0; i < NumPart; i++)
	    if(P[i].Type == j)
	      P[i].w.density_sum = P[i].u.DM_Density;

	}
    }

  /* do second loop: now calculate all density contributions */
  for(j = 0; j < 6; j++)
    {
      if((1 << j) & (DENSITY_SPLIT_BY_TYPE))
	{
	  force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);

	  /* build index list of particles of selectes species */
	  d = (struct unbind_data *) mymalloc(count[j] * sizeof(struct unbind_data));
	  for(i = 0, n = 0; i < NumPart; i++)
	    if(P[i].Type == j)
	      d[n++].index = i;

	  t0 = second();
	  if(ThisTask == 0)
	    printf("Tree construction for species %d (%d).\n", j, countall[j]);

	  CPU_Step[CPU_FOF] += measure_time();

	  force_treebuild(count[j], d);

	  myfree(d);

	  t1 = second();
	  if(ThisTask == 0)
	    printf("tree build for species %d took %g sec\n", j, timediff(t0, t1));

	  /* let's determine the local densities */
	  t0 = second();
	  for(i = 0; i < 6; i++)
	    if((1 << i) & (DENSITY_SPLIT_BY_TYPE))
	      if(j != i)
		{
		  if(countall[i] > All.DesNumNgb)
		    {
		      if(ThisTask == 0)
			printf("calculating density contribution of species %d to species %d\n", j, i);
		      subfind_density(-(i + 1));
		    }
		}
	  t1 = second();
	  if(ThisTask == 0)
	    printf("density() of species %d took %g sec\n", j, timediff(t0, t1));

	  force_treefree();

	  /* let's sum up density contribution */
	  for(i = 0; i < NumPart; i++)
	    if((1 << P[i].Type) & (DENSITY_SPLIT_BY_TYPE))
	      if(j != P[i].Type)
		if(countall[P[i].Type] > All.DesNumNgb)
		  P[i].w.density_sum += P[i].u.DM_Density;
	}
    }

  for(i = 0; i < NumPart; i++)
    {
      P[i].u.DM_Density = P[i].w.density_sum;

      if(P[i].Type == 0)
	P[i].w.int_energy = DMAX(All.MinEgySpec,
				 SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].d.Density * a3inv,
								      GAMMA_MINUS1));
      else
	P[i].w.int_energy = 0;

    }
#else
  force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);

  t0 = second();
  if(ThisTask == 0)
    printf("Tree construction.\n");

  CPU_Step[CPU_FOF] += measure_time();

  force_treebuild(NumPart, NULL);

  t1 = second();
  if(ThisTask == 0)
    printf("tree build took %g sec\n", timediff(t0, t1));

  /* let's determine the local dark matter densities */
  t0 = second();
  subfind_setup_smoothinglengths();
  subfind_density();
  t1 = second();
  if(ThisTask == 0)
    printf("dark matter density() took %g sec\n", timediff(t0, t1));

  force_treefree();
#endif /* DENSITY_SPLIT_BY_TYPE */

  if(DumpFlag)
    {
      /* let's save the densities to a file (for making images) */
      t0 = second();
      subfind_save_densities(num);
      t1 = second();
      if(ThisTask == 0)
	printf("saving densities took %g sec\n", timediff(t0, t1));
    }

#ifdef ONLY_PRODUCE_HSML_FILES
  return;
#endif

  /* count how many groups we have that should be done collectively */
  limit = 0.6 * All.TotNumPart / NTask;

  for(i = 0, ncount = 0; i < Ngroups; i++)
    if(Group[i].Len >= limit)
      ncount++;
  MPI_Allreduce(&ncount, &Ncollective, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      printf("\nNumber of FOF halos treated with collective SubFind code = %d\n", Ncollective);
      printf("(the adopted size-limit for the collective algorithm was %d particles.)\n", limit);
      printf("the other %d FOF halos are treated in parallel with serial code\n\n", TotNgroups - Ncollective);
    }

  /*  to decide on which task a group should be:
   *  if   GrNr <= Ncollective:  collective groupfinding.
   *  the task where the group info is put is TaskNr = (GrNr - 1) % NTask
   */

  /* now we distribute the particles such that small groups are assigned in
   *  total to certain CPUs, and big groups are left where they are
   */

  t0 = second();

  for(i = 0; i < NumPart; i++)
    {
      if(P[i].GrNr > Ncollective && P[i].GrNr <= TotNgroups)	/* particle is in small group */
	P[i].targettask = (P[i].GrNr - 1) % NTask;
      else
	P[i].targettask = ThisTask;
    }

  subfind_exchange();		/* distributes gas particles as well if needed */

  t1 = second();
  if(ThisTask == 0)
    printf("subfind_exchange()() took %g sec\n", timediff(t0, t1));

  subfind_distribute_groups();

  qsort(Group, Ngroups, sizeof(struct group_properties), fof_compare_Group_GrNr);

  for(i = 0; i < NumPart; i++)
    if(P[i].GrNr > Ncollective && P[i].GrNr <= TotNgroups)
      if(((P[i].GrNr - 1) % NTask) != ThisTask)
	{
	  printf("i=%d %d task=%d\n", i, P[i].GrNr, ThisTask);
	  endrun(87);
	}

  /* lets estimate the maximum number of substructures we need to store on the local CPU */
  for(i = 0, nlocid = 0; i < Ngroups; i++)
    nlocid += Group[i].Len;

  MaxNsubgroups = nlocid / All.DesLinkNgb;	/* this is a quite conservative upper limit */
  Nsubgroups = 0;
  SubGroup = (struct subgroup_properties *) mymalloc(MaxNsubgroups * sizeof(struct subgroup_properties));

  for(i = 0; i < NumPart; i++)
    P[i].SubNr = (1 << 30);	/* default */

  /* we begin by applying the collective version of subfind to distributed groups */
  t0 = second();
  for(GrNr = 1; GrNr <= Ncollective; GrNr++)
    subfind_process_group_collectively(num);
  t1 = second();
  if(ThisTask == 0)
    printf("processing of collective halos took %g sec\n", timediff(t0, t1));

#ifdef SUBFIND_COLLECTIVE_STAGE1
  if(ThisTask == 0)
    printf("stage 1 ended\n");
  endrun(0);
#endif

  for(i = 0; i < NumPart; i++)
    {
      P[i].origindex = i;
      P[i].origintask = ThisTask;
    }

  t0 = second();
  qsort(P, NumPart, sizeof(struct particle_data), subfind_compare_P_GrNr_DM_Density);
  t1 = second();
  if(ThisTask == 0)
    printf("sort of local particles()() took %g sec\n", timediff(t0, t1));

  /* now we have the particles of groups consecutively, but SPH particles are
     not aligned. They can however be accessed via SphP[P[i].originindex] */

  /* let's count how many local particles we have in small groups */
  for(i = 0, nlocid = 0; i < NumPart; i++)
    if(P[i].GrNr > Ncollective && P[i].GrNr <= Ngroups)	/* particle is in small group */
      nlocid++;

  if(ThisTask == 0)
    printf("contructing tree for serial subfind of local groups\n");

  subfind_loctree_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);

  if(ThisTask == 0)
    printf("Start to do local groups with serial subfind algorithm\n");

  t0 = second();

  /* we now apply a serial version of subfind to the local groups */
  for(gr = 0, offset = 0; gr < Ngroups; gr++)
    {
      if(Group[gr].GrNr > Ncollective)
	{
	  if(((Group[gr].GrNr - 1) % NTask) == ThisTask)
	    offset = subfind_process_group_serial(gr, offset);
	}
    }

  t1 = second();
  if(ThisTask == 0)
    printf("processing of local groups took took %g sec\n", timediff(t0, t1));

  subfind_loctree_treefree();

  /* bringing back particles in original positions, such that gas particles are aligned */
  t0 = second();
  qsort(P, NumPart, sizeof(struct particle_data), subfind_compare_P_origindex);
  t1 = second();
  if(ThisTask == 0)
    printf("unsorting of local particles()() took %g sec\n", timediff(t0, t1));

  GrNr = -1;			/* to ensure that domain decomposition acts normally again */

  /* now determine the remaining spherical overdensity values for the non-local groups */

  domain_free_trick();

  CPU_Step[CPU_FOF] += measure_time();

#ifdef DENSITY_SPLIT_BY_TYPE
  printf("Task %d: testing particles ...\n", ThisTask);
  for(i = 0; i < NumPart; i++)
    {
      if(P[i].origintask != ThisTask)
	printf("Task %d: Holding particle of task %d !\n", ThisTask, P[i].origintask);
      if(P[i].origindex != i)
	printf("Task %d: Particles is in wrong position (is=%d, was=%d) !\n", ThisTask, i, P[i].origindex);
    }
#endif

  All.DoDynamicUpdate = 0;
  domain_Decomposition();

  force_treebuild(NumPart, NULL);

  /* compute spherical overdensities for FOF groups */
  t0 = second();

  subfind_overdensity();

  t1 = second();
  if(ThisTask == 0)
    printf("determining spherical overdensity masses took %g sec\n", timediff(t0, t1));

  /* determine which halos are contaminated by boundary particles */
  t0 = second();

  subfind_contamination();

  t1 = second();
  if(ThisTask == 0)
    printf("determining contamination of halos took %g sec\n", timediff(t0, t1));

  force_treefree();
  domain_free();

  domain_allocate_trick();

  /* now assemble final output */
  subfind_save_final(num);

  tend = second();

  if(ThisTask == 0)
    printf("\nFinished with SUBFIND.  (total time=%g sec)\n\n", timediff(tstart, tend));

  myfree(SubGroup);

  CPU_Step[CPU_FOF] += measure_time();
}

void subfind_save_final(int num)
{
  int i, j, totsubs, masterTask, groupTask, nprocgroup;
  char buf[1000];
  double t0, t1;

  /* prepare list of ids with assigned group numbers */

  parallel_sort(Group, Ngroups, sizeof(struct group_properties), fof_compare_Group_GrNr);
  parallel_sort(SubGroup, Nsubgroups, sizeof(struct subgroup_properties),
		subfind_compare_SubGroup_GrNr_SubNr);

  ID_list = mymalloc(sizeof(struct id_list) * NumPart);

  for(i = 0, Nids = 0; i < NumPart; i++)
    {
      if(P[i].GrNr <= TotNgroups)
	{
	  ID_list[Nids].GrNr = P[i].GrNr;
	  ID_list[Nids].SubNr = P[i].SubNr;
	  ID_list[Nids].BindingEgy = P[i].v.DM_BindingEnergy;
	  ID_list[Nids].ID = P[i].ID;
#ifdef SUBFIND_SAVE_PARTICLELISTS
	  for(j = 0; j < 3; j++)
	    {
	      ID_list[Nids].Pos[j] = P[i].Pos[j];
	      ID_list[Nids].Vel[j] = P[i].Vel[j];
	    }
	  ID_list[Nids].Type = P[i].Type;
#ifdef STELLARAGE
          ID_list[Nids].Mass = P[i].Mass;
          if(P[i].Type == 4)
            ID_list[Nids].StellarAge = P[i].StellarAge;
          else
            ID_list[Nids].StellarAge = 0;
#endif
#endif
	  Nids++;
	}
    }

  parallel_sort(ID_list, Nids, sizeof(struct id_list), subfind_compare_ID_list);

  MPI_Allreduce(&Nsubgroups, &TotNsubgroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

  /* fill in the FirstSub-values */
  for(i = 0, totsubs = 0; i < Ngroups; i++)
    {
      if(i > 0)
	Group[i].FirstSub = Group[i - 1].FirstSub + Group[i - 1].Nsubs;
      else
	Group[i].FirstSub = 0;
      totsubs += Group[i].Nsubs;
    }

  MPI_Allgather(&totsubs, 1, MPI_INT, Send_count, 1, MPI_INT, MPI_COMM_WORLD);
  for(j = 1, Send_offset[0] = 0; j < NTask; j++)
    Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];

  for(i = 0; i < Ngroups; i++)
    Group[i].FirstSub += Send_offset[ThisTask];

  MPI_Allgather(&Nids, 1, MPI_INT, Send_count, 1, MPI_INT, MPI_COMM_WORLD);
  for(j = 1, Send_offset[0] = 0; j < NTask; j++)
    Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];

  if(ThisTask == 0)
    {
      sprintf(buf, "%s/groups_%03d", All.OutputDir, num);
      mkdir(buf, 02755);
    }
  MPI_Barrier(MPI_COMM_WORLD);

  if(NTask < All.NumFilesWrittenInParallel)
    {
      printf
	("Fatal error.\nNumber of processors must be a smaller or equal than `NumFilesWrittenInParallel'.\n");
      endrun(241931);
    }

  t0 = second();

  nprocgroup = NTask / All.NumFilesWrittenInParallel;
  if((NTask % All.NumFilesWrittenInParallel))
    nprocgroup++;
  masterTask = (ThisTask / nprocgroup) * nprocgroup;
  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))	/* ok, it's this processor's turn */
	subfind_save_local_catalogue(num);
      MPI_Barrier(MPI_COMM_WORLD);	/* wait inside the group */
    }

  t1 = second();

  if(ThisTask == 0)
    {
      printf("Subgroup catalogues saved. took = %g sec\n", timediff(t0, t1));
      fflush(stdout);
    }

  myfree(ID_list);
}

void subfind_save_local_catalogue(int num)
{
  FILE *fd;
  char buf[500], fname[500];
  float *mass, *pos, *vel, *spin;

#ifdef SAVE_MASS_TAB
  float *masstab;
#endif
  int i, j, *len;
  MyIDType *ids;

  sprintf(fname, "%s/groups_%03d/%s_%03d.%d", All.OutputDir, num, "subhalo_tab", num, ThisTask);
  strcpy(buf, fname);
  if(!(fd = fopen(buf, "w")))
    {
      printf("can't open file `%s`\n", buf);
      endrun(1183);
    }

  my_fwrite(&Ngroups, sizeof(int), 1, fd);
  my_fwrite(&TotNgroups, sizeof(int), 1, fd);
  my_fwrite(&Nids, sizeof(int), 1, fd);
  my_fwrite(&TotNids, sizeof(long long), 1, fd);
  my_fwrite(&NTask, sizeof(int), 1, fd);
  my_fwrite(&Nsubgroups, sizeof(int), 1, fd);
  my_fwrite(&TotNsubgroups, sizeof(int), 1, fd);

  /* group len */
  len = mymalloc(Ngroups * sizeof(int));
  for(i = 0; i < Ngroups; i++)
    len[i] = Group[i].Len;
  my_fwrite(len, Ngroups, sizeof(int), fd);
  myfree(len);

  /* offset into id-list */
  len = mymalloc(Ngroups * sizeof(int));
  for(i = 0; i < Ngroups; i++)
    len[i] = Group[i].Offset;
  my_fwrite(len, Ngroups, sizeof(int), fd);
  myfree(len);

  /* mass */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].Mass;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* location (potential minimum) */
  pos = mymalloc(Ngroups * 3 * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    for(j = 0; j < 3; j++)
      pos[i * 3 + j] = Group[i].Pos[j];
  my_fwrite(pos, Ngroups, 3 * sizeof(float), fd);
  myfree(pos);

  /* M_Mean200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].M_Mean200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* R_Mean200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].R_Mean200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* M_Crit200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].M_Crit200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* R_Crit200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].R_Crit200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* M_TopHat200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].M_TopHat200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* R_TopHat200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].R_TopHat200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

#ifdef SO_VEL_DISPERSIONS
  /* VelDisp_Mean200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].VelDisp_Mean200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* VelDisp_Crit200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].VelDisp_Crit200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* VelDisp_TopHat200 */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].VelDisp_TopHat200;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);
#endif

  /* contamination particle count */
  len = mymalloc(Ngroups * sizeof(int));
  for(i = 0; i < Ngroups; i++)
    len[i] = Group[i].ContaminationLen;
  my_fwrite(len, Ngroups, sizeof(int), fd);
  myfree(len);

  /* contamination mass */
  mass = mymalloc(Ngroups * sizeof(float));
  for(i = 0; i < Ngroups; i++)
    mass[i] = Group[i].ContaminationMass;
  my_fwrite(mass, Ngroups, sizeof(float), fd);
  myfree(mass);

  /* number of substructures in FOF group  */
  len = mymalloc(Ngroups * sizeof(int));
  for(i = 0; i < Ngroups; i++)
    len[i] = Group[i].Nsubs;
  my_fwrite(len, Ngroups, sizeof(int), fd);
  myfree(len);

  /* first substructure in FOF group  */
  len = mymalloc(Ngroups * sizeof(int));
  for(i = 0; i < Ngroups; i++)
    len[i] = Group[i].FirstSub;
  my_fwrite(len, Ngroups, sizeof(int), fd);
  myfree(len);

  /* ------------------------------ */

  /* Len of substructure  */
  len = mymalloc(Nsubgroups * sizeof(int));
  for(i = 0; i < Nsubgroups; i++)
    len[i] = SubGroup[i].Len;
  my_fwrite(len, Nsubgroups, sizeof(int), fd);
  myfree(len);

  /* offset of substructure  */
  len = mymalloc(Nsubgroups * sizeof(int));
  for(i = 0; i < Nsubgroups; i++)
    len[i] = SubGroup[i].Offset;
  my_fwrite(len, Nsubgroups, sizeof(int), fd);
  myfree(len);

  /* parent of substructure  */
  len = mymalloc(Nsubgroups * sizeof(int));
  for(i = 0; i < Nsubgroups; i++)
    len[i] = SubGroup[i].SubParent;
  my_fwrite(len, Nsubgroups, sizeof(int), fd);
  myfree(len);

  /* Mass of substructure  */
  mass = mymalloc(Nsubgroups * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    mass[i] = SubGroup[i].Mass;
  my_fwrite(mass, Nsubgroups, sizeof(float), fd);
  myfree(mass);

  /* Pos of substructure  */
  pos = mymalloc(Nsubgroups * 3 * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    for(j = 0; j < 3; j++)
      pos[i * 3 + j] = SubGroup[i].Pos[j];
  my_fwrite(pos, Nsubgroups, 3 * sizeof(float), fd);
  myfree(pos);

  /* Vel of substructure  */
  vel = mymalloc(Nsubgroups * 3 * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    for(j = 0; j < 3; j++)
      vel[i * 3 + j] = SubGroup[i].Vel[j];
  my_fwrite(vel, Nsubgroups, 3 * sizeof(float), fd);
  myfree(vel);

  /* Center of mass of substructure  */
  pos = mymalloc(Nsubgroups * 3 * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    for(j = 0; j < 3; j++)
      pos[i * 3 + j] = SubGroup[i].CM[j];
  my_fwrite(pos, Nsubgroups, 3 * sizeof(float), fd);
  myfree(pos);

  /* Spin of substructure  */
  spin = mymalloc(Nsubgroups * 3 * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    for(j = 0; j < 3; j++)
      spin[i * 3 + j] = SubGroup[i].Spin[j];
  my_fwrite(spin, Nsubgroups, 3 * sizeof(float), fd);
  myfree(spin);

  /* velocity dispesion  */
  mass = mymalloc(Nsubgroups * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    mass[i] = SubGroup[i].SubVelDisp;
  my_fwrite(mass, Nsubgroups, sizeof(float), fd);
  myfree(mass);

  /* maximum circular velocity  */
  mass = mymalloc(Nsubgroups * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    mass[i] = SubGroup[i].SubVmax;
  my_fwrite(mass, Nsubgroups, sizeof(float), fd);
  myfree(mass);

  /* radius of maximum circular velocity  */
  mass = mymalloc(Nsubgroups * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    mass[i] = SubGroup[i].SubVmaxRad;
  my_fwrite(mass, Nsubgroups, sizeof(float), fd);
  myfree(mass);

  /* radius of half the mass  */
  mass = mymalloc(Nsubgroups * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    mass[i] = SubGroup[i].SubHalfMass;
  my_fwrite(mass, Nsubgroups, sizeof(float), fd);
  myfree(mass);

  /* ID of most bound particle  */
  ids = mymalloc(Nsubgroups * sizeof(MyIDType));
  for(i = 0; i < Nsubgroups; i++)
    ids[i] = SubGroup[i].SubMostBoundID;
  my_fwrite(ids, Nsubgroups, sizeof(MyIDType), fd);
  myfree(ids);

  /* GrNr of substructure  */
  len = mymalloc(Nsubgroups * sizeof(int));
  for(i = 0; i < Nsubgroups; i++)
    len[i] = SubGroup[i].GrNr;
  my_fwrite(len, Nsubgroups, sizeof(int), fd);
  myfree(len);

  /* Masstab of substructure  */
#ifdef SAVE_MASS_TAB
  masstab = mymalloc(Nsubgroups * 6 * sizeof(float));
  for(i = 0; i < Nsubgroups; i++)
    for(j = 0; j < 6; j++)
      masstab[i * 6 + j] = SubGroup[i].MassTab[j];
  my_fwrite(masstab, Nsubgroups, 6 * sizeof(float), fd);
  myfree(masstab);
#endif

  fclose(fd);

#ifdef SUBFIND_SAVE_PARTICLELISTS
  sprintf(buf, "%s/groups_%03d/%s_%03d.%d", All.OutputDir, num, "subhalo_posvel", num, ThisTask);
  if(!(fd = fopen(buf, "w")))
    {
      printf("can't open file `%s`\n", buf);
      endrun(1184);
    }

  my_fwrite(&Ngroups, sizeof(int), 1, fd);
  my_fwrite(&TotNgroups, sizeof(int), 1, fd);
  my_fwrite(&Nids, sizeof(int), 1, fd);
  my_fwrite(&TotNids, sizeof(long long), 1, fd);
  my_fwrite(&NTask, sizeof(int), 1, fd);
  my_fwrite(&Send_offset[ThisTask], sizeof(int), 1, fd);
  my_fwrite(&All.Time, sizeof(double), 1, fd);

  float *posdata;
  double a3inv;
  char *types;

  if(All.ComovingIntegrationOn)
    a3inv = 1.0 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;

  posdata = (float *) mymalloc(3 * Nids * sizeof(float));

  for(i = 0; i < Nids; i++)
    for(j = 0; j < 3; j++)
      posdata[i * 3 + j] = ID_list[i].Pos[j];

  my_fwrite(posdata, 3 * sizeof(float), Nids, fd);

  for(i = 0; i < Nids; i++)
    for(j = 0; j < 3; j++)
      posdata[i * 3 + j] = ID_list[i].Vel[j] * sqrt(a3inv);

  my_fwrite(posdata, 3 * sizeof(float), Nids, fd);

  types = (char *) posdata;

  for(i = 0; i < Nids; i++)
    types[i] = ID_list[i].Type;

  my_fwrite(types, sizeof(char), Nids, fd);

#ifdef STELLARAGE
  for(i = 0; i < Nids; i++)
    posdata[i] = ID_list[i].Mass;

  my_fwrite(posdata, sizeof(float), Nids, fd);

  for(i = 0; i < Nids; i++)
    posdata[i] = ID_list[i].StellarAge;

  my_fwrite(posdata, sizeof(float), Nids, fd);
#endif

  myfree(posdata);

  fclose(fd);
#endif

  ids = (MyIDType *) ID_list;

  for(i = 0; i < Nids; i++)
    ids[i] = ID_list[i].ID;

  sprintf(buf, "%s/groups_%03d/%s_%03d.%d", All.OutputDir, num, "subhalo_ids", num, ThisTask);
  if(!(fd = fopen(buf, "w")))
    {
      printf("can't open file `%s`\n", buf);
      endrun(1184);
    }

  my_fwrite(&Ngroups, sizeof(int), 1, fd);
  my_fwrite(&TotNgroups, sizeof(int), 1, fd);
  my_fwrite(&Nids, sizeof(int), 1, fd);
  my_fwrite(&TotNids, sizeof(long long), 1, fd);
  my_fwrite(&NTask, sizeof(int), 1, fd);
  my_fwrite(&Send_offset[ThisTask], sizeof(int), 1, fd);
  my_fwrite(ids, sizeof(MyIDType), Nids, fd);

  fclose(fd);
}

int subfind_compare_ID_list(const void *a, const void *b)
{
  if(((struct id_list *) a)->GrNr < ((struct id_list *) b)->GrNr)
    return -1;

  if(((struct id_list *) a)->GrNr > ((struct id_list *) b)->GrNr)
    return +1;

  if(((struct id_list *) a)->SubNr < ((struct id_list *) b)->SubNr)
    return -1;

  if(((struct id_list *) a)->SubNr > ((struct id_list *) b)->SubNr)
    return +1;

  if(((struct id_list *) a)->BindingEgy < ((struct id_list *) b)->BindingEgy)
    return -1;

  if(((struct id_list *) a)->BindingEgy > ((struct id_list *) b)->BindingEgy)
    return +1;

  return 0;
}

int subfind_compare_SubGroup_GrNr_SubNr(const void *a, const void *b)
{
  if(((struct subgroup_properties *) a)->GrNr < ((struct subgroup_properties *) b)->GrNr)
    return -1;

  if(((struct subgroup_properties *) a)->GrNr > ((struct subgroup_properties *) b)->GrNr)
    return +1;

  if(((struct subgroup_properties *) a)->SubNr < ((struct subgroup_properties *) b)->SubNr)
    return -1;

  if(((struct subgroup_properties *) a)->SubNr > ((struct subgroup_properties *) b)->SubNr)
    return +1;

  return 0;
}

int subfind_compare_P_GrNr_DM_Density(const void *a, const void *b)
{
  if(((struct particle_data *) a)->GrNr < (((struct particle_data *) b)->GrNr))
    return -1;

  if(((struct particle_data *) a)->GrNr > (((struct particle_data *) b)->GrNr))
    return +1;

  if(((struct particle_data *) a)->u.DM_Density > (((struct particle_data *) b)->u.DM_Density))
    return -1;

  if(((struct particle_data *) a)->u.DM_Density < (((struct particle_data *) b)->u.DM_Density))
    return +1;

  return 0;
}

int subfind_compare_P_origindex(const void *a, const void *b)
{
  if(((struct particle_data *) a)->origindex < (((struct particle_data *) b)->origindex))
    return -1;

  if(((struct particle_data *) a)->origindex > (((struct particle_data *) b)->origindex))
    return +1;

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#ifndef PVIZUTIL_H
#define PVIZUTIL_H

#include <math.h>
#include <string>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include <H5Cpp.h>

using namespace std;
using namespace H5;

class PvizUtil
{
public:
	static void split(const int n, const int nsplit, int* counts, int* offsets)
	{
		if (nsplit == 0)
			throw "Cannot split by zero";

		// i = 0
		offsets[0] = 0;
		counts[0] = ceil(n/(double)nsplit);

		// i > 0
		for (int i = 1; i < nsplit; i++)
		{
			offsets[i] = offsets[i-1] + counts[i-1];
			counts[i] = ceil((n - offsets[i])/(double)(nsplit - i));
		}
	};

	static bool H5CheckDataset(const std::string &filename,
							   const std::string &datasetname)
	{
		bool rtn = true;
		try {
			H5::H5File file(filename, H5F_ACC_RDONLY );
			H5::DataSet dataset = file.openDataSet(datasetname);
		}
		catch (H5::FileIException e)
		{
			rtn = false;
		}

		return rtn;
	};

	static void H5GetDatasetDim(const std::string &filename,
								const std::string &datasetname,
								const int &rank,
								hsize_t *dim)
	{
		H5::H5File file( filename, H5F_ACC_RDONLY );
		H5::DataSet dataset = file.openDataSet(datasetname);

		/*
		 * Get dataspace of the dataset.
		 */
		H5::DataSpace dataspace = dataset.getSpace();

		/*
		 * Get the number of dimensions in the dataspace.
		 */
		if (rank != dataspace.getSimpleExtentNdims())
		{
			throw std::invalid_argument("Rank mismatch!!");
		}

		dataspace.getSimpleExtentDims(dim, NULL);
	};

	static void H5ReadByOffset(const std::string &filename,
							   const std::string &datasetname,
							   const H5::DataType &mem_type,
							   const size_t &rowOffset,
							   const size_t &size,
							   void *m)
	{
		H5::H5File file(filename, H5F_ACC_RDONLY);
		H5::DataSet dataset = file.openDataSet(datasetname);
		H5::DataSpace dataspace = dataset.getSpace();

		checkDataset(dataset, mem_type, rowOffset, size, m);

		hsize_t      offset[1];   // hyperslab offset in the file
		hsize_t      count[1];    // size of the hyperslab in the file
		offset[0] = rowOffset;
		count[0]  = size;

		dataspace.selectHyperslab(H5S_SELECT_SET, count, offset);

		H5::DataSpace memspace(1, count);

		dataset.read(m, mem_type, memspace, dataspace);
	};

	static void H5ReadByRow(const std::string &filename,
							const std::string &datasetname,
							const H5::DataType &mem_type,
							const size_t &rowOffset,
							const size_t &size1,
							const size_t &size2,
							void *m)
	{
		H5::H5File file(filename, H5F_ACC_RDONLY);
		H5::DataSet dataset = file.openDataSet(datasetname);
		H5::DataSpace dataspace = dataset.getSpace();

		checkDataset(dataset, mem_type, rowOffset, size1, size2, m);

		hsize_t      offset[2];   // hyperslab offset in the file
		hsize_t      count[2];    // size of the hyperslab in the file
		offset[0] = rowOffset;
		offset[1] = 0;
		count[0]  = size1;
		count[1]  = size2;

		dataspace.selectHyperslab(H5S_SELECT_SET, count, offset);

		H5::DataSpace memspace(2, count);

		dataset.read(m, mem_type, memspace, dataspace);
	};

	static void H5WriteByRow(const std::string &filename,
							 const std::string &datasetname,
							 const H5::DataType &mem_type,
							 const size_t &rowOffset,
							 const size_t &size1,
							 const size_t &size2,
							 const gsl_matrix *m)
	{
		H5::H5File file(filename, H5F_ACC_RDWR);
		H5::DataSet dataset = file.openDataSet(datasetname);
		H5::DataSpace dataspace = dataset.getSpace();

		checkDataset(dataset, mem_type, rowOffset, size1, size2, m);

		hsize_t      offset[2];   // hyperslab offset in the file
		hsize_t      count[2];    // size of the hyperslab in the file
		offset[0] = rowOffset;
		offset[1] = 0;
		count[0]  = size1;
		count[1]  = size2;

		dataspace.selectHyperslab(H5S_SELECT_SET, count, offset);

		H5::DataSpace memspace(2, count);

		dataset.write(m, mem_type, memspace, dataspace);
	};

protected:
	static void checkDataset(const H5::DataSet& dataset,
							 const H5::DataType &mem_type,
							 const size_t rowOffset,
							 const size_t &size1,
							 const size_t &size2,
							 const void *m)
	{
		// Check datatype
		/*
		H5T_class_t class_t = dataset.getTypeClass();

		if ((class_t != H5T_FLOAT) && (class_t != H5T_INTEGER))
		{
			throw std::runtime_error("Type mismatch!!");
		}
		*/

		H5::DataType dtype = dataset.getDataType();
		if (!(dtype == mem_type))
		{
			throw std::runtime_error("Data type mismatch!!");
		}

		H5::DataSpace dataspace = dataset.getSpace();

		// Check rank
		if (dataspace.getSimpleExtentNdims() != 2)
		{
			throw std::runtime_error("Rank mismatch!!");
		}

		// Check dimension
		hsize_t dims_out[2];
		dataspace.getSimpleExtentDims(dims_out, NULL);
		if (dims_out[0] < size1)
		{
			throw std::runtime_error("Dimension mismatch!!");
		}

		if (dims_out[1] != size2)
		{
			throw std::runtime_error("Dimension mismatch!!");
		}

		// Check parameter
		if (dims_out[0] <= rowOffset)
		{
			throw std::runtime_error("Offset is out of range!!");
		}

		if (dims_out[0] < rowOffset + size1)
		{
			throw std::runtime_error("Count is out of range!!");
		}
	};

	static void checkDataset(const H5::DataSet& dataset,
							 const H5::DataType &mem_type,
							 const size_t rowOffset,
							 const size_t &size,
							 const void *m)
	{
		// Check datatype
		/*
		 H5T_class_t class_t = dataset.getTypeClass();

		 if ((class_t != H5T_FLOAT) && (class_t != H5T_INTEGER))
		 {
		 throw std::runtime_error("Type mismatch!!");
		 }
		 */

		H5::DataType dtype = dataset.getDataType();
		if (!(dtype == mem_type))
		{
			throw std::runtime_error("Data type mismatch!!");
		}

		H5::DataSpace dataspace = dataset.getSpace();

		// Check rank
		if (dataspace.getSimpleExtentNdims() != 1)
		{
			throw std::runtime_error("Rank mismatch!!");
		}

		// Check dimension
		hsize_t dims_out[1];
		dataspace.getSimpleExtentDims(dims_out, NULL);
		if (dims_out[0] < size)
		{
			throw std::runtime_error("Dimension mismatch!!");
		}

		// Check parameter
		if (dims_out[0] <= rowOffset)
		{
			throw std::runtime_error("Offset is out of range!!");
		}

		if (dims_out[0] < rowOffset + size)
		{
			throw std::runtime_error("Count is out of range!!");
		}
	};
};

#endif //PVIZUTIL_H

ext: h


text:
/**
 * @file      Sampler.h
 *
 * @brief     This private head contains the definition of the Sampler type
 *
 * @author    Stefan Reinhold
 * @copyright Copyright (C) 2018 Stefan Reinhold  -- All Rights Reserved.
 *            You may use, distribute and modify this code under the terms of
 *            the AFL 3.0 license; see LICENSE for full license details.
 */

#pragma once

#include "MeasurementModel.h"
#include "VoxelVolume.h"

#include <Eigen/Core>
#include <gsl/gsl>

namespace CortidQCT {

#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wpadded"
class VolumeSampler {

public:
  VoxelVolume::ValueType outside;

  inline explicit VolumeSampler(
      VoxelVolume const &vol,
      VoxelVolume::ValueType outside_ = VoxelVolume::ValueType{0}) noexcept
      : outside{outside_}, volume_{vol} {}

  template <class Derived, class DerivedOut>
  inline void operator()(Eigen::MatrixBase<Derived> const &positions,
                         Eigen::MatrixBase<DerivedOut> &values,
                         VoxelVolume::ValueType slope,
                         VoxelVolume::ValueType intercept) const {
    using Eigen::Vector3f;

    Expects(values.rows() == positions.rows());
    Expects(positions.cols() == 3);

    Vector3f const scale{1.f / volume_.voxelSize().width,
                         1.f / volume_.voxelSize().height,
                         1.f / volume_.voxelSize().depth};

    volume_.withUnsafeDataPointer([this, &positions, &values, scale, slope,
                                   intercept](auto const *ptr) {

#pragma omp parallel for
      for (Eigen::Index i = 0; i < positions.rows(); ++i) {

        values(i) = this->interpolate(
                        (positions.row(i).array() * scale.array().transpose())
                            .matrix()
                            .transpose(),
                        gsl::make_not_null(ptr)) *
                        slope +
                    intercept;
      }
    });
  }

  template <class Derived>
  inline Eigen::Matrix<VoxelVolume::ValueType, Eigen::Dynamic, 1>
  operator()(Eigen::MatrixBase<Derived> const &positions,
             VoxelVolume::ValueType slope,
             VoxelVolume::ValueType intercept) const {
    Eigen::Matrix<VoxelVolume::ValueType, Eigen::Dynamic, 1> values(
        positions.rows());

    operator()(positions, values, slope, intercept);

    return values;
  }

private:
  template <class Derived>
  inline VoxelVolume::ValueType
  at(Eigen::MatrixBase<Derived> const &pos,
     gsl::not_null<VoxelVolume::ValueType const *> ptr) const {
    auto const &size = volume_.size();

    auto const isInside = (pos.array() >= 0).all() &&
                          pos(0) < gsl::narrow_cast<int>(size.width) &&
                          pos(1) < gsl::narrow_cast<int>(size.height) &&
                          pos(2) < gsl::narrow_cast<int>(size.depth);

    auto const index =
        gsl::narrow_cast<std::size_t>(pos(2)) * size.width * size.height +
        gsl::narrow_cast<std::size_t>(pos(1)) * size.width +
        gsl::narrow_cast<std::size_t>(pos(0));

    return isInside ? ptr.get()[index] : outside;
  }

  template <class Derived>
  inline auto
  interpolate(Eigen::MatrixBase<Derived> const &pos,
              gsl::not_null<VoxelVolume::ValueType const *> ptr) const {

    using Value = VoxelVolume::ValueType;
    using Eigen::Vector3f;
    using Eigen::Vector3i;
    using gsl::make_not_null;

    Vector3i const x0 = pos.array().floor().matrix().template cast<int>();
    Vector3i const x1 = pos.array().ceil().matrix().template cast<int>();
    Vector3f const xd = pos - x0.cast<float>();
    Vector3f const xn = Vector3f::Ones() - xd;

    // get values of lattice points
    Value const c000 = at(x0, ptr);
    Value const c001 = at(Vector3i{x0(0), x0(1), x1(2)}, ptr);
    Value const c010 = at(Vector3i{x0(0), x1(1), x0(2)}, ptr);
    Value const c011 = at(Vector3i{x0(0), x1(1), x1(2)}, ptr);
    Value const c100 = at(Vector3i{x1(0), x0(1), x0(2)}, ptr);
    Value const c101 = at(Vector3i{x1(0), x0(1), x1(2)}, ptr);
    Value const c110 = at(Vector3i{x1(0), x1(1), x0(2)}, ptr);
    Value const c111 = at(x1, ptr);

    auto const c00 = c000 * xn(0) + c100 * xd(0);
    auto const c01 = c001 * xn(0) + c101 * xd(0);
    auto const c10 = c010 * xn(0) + c110 * xd(0);
    auto const c11 = c011 * xn(0) + c111 * xd(0);

    auto const c0 = c00 * xn(1) + c10 * xd(1);
    auto const c1 = c01 * xn(1) + c11 * xd(1);

    auto const c = c0 * xn(2) + c1 * xd(2);

    return c;
  }

  VoxelVolume const &volume_;
};

class ModelSampler {

public:
  inline explicit ModelSampler(MeasurementModel const &model) noexcept
      : model_(model) {}

  template <class DerivedIn, class VectorOut>
  inline void operator()(Eigen::MatrixBase<DerivedIn> const &positions,
                         float offset, VectorOut &&values) const {
    using Eigen::Dynamic;
    using Eigen::Matrix;
    using Eigen::Vector3f;

    Expects(values.rows() == positions.rows());
    Expects(positions.cols() == 4);
    Expects(values.cols() == 1);

    Vector3f const min{model_.samplingRange.min, model_.densityRange.min,
                       model_.angleRange.min};
    Vector3f const scale{1.f / model_.samplingRange.stride,
                         1.f / model_.densityRange.stride,
                         1.f / model_.angleRange.stride};

    for (auto &&label : model_.labels()) {
      model_.withUnsafeDataPointer(label, [this, &positions, &values, min,
                                           scale, offset,
                                           label](float const *ptr) {
        for (Eigen::Index i = 0; i < positions.rows(); ++i) {
          if (static_cast<MeasurementModel::Label>(positions(i, 3)) != label) {
            continue;
          }

          Vector3f const position{positions(i, 0) + offset, positions(i, 1),
                                  positions(i, 2)};
          values(i) = this->interpolate(
              ((position - min).array() * scale.array()).matrix(), ptr);
        }
      });
    }
  }

  template <class Derived>
  inline Eigen::Matrix<float, Eigen::Dynamic, 1>
  operator()(Eigen::MatrixBase<Derived> const &positions, float offset) const {
    Eigen::Matrix<float, Eigen::Dynamic, 1> values(positions.rows());

    operator()(positions, offset, values);

    return values;
  }

private:
  inline float at(int x, int y, int z, Eigen::Vector3i const &sizeI,
                  float const *ptr) const {

    auto const index = x * sizeI.y() * sizeI.z() + z * sizeI.y() + y;
    return ptr[index];
  }

  template <class Derived>
  inline float interpolate(Eigen::MatrixBase<Derived> const &pos,
                           float const *ptr) const {
    using Eigen::Vector3f;
    using Eigen::Vector3i;
    using std::clamp;
    using std::log;

    Vector3i const sizeI{static_cast<int>(model_.samplingRange.numElements()),
                         static_cast<int>(model_.densityRange.numElements()),
                         static_cast<int>(model_.angleRange.numElements())};

    auto const x00 = clamp(static_cast<int>(pos(0)), 0, sizeI.x() - 1);
    auto const x01 = clamp(static_cast<int>(pos(2) + 0.5f), 0, sizeI.z() - 1);
    auto const x0 = static_cast<int>(pos(1));
    auto const x1 = x0 + 1;
    auto const xd = pos(1) - static_cast<float>(x0);
    auto const xn = 1.f - xd;

    auto const c0 = at(x00, clamp(x0, 0, sizeI.y() - 1), x01, sizeI, ptr);
    auto const c1 = at(x00, clamp(x1, 0, sizeI.y() - 1), x01, sizeI, ptr);

    auto const c = c0 * xn + c1 * xd;

    return log(c);
  }

  MeasurementModel const &model_;
};
#pragma clang diagnostic pop

} // namespace CortidQCT

max_stars_repo_name: ithron/CortidQCT

max_issues_repo_name: ithron/CortidQCT

max_forks_repo_name: ithron/CortidQCT


text:
#ifndef __GSL_PERMUTE_H__
#define __GSL_PERMUTE_H__

#include <gsl/gsl_permute_complex_long_double.h>
#include <gsl/gsl_permute_complex_double.h>
#include <gsl/gsl_permute_complex_float.h>

#include <gsl/gsl_permute_long_double.h>
#include <gsl/gsl_permute_double.h>
#include <gsl/gsl_permute_float.h>

#include <gsl/gsl_permute_ulong.h>
#include <gsl/gsl_permute_long.h>

#include <gsl/gsl_permute_uint.h>
#include <gsl/gsl_permute_int.h>

#include <gsl/gsl_permute_ushort.h>
#include <gsl/gsl_permute_short.h>

#include <gsl/gsl_permute_uchar.h>
#include <gsl/gsl_permute_char.h>

#endif /* __GSL_PERMUTE_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
#ifdef FUTILITY_HAVE_PETSC
#include <petsc.h>
#include <petscsys.h>
#include <petscfix.h>
#include <petscksp.h>
#if ((PETSC_VERSION_MAJOR<=3) && (PETSC_VERSION_MINOR<=5))
#include <petsc-private/fortranimpl.h>
#include <petsc-private/pcimpl.h>
#else
#include <petsc/private/fortranimpl.h>
#include <petsc/private/pcmgimpl.h>
#endif

#ifdef PETSC_USE_POINTER_CONVERSION
extern void *PetscToPointer(void*);
extern int PetscFromPointer(void *);
extern void PetscRmPointer(void*);
#else
#define PetscToPointer(a) (*(PetscFortranAddr *)(a))
#define PetscFromPointer(a) (PetscFortranAddr)(a)
#define PetscRmPointer(a)
#endif

#if ((PETSC_VERSION_MAJOR<=3) && (PETSC_VERSION_MINOR<=5))
//This is taken from petsc-3.5/src/ksp/pc/impls/mg/mgimpl.h
//There is no version of petsc/private/pcmgimpl.h that is includable in older
//  versions of PETSc, so we just have to manually redefine the typdefs defined
//  in mgimpl.h here.
typedef struct {
  PetscInt cycles;                             /* Type of cycle to run: 1 V 2 W */
  PetscInt level;                              /* level = 0 coarsest level */
  PetscInt levels;                             /* number of active levels used */
  Vec      b;                                  /* Right hand side */
  Vec      x;                                  /* Solution */
  Vec      r;                                  /* Residual */

  PetscErrorCode (*residual)(Mat,Vec,Vec,Vec);

  Mat           A;                             /* matrix used in forming residual*/
  KSP           smoothd;                       /* pre smoother */
  KSP           smoothu;                       /* post smoother */
  Mat           interpolate;
  Mat           restrct;                       /* restrict is a reserved word in C99 and on Cray */
  Vec           rscale;                        /* scaling of restriction matrix */
  PetscLogEvent eventsmoothsetup;              /* if logging times for each level */
  PetscLogEvent eventsmoothsolve;
  PetscLogEvent eventresidual;
  PetscLogEvent eventinterprestrict;
} PC_MG_Levels;
typedef struct {
  PCMGType  am;                               /* Multiplicative, additive or full */
  PetscInt  cyclesperpcapply;                 /* Number of cycles to use in each PCApply(), multiplicative only*/
  PetscInt  maxlevels;                        /* total number of levels allocated */
  PetscInt  galerkin;                         /* use Galerkin process to compute coarser matrices, 0=no, 1=yes, 2=yes but computed externally */
  PetscBool usedmfornumberoflevels;           /* sets the number of levels by getting this information out of the DM */

  PetscInt     nlevels;
  PC_MG_Levels **levels;
  PetscInt     default_smoothu;               /* number of smooths per level if not over-ridden */
  PetscInt     default_smoothd;               /*  with calls to KSPSetTolerances() */
  PetscReal    rtol,abstol,dtol,ttol;         /* tolerances for when running with PCApplyRichardson_MG */

  void          *innerctx;                   /* optional data for preconditioner, like PCEXOTIC that inherits off of PCMG */
  PetscLogStage stageApply;
} PC_MG;
#endif

#ifdef PETSC_HAVE_FORTRAN_CAPS
#define pcmgsoftreset_ PCMGSOFTRESET
#elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) && !defined(FORTRANDOUBLEUNDERSCORE)
#define pcmgsoftreset_ pcmgsoftreset
#endif

#undef __FUNCT__
#define __FUNCT__ "PCMGSoftReset"

PetscErrorCode PCMGSoftReset(PC pc)
{
  PC_MG          *mg        = (PC_MG*)pc->data;
  PC_MG_Levels   **mglevels = mg->levels;
  PetscErrorCode ierr;
  PetscInt       i,n;

  PetscFunctionBegin;
  if (mglevels) {
    n = mglevels[0]->levels;
    for (i=0; i<n; i++) {
      ierr = MatDestroy(&mglevels[i]->A);CHKERRQ(ierr);
      if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
        ierr = KSPReset(mglevels[i]->smoothd);CHKERRQ(ierr);
      }
      ierr = KSPReset(mglevels[i]->smoothu);CHKERRQ(ierr);
    }
  }
  //This may not work if DM's are used to define the MG structure.
  pc->setupcalled = PETSC_FALSE;
  PetscFunctionReturn(0);
}

PETSC_EXTERN void PETSC_STDCALL  pcmgsoftreset_(PC pc, int *__ierr ){
*__ierr = PCMGSoftReset(
  (PC)PetscToPointer((pc) ));
}
#endif

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
/* cdf/poisson.c
 *
 * Copyright (C) 2004 Jason H. Stover.
 *
 * 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.
 */

/*
 * Computes the cumulative distribution function for a Poisson
 * random variable. For a Poisson random variable X with parameter
 * mu,
 *
 *          Pr( X <= k ) = Pr( Y >= p )
 *
 * where Y is a gamma random variable with parameters k+1 and 1.
 *
 * Reference:
 *
 * W. Feller, "An Introduction to Probability and Its
 * Applications," volume 1. Wiley, 1968. Exercise 46, page 173,
 * chapter 6.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>

#include "error.h"

/*
 * Pr (X <= k) for a Poisson random variable X.
 */

double
gsl_cdf_poisson_P (const unsigned int k, const double mu)
{
  double P;
  double a;

  if (mu <= 0.0)
    {
      CDF_ERROR ("mu <= 0", GSL_EDOM);
    }

  a = (double) k + 1.0;
  P = gsl_cdf_gamma_Q (mu, a, 1.0);

  return P;
}

/*
 * Pr ( X > k ) for a Possion random variable X.
 */

double
gsl_cdf_poisson_Q (const unsigned int k, const double mu)
{
  double Q;
  double a;

  if (mu <= 0.0)
    {
      CDF_ERROR ("mu <= 0", GSL_EDOM);
    }

  a = (double) k + 1.0;
  Q = gsl_cdf_gamma_P (mu, a, 1.0);

  return Q;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include "scripts.h"
#include "utility/fileutils.h"
#include "lua.h"
#include "lauxlib.h"
#include "subsystems.h"
#include "rendering/rendering.h"
#include <gsl/gsl_math.h>

/**
 * Retrieve the integer value in the table at the top of the stack
 * associated with the specified key.
 */
sd_result get_int_field(lua_State *L, const char *key, int *out) {
  lua_getfield(L, -1, key);
  if (!lua_isnumber(L, -1)) {
    return SD_SCRIPT_TYPE_ERROR;
  }
  *out = (int)lua_tonumber(L, -1);
  lua_pop(L, 1);
  return SD_OK;
}

/**
 * Retrieve the double value in the table at the top of the stack
 * associated with the specified key.
 */
sd_result get_double_field(lua_State *L, const char *key, double *out) {
  lua_getfield(L, -1, key);
  if (!lua_isnumber(L, -1)) {
    return SD_SCRIPT_TYPE_ERROR;
  }
  *out = lua_tonumber(L, -1);
  lua_pop(L, 1);
  return SD_OK;
}

/**
 * Retrieve the string value in the table at the top of the stack
 * associated with the specified key.
 */
sd_result get_string_field(lua_State *L, const char *key, char *out, size_t out_size) {
  lua_getfield(L, -1, key);
  if (!lua_isstring(L, -1)) {
    return SD_SCRIPT_TYPE_ERROR;
  }
  size_t max = GSL_MIN(lua_rawlen(L, -1), out_size);
  strncpy(out, lua_tostring(L, -1), max);
  lua_pop(L, 1);
  return SD_OK;
}

/**
 * Retrieve the boolean value in the table at the top of the stack
 * associated with the specified key.
 */
sd_result get_bool_field(lua_State *L, const char *key, char *out) {
  lua_getfield(L, -1, key);
  if (!lua_isboolean(L, -1)) {
    return SD_SCRIPT_TYPE_ERROR;
  }
  *out = lua_toboolean(L, -1);
  lua_pop(L, 1);
  return SD_OK;
}

sd_result read_init_info(const char *resource_path, initialization_params *out) {
  lua_State *L = luaL_newstate();
  bzero(out, sizeof(initialization_params));
  window_creation_info *winfo = calloc(1, sizeof(window_creation_info));
  out->window_info = winfo;

  char *file_path = NULL;
  fileutils.get_full_path(resource_path, &file_path);
  if (luaL_loadfile(L, file_path) || lua_pcall(L, 0, 0, 0)) {
    printf("Problem loading file %s\n", file_path);
    return SD_GEN_ERROR;
  }
  lua_getglobal(L, "timer");
  if (!lua_isstring(L, -1)) {
    printf("Timer type is not a string\n");
    return SD_SCRIPT_TYPE_ERROR;
  }
  const char *timer_type = lua_tostring(L,-1);
  if (strncasecmp("glfw", timer_type, lua_rawlen(L, -1))) {
    printf("glfw timer specified\n");
  };
  out->timer_function = glfwGetTime;
  lua_pop(L, 1);

  lua_getglobal(L, "window");
  if (!lua_istable(L, -1)) {
    printf("Window description is not a table\n");
    return SD_SCRIPT_TYPE_ERROR;
  }

  get_int_field(L, "w", &winfo->width);
  get_int_field(L, "h", &winfo->height);
  get_bool_field(L, "full_screen", &winfo->full_screen);
  get_string_field(L, "title", winfo->window_title,
                  sizeof(winfo->window_title));

  // Sub table
  lua_getfield(L, -1, "clear_color");
  double dbuffer;
  get_double_field(L, "R", &dbuffer);
  winfo->clear_color[0] = (GLfloat)dbuffer;
  get_double_field(L, "G", &dbuffer);
  winfo->clear_color[1] = (GLfloat)dbuffer;
  get_double_field(L, "B", &dbuffer);
  winfo->clear_color[2] = (GLfloat)dbuffer;
  get_double_field(L, "A", &dbuffer);
  winfo->clear_color[3] = (GLfloat)dbuffer;
  lua_pop(L, 1);

  return SD_OK;
}

scripts_struct const scripts = {
  .read_init_info=read_init_info
};


text:
#ifndef CWANNIER_PARSESCF_H
#define CWANNIER_PARSESCF_H

#include <stdbool.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include "bstrlib/bstrlib.h"

#define CWANNIER_PARSESCF_OK 0
#define CWANNIER_PARSESCF_ERR 1

int ParseSCF(char *filePath, double *num_electrons, double *alat, gsl_matrix *R);

void get_b_row(gsl_matrix *R, int row, bstring line);

#endif // CWANNIER_PARSESCF_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
#pragma once
#include <memory>
#include "halley/resources/resource.h"
#include "halley/utils/utils.h"
#include <gsl/gsl>
#include "halley/core/resources/resource_collection.h"

namespace Halley
{
	class ResourceLoader;

	class BinaryFile : public Resource
	{
	public:
		BinaryFile();
		explicit BinaryFile(const Bytes& data);
        explicit BinaryFile(Bytes&& data);
		explicit BinaryFile(gsl::span<const gsl::byte> data);
		explicit BinaryFile(std::unique_ptr<ResourceDataStream> stream);

		static std::unique_ptr<BinaryFile> loadResource(ResourceLoader& loader);
		constexpr static AssetType getAssetType() { return AssetType::BinaryFile; }
		void reload(Resource&& resource) override;

		const Bytes& getBytes() const;
		Bytes& getBytes();
		gsl::span<const gsl::byte> getSpan() const;
		gsl::span<gsl::byte> getSpan();

		std::shared_ptr<ResourceDataStream> getStream() const;

	private:
		Bytes data;
		std::shared_ptr<ResourceDataStream> stream;
		bool streaming = false;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
/*
 *  lowly.h
 *  low-ell likelihood in the pixel domain
 *
 *  JF Cardoso and Simon Prunet.  January 2008.
 *
 */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_int.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>

#ifdef HAS_PMC
#include "invgamma.h"
#endif
#include "lowly_common.h"
#include "chealpix.h"

#ifndef __LOWLY__
#define __LOWLY__

// llll = low ell log-likelihood
typedef struct  {
  long npix_seen ;
  long npix,nside ;
  double * X_seen  ;
  double * F_seen  ;
  double * NoiseVar ;
  double * CosMat ;
  long lmax ;
  double * Wl  ;
  double * q_buf;
  double * X_seen_buf;
  double * F_seen_buf;
  double * covmat_buf;
  double * buffer;
  int isLog;
  int partial_nindex;
  int* partial_index;
  double *partial_cl;
  double v2;
  STRUCT_TIMER_DEFS;
} llll_struct ;

typedef llll_struct llll_partial_struct;

// create and initialize the structure holding the material to compute the likelihood
llll_struct* init_lowly(
                        double *X,          // map size npix
                        unsigned char *Mask,       // mask size npix
                        double* Foreground, // if !=NULL a foreground map
                        long npix,
                        char* ordering,     // ring or nested
                        double noisevar,    // microK/pixel
                        double v2,          // sigma^2 for foreground marginalization
                        double *Wl,         // wl Cl unit lmax+1 values
                        double lmax,
                        error **err );

void free_lowly(llll_struct** self);

// return the low-ell log-likelihood for spectrum C_l
double lowly_lkl(
    void * llll_str ,
    double * C_l ,
    error **err
    ) ;

// return the low-ell log-likelihood for log(Cl)
double lowly_lkl_logcl(
    void * llll_str,
    double * logCl,
    error **err
    );

llll_struct * init_lowly_log(llll_struct* self,error **err);

// These are wrappers on top of lowly_lkl and lowly_lkl_logcl
// in case only a subset of Cls are varied as parameters
double lowly_partial_lkl(void *llll_str, double *cl, error **err);
double lowly_partial_lkl_logcl(void *llll_str, double *logcl, error **err);

void free_lowly_partial(llll_partial_struct **elf);
llll_partial_struct* init_lowly_partial( llll_struct *orig,int n,int* index, double* cl,error **err) ;
llll_partial_struct* init_lowly_partial_logcl( llll_struct *orig,int n,int* index, double* cl,error **err) ;
double lowly_lkl_with_noise(void *llll_str, double *cl, error **err);

// ------ internals.  Used for computation ---------------------

// assign and compute the npix_seen * npix_seen matrix of the cosines between seen pixels
double * build_cosine_matrix (
    const long nside,
    const long * pixel_indices ,
    const long npix_seen ,
    const int ordering,
    error **err
    ) ;

// assign and compute \sum_l W_l C_l P_l(cos) (2l+1)/4pi
double * build_cov_matrix (
    double *orig,
    const double * CosMat ,
    const double * NoiseVar,
    const long npix_seen ,
    const long lmax ,
    const double * q,
   error **err
    ) ;

// assign and compute inverse covariance matrix, same inputs as above
double * build_invcov_matrix (
    double *orig,
    const double * cosmat,
    const double * noisevar,
    const long npix_seen,
    const long lmax,
    const double *q,
    error **err
    ) ;

// assign and compute W_l P_l(cos) (2l+1)/4pi
double * build_dcov_matrix (
    double *orig,
    const double *cosmat,
    const long npix_seen,
    const long lmax,
    const double wl,
    error **err
    );

// compute signal covariance for one value of z=cos(theta_ij)
double scalCov (
    const double z,
    const long lmax,
    const double *q
    ) ;

// compute signal covariance derivative (with respect to C_l) for one value of z=cos(theta_ij)
double scaldCov (
    const double z,
    const double wl, // = Wl[ell]
    const long ell
    ) ;

// tentative computation of signal covariance matrix using matrix (blas) algebra
double * matCov(
    const double * cosmat,
    const long npix_seen,
    const long lmax,
    const double * q);

// Computes gradient of likelihood with respect to C_l
double * grad_lowly_lkl (
    double * orig,
    llll_struct * self,
    const double * covmat,
    const double * invcovmat,
    error **err
    ) ;

// Computes diagonal of fisher matrix, with respect to C_l
double * diag_fisher_lowly_lkl (
    double *orig,
    llll_struct * self,
    //const double * invcovmat,
    const double *q,
    error **err) ;

// Computes peak of posterior by pseudo-Newton iteration
double * pseudo_newton(
    llll_struct * self,
    double * cl,
    long niter,
    error **err
    ) ;

#ifdef HAS_PMC
typedef struct {
  invgamma * _ivg;
  double *data,*wl,*N;
  size_t ndim;
} invgamma_prop;

invgamma_prop* invgamma_prop_init(size_t ndim,int* ell,double* clhat, double* fsky, double* wl,double noisevar,long nside,error **err);
void free_invgamma_prop(invgamma_prop **pelf);
double invgamma_prop_log_pdf(invgamma_prop *self, double *x, error **err);
long simulate_invgamma_prop(pmc_simu *psim,invgamma_prop *proposal, gsl_rng * r,parabox *pb,error **err);
size_t get_importance_weight_invgamma_prop(pmc_simu *psim, const invgamma_prop *m,
                                           posterior_log_pdf_func *posterior_log_pdf,
                                           void *extra, error **err);
void update_invgamma_prop(invgamma_prop *g,pmc_simu *psim,error **err);
invgamma_prop* invgamma_prop_init_new(size_t ndim,int* ell,double* clhat, double* fsky, double* wl,double* N,error **err);
#endif

void computeAlphaBeta(double *alpha, double *beta,int ell, double clhat, double fsky, double wl, double Nl);
double cl2xl(double cl, double wl, double Nl);

#define lowly_base           -10 + clowly_base
#define lowly_outofrange     -2 + lowly_base
#define lowly_negcl          -4 + lowly_base
#define lowly_fish_negative  -5 + lowly_base
#define lowly_spec_negative  -6 + lowly_base

#endif


text:
/* See compilation notes in simple_cplx.c*/
#include <complex.h> //nice names for C's complex types
#include <gsl/gsl_vector.h> //gsl_vector_complex

gsl_vector_complex *cvec_dot_gslcplx(gsl_vector_complex *v, gsl_complex x);
gsl_vector_complex *vec_dot_gslcplx(gsl_vector *v, gsl_complex x);
gsl_vector_complex *cvec_dot_c(gsl_vector_complex *v, complex double x);
gsl_vector_complex *vec_dot_c(gsl_vector *v, complex double x);
void gsl_vector_complex_print(gsl_vector_complex *v);

#define gsl_cplx_from_c99(x) (gsl_complex){.dat= {creal(x), cimag(x)}}

complex double ddot (complex double x, complex double y);

#define dot(x,y)  _Generic((x),                                 \
                    gsl_vector*: dot_given_vec(y),              \
                    gsl_vector_complex*: dot_given_cplx_vec(y), \
                    default: ddot)((x),(y))

#define dot_given_vec(y) _Generic((y),                          \
                        gsl_complex: vec_dot_gslcplx,           \
                        default: vec_dot_c)

#define dot_given_cplx_vec(y) _Generic((y),                     \
                        gsl_complex: cvec_dot_gslcplx,          \
                        default: cvec_dot_c)

max_stars_repo_name: catalingheorghe/learning-activities

max_issues_repo_name: catalingheorghe/learning-activities

max_forks_repo_name: catalingheorghe/learning-activities


text:
#ifndef SCANMATCHERPROCESSOR_H
#define SCANMATCHERPROCESSOR_H

#include <gmapping/log/sensorlog.h>
#include <gmapping/sensor/sensor_range/rangesensor.h>
#include <gmapping/sensor/sensor_range/rangereading.h>
//#include <gsl/gsl_eigen.h>
#include <gmapping/scanmatcher/scanmatcher.h>

namespace GMapping {

class ScanMatcherProcessor
{
	public:
  ScanMatcherProcessor(const ScanMatcherMap& m);
  ScanMatcherProcessor (double xmin, double ymin, double xmax, double ymax, double delta, double patchdelta);
		virtual ~ScanMatcherProcessor ();
		virtual void processScan(const RangeReading & reading);
		void setSensorMap(const SensorMap& smap, std::string sensorName="FLASER");
		void init();
		void setMatchingParameters
			(double urange, double range, double sigma, int kernsize, double lopt, double aopt, int iterations, bool computeCovariance=false);
		void setRegistrationParameters(double regScore, double critScore);
		OrientedPoint getPose() const;
		inline const ScanMatcherMap& getMap() const {return m_map;}
		inline ScanMatcher& matcher() {return m_matcher;}
		inline void setmaxMove(double mmove){m_maxMove=mmove;}
		bool useICP;
	protected:
		ScanMatcher m_matcher;
		bool m_computeCovariance;
		bool m_first;
		SensorMap m_sensorMap;
		double m_regScore, m_critScore;
		unsigned int m_beams;
		double m_maxMove;
		//state
		ScanMatcherMap m_map;
		OrientedPoint m_pose;
		OrientedPoint m_odoPose;
		int  m_count;
		//gsl_eigen_symmv_workspace * m_eigenspace;
};

};

#endif

max_stars_repo_name: HouCongCN/Gmapping-Code

max_issues_repo_name: HouCongCN/Gmapping-Code

max_forks_repo_name: HouCongCN/Gmapping-Code


text:
#include <math.h>
#include <stdio.h>
#include <complex.h>

#include <lapacke.h>
#include <libhades.h>

#include "unittest.h"

int test_trace()
{
    matrix_t *A, *B;
    matrix_complex_t *C, *D;
    double trace;
    complex_t trace_cplx;
    unittest_t test;
    unittest_init(&test, "matrix_mult", "Test trace");

    A = matrix_alloc(3,3);
    B = matrix_alloc(3,3);

    matrix_set(A, 0,0, 1);
    matrix_set(A, 0,1, 2);
    matrix_set(A, 0,2, 3);
    matrix_set(A, 1,0, 4);
    matrix_set(A, 1,1, 5);
    matrix_set(A, 1,2, 6);
    matrix_set(A, 2,0, 7);
    matrix_set(A, 2,1, 8);
    matrix_set(A, 2,2, 9);

    matrix_set(B, 0,0, 9);
    matrix_set(B, 0,1, 8);
    matrix_set(B, 0,2, 7);
    matrix_set(B, 1,0, 6);
    matrix_set(B, 1,1, 5);
    matrix_set(B, 1,2, 4);
    matrix_set(B, 2,0, 3);
    matrix_set(B, 2,1, 2);
    matrix_set(B, 2,2, 1);

    trace = matrix_trace_AB(A, B);
    AssertEqual(&test, trace, 189);

    AssertEqual(&test, matrix_trace(A), 15);
    AssertEqual(&test, matrix_trace(B), 15);

    matrix_free(A);
    matrix_free(B);

    C = matrix_complex_alloc(3,3);
    D = matrix_complex_alloc(3,3);

    matrix_set(C, 0,0, CPLX(0,1));
    matrix_set(C, 0,1, CPLX(-2,0));
    matrix_set(C, 0,2, CPLX(3,1));
    matrix_set(C, 1,0, CPLX(-2,0));
    matrix_set(C, 1,1, CPLX(1,0));
    matrix_set(C, 1,2, CPLX(0,1));
    matrix_set(C, 2,0, CPLX(0,-2));
    matrix_set(C, 2,1, CPLX(1,3));
    matrix_set(C, 2,2, CPLX(5,0));

    matrix_set(D, 0,0, CPLX(0,-2));
    matrix_set(D, 0,1, CPLX(0,1));
    matrix_set(D, 0,2, CPLX(2,-2));
    matrix_set(D, 1,0, CPLX(1,0));
    matrix_set(D, 1,1, CPLX(5,0));
    matrix_set(D, 1,2, CPLX(0,1));
    matrix_set(D, 2,0, CPLX(-1,0));
    matrix_set(D, 2,1, CPLX(-4,-4));
    matrix_set(D, 2,2, CPLX(0,1));

    trace_cplx = matrix_complex_trace(C);
    AssertEqual(&test, CREAL(trace_cplx), 6);
    AssertEqual(&test, CIMAG(trace_cplx), 1);
    trace_cplx = matrix_complex_trace(D);
    AssertEqual(&test, CREAL(trace_cplx), 5);
    AssertEqual(&test, CIMAG(trace_cplx), -1);

    trace = matrix_trace_complex_AB_real(C, D);
    AssertEqual(&test, trace, -1);

    trace_cplx = matrix_trace_complex_AB(C, D);
    AssertEqual(&test, CREAL(trace_cplx), -1);
    AssertEqual(&test, CIMAG(trace_cplx), -5);

    matrix_complex_free(C);
    matrix_complex_free(D);

    return test_results(&test, stderr);
}

int test_matrix_mult()
{
    matrix_t *A, *B, *C;
    unittest_t test;
    unittest_init(&test, "matrix_mult", "Test matrix multiplication");

    A = matrix_alloc(2,3);
    B = matrix_alloc(3,2);

    matrix_set(A, 0,0, 3);
    matrix_set(A, 0,1, 2);
    matrix_set(A, 0,2, 1);
    matrix_set(A, 1,0, 1);
    matrix_set(A, 1,1, 0);
    matrix_set(A, 1,2, 2);

    matrix_set(B, 0,0, 1);
    matrix_set(B, 0,1, 2);
    matrix_set(B, 1,0, 0);
    matrix_set(B, 1,1, 1);
    matrix_set(B, 2,0, 4);
    matrix_set(B, 2,1, 0);

    C = matrix_mult(A,B,1,NULL);

    matrix_free(A);
    matrix_free(B);

    AssertEqual(&test, matrix_get(C,0,0), 7);
    AssertEqual(&test, matrix_get(C,0,1), 8);
    AssertEqual(&test, matrix_get(C,1,0), 9);
    AssertEqual(&test, matrix_get(C,1,1), 2);

    matrix_free(C);

    return test_results(&test, stderr);
}

int test_matrix_zeros()
{
    int m,n;
    matrix_t *A;
    matrix_complex_t *C;
    unittest_t test;

    unittest_init(&test, "matrix_zeros", "Test function matrix_zeros");

    for(m = 6; m < 20; m++)
        for(n = 1; n < 8; n++)
        {
            int im, in;
            A = matrix_zeros(m,n, NULL);
            C = matrix_complex_zeros(m,n, NULL);
            for(im = 0; im < m; im++)
                for(in = 0; in < n; in++)
                {
                    complex_t c = matrix_get(C,im,in);
                    AssertEqual(&test, matrix_get(A,im,in),0);
                    AssertEqual(&test, CREAL(c), 0);
                    AssertEqual(&test, CIMAG(c), 0);
                }

            matrix_free(A);
            matrix_complex_free(C);
        }

    return test_results(&test, stderr);
}

int test_zeros()
{
    matrix_t *A;
    matrix_complex_t *B;
    unittest_t test;
    unittest_init(&test, "matrix_zeros", "Test zeros");
    size_t rows = 100, columns = 150;

    A = matrix_zeros        (rows, columns, NULL);
    B = matrix_complex_zeros(rows, columns, NULL);

    for(size_t i = 0; i < rows; i++)
        for(size_t j = 0; j < columns; j++)
        {
            AssertEqual(&test, matrix_get(A, i,j), 0);
            AssertEqual(&test, matrix_get(B, i,j), 0);
        }

    matrix_free(A);
    matrix_complex_free(B);

    return test_results(&test, stderr);
}

int test_matrix_kron()
{
    matrix_t *A, *B, *AB;
    matrix_complex_t *C, *D, *CD;
    matrix_t *E, *EE;
    unittest_t test;
    unittest_init(&test, "matrix_kron", "Test function matrix_kron");

    A = matrix_alloc(3,2);
    B = matrix_alloc(2,2);

    matrix_set(A, 0,0, 1);
    matrix_set(A, 0,1, 2);
    matrix_set(A, 1,0, 3);
    matrix_set(A, 1,1, 4);
    matrix_set(A, 2,0, 5);
    matrix_set(A, 2,1, 6);

    matrix_set(B, 0,0, 7);
    matrix_set(B, 0,1, 8);
    matrix_set(B, 1,0, 9);
    matrix_set(B, 1,1, 0);

    AB = matrix_kron(A, B, NULL);

    matrix_free(A);
    matrix_free(B);

    AssertEqual(&test, matrix_get(AB, 0,0), 7);
    AssertEqual(&test, matrix_get(AB, 0,1), 8);
    AssertEqual(&test, matrix_get(AB, 0,2), 14);
    AssertEqual(&test, matrix_get(AB, 0,3), 16);

    AssertEqual(&test, matrix_get(AB, 1,0), 9);
    AssertEqual(&test, matrix_get(AB, 1,1), 0);
    AssertEqual(&test, matrix_get(AB, 1,2), 18);
    AssertEqual(&test, matrix_get(AB, 1,3), 0);

    AssertEqual(&test, matrix_get(AB, 2,0), 21);
    AssertEqual(&test, matrix_get(AB, 2,1), 24);
    AssertEqual(&test, matrix_get(AB, 2,2), 28);
    AssertEqual(&test, matrix_get(AB, 2,3), 32);

    AssertEqual(&test, matrix_get(AB, 3,0), 27);
    AssertEqual(&test, matrix_get(AB, 3,1), 0);
    AssertEqual(&test, matrix_get(AB, 3,2), 36);
    AssertEqual(&test, matrix_get(AB, 3,3), 0);

    AssertEqual(&test, matrix_get(AB, 4,0), 35);
    AssertEqual(&test, matrix_get(AB, 4,1), 40);
    AssertEqual(&test, matrix_get(AB, 4,2), 42);
    AssertEqual(&test, matrix_get(AB, 4,3), 48);

    AssertEqual(&test, matrix_get(AB, 5,0), 45);
    AssertEqual(&test, matrix_get(AB, 5,1), 0);
    AssertEqual(&test, matrix_get(AB, 5,2), 54);
    AssertEqual(&test, matrix_get(AB, 5,3), 0);

    matrix_free(AB);

    C = matrix_complex_alloc(2,2);
    D = matrix_complex_alloc(2,2);

    matrix_set(C, 0,0, CPLX(1,1));
    matrix_set(C, 0,1, CPLX(2,-2));
    matrix_set(C, 1,0, CPLX(5,0));
    matrix_set(C, 1,1, CPLX(0,-1));

    matrix_set(D, 0,0, CPLX(0,-1));
    matrix_set(D, 0,1, CPLX(2,0));
    matrix_set(D, 1,0, CPLX(5,0));
    matrix_set(D, 1,1, CPLX(0,7));

    CD = matrix_complex_alloc(4,4);
    matrix_complex_kron(C, D, CD);

    matrix_complex_free(C);
    matrix_complex_free(D);

    AssertEqual(&test, CREAL(matrix_get(CD, 0,0)), 1);
    AssertEqual(&test, CIMAG(matrix_get(CD, 0,0)), -1);
    AssertEqual(&test, CREAL(matrix_get(CD, 0,1)), 2);
    AssertEqual(&test, CIMAG(matrix_get(CD, 0,1)), 2);
    AssertEqual(&test, CREAL(matrix_get(CD, 0,2)), -2);
    AssertEqual(&test, CIMAG(matrix_get(CD, 0,2)), -2);
    AssertEqual(&test, CREAL(matrix_get(CD, 0,3)), 4);
    AssertEqual(&test, CIMAG(matrix_get(CD, 0,3)), -4);

    AssertEqual(&test, CREAL(matrix_get(CD, 1,0)), 5);
    AssertEqual(&test, CIMAG(matrix_get(CD, 1,0)), 5);
    AssertEqual(&test, CREAL(matrix_get(CD, 1,1)), -7);
    AssertEqual(&test, CIMAG(matrix_get(CD, 1,1)), 7);
    AssertEqual(&test, CREAL(matrix_get(CD, 1,2)), 10);
    AssertEqual(&test, CIMAG(matrix_get(CD, 1,2)), -10);
    AssertEqual(&test, CREAL(matrix_get(CD, 1,3)), 14);
    AssertEqual(&test, CIMAG(matrix_get(CD, 1,3)), 14);

    AssertEqual(&test, CREAL(matrix_get(CD, 2,0)), 0);
    AssertEqual(&test, CIMAG(matrix_get(CD, 2,0)), -5);
    AssertEqual(&test, CREAL(matrix_get(CD, 2,1)), 10);
    AssertEqual(&test, CIMAG(matrix_get(CD, 2,1)), 0);
    AssertEqual(&test, CREAL(matrix_get(CD, 2,2)), -1);
    AssertEqual(&test, CIMAG(matrix_get(CD, 2,2)), 0);
    AssertEqual(&test, CREAL(matrix_get(CD, 2,3)), 0);
    AssertEqual(&test, CIMAG(matrix_get(CD, 2,3)), -2);

    AssertEqual(&test, CREAL(matrix_get(CD, 3,0)), 25);
    AssertEqual(&test, CIMAG(matrix_get(CD, 3,0)), 0);
    AssertEqual(&test, CREAL(matrix_get(CD, 3,1)), 0);
    AssertEqual(&test, CIMAG(matrix_get(CD, 3,1)), 35);
    AssertEqual(&test, CREAL(matrix_get(CD, 3,2)), 0);
    AssertEqual(&test, CIMAG(matrix_get(CD, 3,2)), -5);
    AssertEqual(&test, CREAL(matrix_get(CD, 3,3)), 7);
    AssertEqual(&test, CIMAG(matrix_get(CD, 3,3)), 0);

    matrix_complex_free(CD);

    E = matrix_zeros(2,2, NULL);

    matrix_set(E, 0,0, 1);
    matrix_set(E, 1,1, 2);

    EE = matrix_kron(E, E, NULL);
    matrix_free(E);

    AssertEqual(&test, matrix_get(EE, 0,0), 1);
    AssertEqual(&test, matrix_get(EE, 1,1), 2);
    AssertEqual(&test, matrix_get(EE, 2,2), 2);
    AssertEqual(&test, matrix_get(EE, 3,3), 4);

    for(size_t i = 0; i < 4; i++)
        for(size_t j = 0; j < 4; j++)
            if(i != j)
                AssertEqual(&test, matrix_get(EE, i,j), 0);

    matrix_free(EE);

    return test_results(&test, stderr);
}

int main(int argc, char *argv[])
{
    test_matrix_zeros();
    test_matrix_kron();
    test_matrix_mult();
    test_zeros();
    test_trace();

    return 0;
}

max_stars_repo_name: michael-hartmann/libhades

max_issues_repo_name: michael-hartmann/libhades

max_forks_repo_name: michael-hartmann/libhades


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f1(double x, void* params)
{
    return x*x;
}

double f2(double x, void* params)
{
    return 1 / sqrt(x);
}

void calculateIntegralVoid(double (*f)(double x), double a, double b, int maxInterval, double expected)
{
    FILE* output = fopen("integral.txt", "w");
    int n = 1;
    while (n < maxInterval)
    {
        double dx = (b - a) / (double)n;
        double sum = 0;

        for (int i = 1; i < n; ++i)
        {
            sum += f(a + i * dx);
        }
        sum *= dx;
        fprintf (output,"%d %f\n", n, expected - sum);

        n++;
    }
}

int calculateIntegral(double (*f)(double x), double a, double b, int n, double expected, double error)
{
    double dx = (b - a) / (double)n;
    double sum = 0;

    for (int i = 1; i < n; ++i)
    {
        sum += f(a + i * dx);
    }
    sum *= dx;
    if (fabs(sum - expected) < error)
    {
        printf("Sum with error %.1e error for %d samples: %.18f\n", error, n, sum);
        return 1;
    }
    printf("DIFFERENCE with error %.1e error for %d samples: %.18f\n", error, n, fabs(sum - expected));
    return 1;
}

void tryFindIntergral(double (*f)(double x), double a, double b, double expected, double error, int step)
{
    int interval = step;
    while (!calculateIntegral(f, a, b, interval, expected, error))
    {
        interval += step;
    }
}

void gslIntegral(double (*f)(double x, void* params), double expected)
{
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

    double result, error;
    double alpha = 1.0;

    gsl_function F;
    F.function = f;
    F.params = α

    gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000, w, &result, &error);

    printf ("result          = % .18f\n", result);
    printf ("exact result    = % .18f\n", expected);
    printf ("estimated error = % .18f\n", error);
    printf ("actual error    = % .18f\n", result - expected);
    printf ("intervals       = %zu\n", w->size);

    gsl_integration_workspace_free (w);
}

int main (int argc, char* argv[])
{
    // tryFindIntergral(f2, 0, 1.0, 2.0, 1e-3, 1000);
    // tryFindIntergral(f2, 0, 1.0, 2.0, 1e-4, 215000000);
    // tryFindIntergral(f2, 0, 1.0, 2.0, 1e-5, 100000);
    // tryFindIntergral(f2, 0, 1.0, 2.0, 1e-6, 1000000);

    // calculateIntegralVoid(f2, 0.0, 1.0, 100000, 2.0);

    gslIntegral(f2, 2.0);
    return 0;
}


text:
/******************************************************************************
 *                                                                            *
 * PROBLEM.C                                                                  *
 *                                                                            *
 * INITIAL CONDITIONS FOR BONDI INFLOW                                        *
 *                                                                            *
 ******************************************************************************/

#include "decs.h"

// Local functions
void coord_transform(double *Pr, int i, int j);
double lfish_calc(double rmax) ;

void set_problem_params()
{
}

// Rootfinding for analytic Bondi solution
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
struct params {
  double r, mdot, n, C1, C2;
};
double rfunc(double T, void *params)
{
  struct params *p = (struct params *) params;
  //double r = p->r;
  double r = MY_MAX(2.5, p->r); // Solution breaks inside event horizon
  double mdot = p->mdot;
  double n = p->n;
  double C1 = p->C1;
  double C2 = p->C2;

  double resid = (pow(1. + (1. + n)*T,2)*(1. - 2.*mdot/r + pow(C1,2)/
         (pow(r,4)*pow(T,2*n))) - C2);

  //printf("r = %e T = %e resid = %e\n", r, T, resid);
  return resid;
  //return (pow(1. + (1. + n)*T,2)*(1. - 2.*mdot/r + pow(C1,2)/
  //       (pow(r,4)*pow(T,2*n))) - C2);
}

double C1, C2, T_bondi[N1+2*NG][N2+2*NG], n;

// Adapted from M. Chandra
double get_Tfunc(double T, double r)
{
  return pow(1.+(1.+n)*T,2.)*(1.-2./r+pow(C1/r/r/pow(T,n),2.))-C2;
}

double get_T(double r)
{
  double rtol = 1.e-12;
  double ftol = 1.e-14;
  double Tmin = 0.6*(sqrt(C2) - 1.)/(n + 1);
  double Tmax = pow(C1*sqrt(2./r/r/r),1./n);
  double f0, f1, fh;
  double T0, T1, Th;
  T0 = 0.6*Tmin;
  f0 = get_Tfunc(T0, r);
  T1 = Tmax;
  f1 = get_Tfunc(T1, r);

  if (f0*f1 > 0.) {
    printf("Failed solving for T at r = %e C1 = %e C2 = %e\n", r, C1, C2);
    exit(-1);
  }

  Th = (f1*T0 - f0*T1)/(f1 - f0);
  fh = get_Tfunc(Th, r);
  double epsT = rtol*(Tmin + Tmax);
  while (fabs(Th - T0) > epsT && fabs(Th - T1) > epsT && fabs(fh) > ftol) {
    if (fh*f0 < 0.) {
      T0 = Th;
      f0 = fh;
    } else {
      T1 = Th;
      f1 = fh;
    }

    Th = (f1*T0 - f0*T1)/(f1 - f0);
    fh = get_Tfunc(Th, r);
  }

  return Th;
}

void bl_to_ks(double X[NDIM], double ucon_bl[NDIM], double ucon_ks[NDIM])
{
  double r, th;
  bl_coord(X, &r, &th);

  double trans[NDIM][NDIM];
  DLOOP2 trans[mu][nu] = 0.;
  DLOOP1 trans[mu][mu] = 1.;
  trans[0][1] = 2.*r/(r*r - 2.*r + a*a);
  trans[3][1] = a/(r*r - 2.*r + a*a);

  DLOOP1 ucon_ks[mu] = 0.;
  DLOOP2 ucon_ks[mu] += trans[mu][nu]*ucon_bl[nu];
}

void fourvel_to_prim(double ucon[NDIM], double prim[NVAR],
  struct of_geom *geom)
{
  double alpha, beta[NDIM], gamma;
  alpha = 1.0/sqrt(-geom->gcon[0][0]);
  beta[1] = alpha*alpha*geom->gcon[0][1];
  beta[2] = alpha*alpha*geom->gcon[0][2];
  beta[3] = alpha*alpha*geom->gcon[0][3];
  gamma = ucon[0]*alpha;

  prim[U1] = ucon[1] + beta[1]*gamma/alpha;
  prim[U2] = ucon[2] + beta[2]*gamma/alpha;
  prim[U3] = ucon[3] + beta[3]*gamma/alpha;
}

void set_ut(double ucon[NDIM], struct of_geom *geom)
{
  double AA, BB, CC;

  AA = geom->gcov[0][0];
  BB = 2.*(geom->gcov[0][1]*ucon[1] +
           geom->gcov[0][2]*ucon[2] +
           geom->gcov[0][3]*ucon[3]);
  CC = 1. + geom->gcov[1][1]*ucon[1]*ucon[1] +
            geom->gcov[2][2]*ucon[2]*ucon[2] +
            geom->gcov[3][3]*ucon[3]*ucon[3] +
       2. *(geom->gcov[1][2]*ucon[1]*ucon[2] +
            geom->gcov[1][3]*ucon[1]*ucon[3] +
            geom->gcov[2][3]*ucon[2]*ucon[3]);

  double discr = BB*BB - 4.*AA*CC;
  ucon[0] = (-BB - sqrt(discr))/(2.*AA);
}

void get_prim_bondi(int i, int j, int k, double P[NVAR])
{
  double r, th, X[NDIM];
  coord(i, j, k, CENT, X);
  bl_coord(X, &r, &th);

  while (r < Reh) {
    i++;
    coord(i, j, k, CENT, X);
    bl_coord(X, &r, &th);
  }

  double T = T_bondi[i][j];
  double ur = -C1/(pow(T,n)*pow(r,2));
  double rho = pow(T,n);
  double u = rho*T/(gam - 1.);
  double ucon_bl[NDIM], ucon_ks[NDIM], ucon_mks[NDIM];
  struct of_geom geom_bl, *geom_mks;

  blgset(i, j, &geom_bl);

  DLOOP1 {
    ucon_bl[mu] = 0.;
    ucon_ks[mu] = 0.;
    ucon_mks[mu] = 0.;
  }
  ucon_bl[1] = ur;

  set_ut(ucon_bl, &geom_bl);
  bl_to_ks(X, ucon_bl, ucon_ks);

  double dxdX[NDIM][NDIM], dXdx[NDIM][NDIM];
  set_dxdX(X, dxdX);
  invert(&dxdX[0][0], &dXdx[0][0]);
  DLOOP2 {
    ucon_mks[mu] += dXdx[mu][nu]*ucon_ks[nu];
  }

  geom_mks = get_geometry(i, j, k, CENT);
  fourvel_to_prim(ucon_mks, P, geom_mks);
  //printf("[%i %i %i] r = %e ucons = %e %e %e vr = %e\n",
  //  i, j, k, r, ucon_bl[1], ucon_ks[1], ucon_mks[1], P[U1]);

  P[RHO] = rho;
  P[UU] = u;
  P[B1] = 0.;
  P[B2] = 0.;
  P[B3] = 0.;

  /*if (j == N2/2){
    printf("\nu[1]: %e %e %e\n", ucon_bl[1], ucon_ks[1], ucon_mks[1]);
  printf("T rho ur u = %e %e %e %e\n", T, rho, ur, u);
  printf("P[] = %e %e %e %e\n", P[0], P[1], P[2], P[3]);
  }*/

  /*double gcov[NDIM][NDIM];
  bl_gcov_func(r, th, gcov);
  double ut = -sqrt((-1. - gcov[1][1]*ur*ur)/gcov[0][0]);
  double ucon[NDIM] = {ut, ur, 0, 0};

//  if ( j == NG)
//  printf("ucon[] = %e %e %e %e\n", ucon[0], ucon[1], ucon[2], ucon[3]);

  double trans[NDIM][NDIM];
  // transform to Kerr-Schild
  // Make transform matrix
  memset(trans, 0, 16*sizeof(double));
  for (int mu = 0; mu < NDIM; mu++) {
    trans[mu][mu] = 1.;
  }
  trans[0][1] = 2. * r / (r * r - 2. * r + a * a);
  trans[3][1] = a / (r * r - 2. * r + a * a);

  // Transform from BL to KS coordinates
  double tmp[NDIM];
  for (int mu = 0; mu < NDIM; mu++) tmp[mu] = 0.;
  for (int mu = 0; mu < NDIM; mu++) {
    for (int nu = 0; nu < NDIM; nu++) {
    tmp[mu] += trans[mu][nu]*ucon[nu];
    }
  }
  for (int mu = 0; mu < NDIM; mu++) ucon[mu] = tmp[mu];

  // Transform from KS to MKS coordinates
  set_dxdX(X, trans);
  for (int mu = 0; mu < NDIM; mu++) tmp[mu] = 0.;
  for (int mu = 0; mu < NDIM; mu++) {
    for (int nu = 0; nu < NDIM; nu++) {
    tmp[mu] += trans[mu][nu]*ucon[nu];
    }
  }
  for (int mu = 0; mu < NDIM; mu++) ucon[mu] = tmp[mu];

  struct of_geom *geom = get_geometry(i, j, k, CENT);
  double alpha = 1.0 / sqrt( -geom->gcon[0][0] ) ;
  double gamma = ucon[0] * alpha;

  double beta[NDIM];
  beta[1] = alpha * alpha * geom->gcon[0][1];
  beta[2] = alpha * alpha * geom->gcon[0][2];
  beta[3] = alpha * alpha * geom->gcon[0][3];

  P[RHO] = rho;
  P[UU] = u;
  P[U1] = ucon[1] + beta[1] * gamma / alpha;
  P[U2] = 0.;
  P[U3] = 0.;
  P[B1] = 0.;
  P[B2] = 0.;
  P[B3] = 0.;
  if (r < 5.) P[U1] = 0.;

  //P[U1] *= 0.;*/
}

/*void get_bondi_prim(double r, double P[NDIM])
{
  double T, rho, ur, ut;
  if (r > 2. + 1.e-8) {
    T = get_T(r);
    rho = pow(T, n);
    ur = -C1/r/r/pow(T,n);
    ut = (2.*ur/r+sqrt(ur*ur+1.-2./r))/(1.-2./r);
  } else {
    T = 1.;
    rho = RHOMIN;
    ur = 0.;
    ut = 1.;
  }

  P[RHO] = rho;
  P[UU] = n*rho*T;
  P[U1] = (ur + 2.*ut/(r + 2.))/r;
  P[U2] = 0.;
  P[U3] = 0.;
  P[B1] = 0.;
  P[B2] = 0.;
  P[B3] = 0.;
}*/

void init_prob()
{
  double r, th, X[NDIM];

  double mdot = 1.;
  double rs = 8.;
  n = 1./(gam - 1.);

  // Solution constants
  double uc = sqrt(mdot/(2.*rs));
  double Vc = -sqrt(pow(uc,2)/(1. - 3.*pow(uc,2)));
  double Tc = -n*pow(Vc,2)/((n + 1.)*(n*pow(Vc,2) - 1.));
  C1 = uc*pow(rs,2)*pow(Tc,n);
  C2 = pow(1. + (1. + n)*Tc,2)*(1. - 2.*mdot/rs + pow(C1,2)/
       (pow(rs,4)*pow(Tc,2*n)));

  printf("a = %e Reh = %e\n", a, Reh);

  printf("mdot = %e\n", mdot);
  printf("rs   = %e\n", rs);
  printf("n    = %e\n", n);
  printf("uc   = %e\n", uc);
  printf("Vc   = %e\n", Vc);
  printf("Tc   = %e\n", Tc);
  printf("C1   = %e\n", C1);
  printf("C2   = %e\n", C2);

  ZSLOOP(-NG,N1+NG-1,-NG,N2+NG-1,-NG,N3+NG-1) {
    coord(i, j, k, CENT, X);
    bl_coord(X, &r, &th);

    T_bondi[i][j] = get_T(r);
  }

  ZLOOP {
    get_prim_bondi(i, j, k, P[i][j][k]);
  } // ZSLOOP
}

// Convert Boyer-Lindquist four-velocity to MKS 3-velocity
void coord_transform(double *Pr, int ii, int jj)
{
  double X[NDIM], r, th, ucon[NDIM], trans[NDIM][NDIM], tmp[NDIM];
  double AA, BB, CC, discr;
  double alpha, gamma, beta[NDIM];
  struct of_geom *geom, blgeom;

  coord(ii, jj, 0, CENT, X);
  bl_coord(X, &r, &th);
  blgset(ii, jj, &blgeom);

  ucon[1] = Pr[U1];
  ucon[2] = Pr[U2];
  ucon[3] = Pr[U3];

  AA = blgeom.gcov[0][0];
  BB = 2. * (blgeom.gcov[0][1] * ucon[1] +
       blgeom.gcov[0][2] * ucon[2] +
       blgeom.gcov[0][3] * ucon[3]);
  CC = 1. +
      blgeom.gcov[1][1] * ucon[1] * ucon[1] +
      blgeom.gcov[2][2] * ucon[2] * ucon[2] +
      blgeom.gcov[3][3] * ucon[3] * ucon[3] +
      2. * (blgeom.gcov[1][2] * ucon[1] * ucon[2] +
      blgeom.gcov[1][3] * ucon[1] * ucon[3] +
      blgeom.gcov[2][3] * ucon[2] * ucon[3]);

  discr = BB * BB - 4. * AA * CC;
  ucon[0] = (-BB - sqrt(discr)) / (2. * AA);
  // This is ucon in BL coords

  // transform to Kerr-Schild
  // Make transform matrix
  memset(trans, 0, 16*sizeof(double));
  for (int mu = 0; mu < NDIM; mu++) {
    trans[mu][mu] = 1.;
  }
  trans[0][1] = 2. * r / (r * r - 2. * r + a * a);
  trans[3][1] = a / (r * r - 2. * r + a * a);

  // Transform from BL to KS coordinates
  for (int mu = 0; mu < NDIM; mu++) tmp[mu] = 0.;
  for (int mu = 0; mu < NDIM; mu++) {
    for (int nu = 0; nu < NDIM; nu++) {
     tmp[mu] += trans[mu][nu]*ucon[nu];
    }
  }
  for (int mu = 0; mu < NDIM; mu++) ucon[mu] = tmp[mu];

  // Transform from KS to MKS coordinates
  set_dxdX(X, trans);
  for (int mu = 0; mu < NDIM; mu++) tmp[mu] = 0.;
  for (int mu = 0; mu < NDIM; mu++) {
    for (int nu = 0; nu < NDIM; nu++) {
     tmp[mu] += trans[mu][nu]*ucon[nu];
    }
  }
  for (int mu = 0; mu < NDIM; mu++) ucon[mu] = tmp[mu];

  //ucon[1] *= (1. / (r - R0));
  //ucon[1] /= dr_dx(X[1]);
  //ucon[2] *=
  //    (1. / (M_PI + (1. - hslope) * M_PI * cos(2. * M_PI * X[2])));
  //ucon[2] /= dth_dx(X[2]);

  // Solve for v. Use same u^t, unchanged under KS -> KS'
  geom = get_geometry(ii, jj, 0, CENT) ;
  alpha = 1.0 / sqrt( -geom->gcon[0][0] ) ;
  gamma = ucon[0] * alpha;

  beta[1] = alpha * alpha * geom->gcon[0][1];
  beta[2] = alpha * alpha * geom->gcon[0][2];
  beta[3] = alpha * alpha * geom->gcon[0][3];

  Pr[U1] = ucon[1] + beta[1] * gamma / alpha;
  Pr[U2] = ucon[2] + beta[2] * gamma / alpha;
  Pr[U3] = ucon[3] + beta[3] * gamma / alpha;
}

double lfish_calc(double r)
{
  return (((pow(a, 2) - 2. * a * sqrt(r) + pow(r, 2)) *
     ((-2. * a * r *
       (pow(a, 2) - 2. * a * sqrt(r) +
        pow(r,
      2))) / sqrt(2. * a * sqrt(r) + (-3. + r) * r) +
      ((a + (-2. + r) * sqrt(r)) * (pow(r, 3) + pow(a, 2) *
      (2. + r))) / sqrt(1 + (2. * a) / pow (r, 1.5) - 3. / r)))
    / (pow(r, 3) * sqrt(2. * a * sqrt(r) + (-3. + r) * r) *
       (pow(a, 2) + (-2. + r) * r))
      );
}

void bound_gas_prob_x1r(int i, int j, int k, grid_prim_type P)
{
  get_prim_bondi(i, j, k, P[i][j][k]);
}