text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <time.h>
#include <gsl/gsl_math.h>

#include "cosmocalc.h"

void test_nonlinear_corrfunc(void)
{
  double lnrmin = log(1e-5);
  double lnrmax = log(3e2);
  int N = 10000;
  double r,dlnr = (lnrmax - lnrmin)/N;
  int i;
  FILE *fp;
  double a;

  a = 1.0;
  fp = fopen("xinltest.txt","w");
  for(i=0;i<N;++i)
    {
      r = exp(dlnr*i + lnrmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",r,nonlinear_corrfunc(r,a),linear_corrfunc(r,a));
    }
  fclose(fp);
}

void test_nonlinear_powspec(void)
{
  double lnkmin = log(1e-4);
  double lnkmax = log(1e2);
  double dlnk,k;
  int N = 10000,i;
  FILE *fp;
  double a = 1.0;

  dlnk = (lnkmax - lnkmin)/N;

  a = 1.0;
  fp = fopen("pknltest.txt","w");
  for(i=0;i<N;++i)
    {
      k = exp(dlnk*i + lnkmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",k,nonlinear_powspec(k,a),linear_powspec(k,a));
    }
  fclose(fp);

  a = 1.0/(1.0 + 0.5);
  fp = fopen("pknltest1.txt","w");
  for(i=0;i<N;++i)
    {
      k = exp(dlnk*i + lnkmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",k,nonlinear_powspec(k,a),linear_powspec(k,a));
    }
  fclose(fp);

  a = 1.0/(1.0 + 1.0);
  fp = fopen("pknltest2.txt","w");
  for(i=0;i<N;++i)
    {
      k = exp(dlnk*i + lnkmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",k,nonlinear_powspec(k,a),linear_powspec(k,a));
    }
  fclose(fp);

  a = 1.0/(1.0 + 2.00);
  fp = fopen("pknltest3.txt","w");
  for(i=0;i<N;++i)
    {
      k = exp(dlnk*i + lnkmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",k,nonlinear_powspec(k,a),linear_powspec(k,a));
    }
  fclose(fp);

  a = 1.0/(1.0 + 3.00);
  fp = fopen("pknltest4.txt","w");
  for(i=0;i<N;++i)
    {
      k = exp(dlnk*i + lnkmin);
      fprintf(fp,"%.20e\t%.20e\t%.20e\n",k,nonlinear_powspec(k,a),linear_powspec(k,a));
    }
  fclose(fp);
}

void test_linxi_fftlog(void)
{
  double rmin;
  long i,N = 2048;
  double dlnr;
  FILE *fp;
  double r;
  double a = 1.0;

  double *pk,*xi;
  double k0,L,q,r0,k;
  double kmin = 1e-7;
  double kmax = 1e7;
  double dlnk = log(kmax/kmin)/N;

  L = log(kmax/kmin);
  k0 = exp(log(kmin) + L/2.0);
  q = 0.0;

  pk = (double*)malloc(sizeof(double)*N);
  assert(pk != NULL);
  xi = (double*)malloc(sizeof(double)*N);
  assert(xi != NULL);

  for(i=0;i<N;++i)
    {
      k = exp(i*dlnk + log(kmin));
      pk[i] = linear_powspec_exact(k,a);
    }

  r0 = 1.0;
  compute_discrete_spherical_fft(pk,N,k0,L,q,xi,&r0);

  if(0)
    {
      for(i=0;i<N;++i)
	xi[i] = xi[i]/pow(2.0*M_PI,3.0);
      compute_discrete_spherical_fft(xi,N,r0,L,q,pk,&k0);

      kmin = exp(log(k0)-L/2.0);
      //fprintf(stderr,"k0 = %e, kmin = %e, L = %f, N = %ld\n",k0,kmin,L,N);
      fp = fopen("linxitest.txt","w");
      for(i=0;i<N;++i)
	{
	  k = exp(i*dlnk + log(kmin));
	  //fprintf(stderr,"k = %e, pk in,out = %e|%e\n",k,linear_powspec(k,a),pk[i]);
	  fprintf(fp,"%.20e \t %.20e \t %.20e \n",k,pk[i],linear_powspec_exact(k,a));
	}
      fclose(fp);

      return;
    }

  rmin = exp(log(r0)-L/2.0);
  dlnr = dlnk;
  fp = fopen("linxitest.txt","w");
  for(i=0;i<N;++i)
    {
      r = exp(i*dlnr + log(rmin));
      xi[i] = xi[i]/pow(2.0*M_PI,3.0);
      //fprintf(stderr,"%.20e \t %.20e \t %.20e \n",r,xi[i],linear_corrfunc_exact(r,a));
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",r,xi[i],linear_corrfunc_exact(r,a));
    }
  fclose(fp);
}

void test_linxi(void)
{
  double rmin=1e-9;
  double rmax=1e5;
  long i,Ntest = 1000;
  double dlnr = log(rmax/rmin)/Ntest;
  FILE *fp;
  double r;
  double a = 1.0;

  fp = fopen("linxitest.txt","w");
  for(i=0;i<Ntest;++i)
    {
      r = exp(i*dlnr + log(rmin));
      //fprintf(stderr,"%.20e \t %.20e \t %.20e \n",r,linear_corrfunc_exact(r,a),linear_corrfunc_exact(r,a));
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",r,linear_corrfunc(r,a),linear_corrfunc_exact(r,a));
    }
  fclose(fp);
}

void test_biasfunc(void)
{
  FILE *fp;
  double m,lgm;
  double mmin = 6.0;
  double mmax = 16.0;
  double dm = 0.01;
  double a = 1.0;

  fp = fopen("btest.txt","w");
  lgm = mmax;
  while(lgm >= mmin)
    {
      m = pow(10.0,lgm);
      fprintf(fp,"%.20e \t %.20e \n",1.686/sigmaMtophat(m,a),bias_function(m,a));
      lgm -= dm;
    }
  fclose(fp);

}

void test_massfunc(void)
{
  FILE *fp;
  double m,lgm;
  double mmin = 10.0;
  double mmax = 16.0;
  double dm = 0.01;
  double a = 1.0;

  fp = fopen("mftest.txt","w");
  lgm = mmax;
  while(lgm >= mmin)
    {
      m = pow(10.0,lgm);
      fprintf(fp,"%.20e \t %.20e \n",m,mass_function(m,a)*m*m/RHO_CRIT/cosmoData.OmegaM);
      lgm -= dm;
    }
  fclose(fp);
}

void test_peakheight(void)
{
  FILE *fp;
  double r;
  double rmin = 1e-40;
  double rmax = 1e3;
  double dr = 0.99;
  double a = 1.0;
  double numin = 0.15;
  double numax = 6.0;
  double dnu = 0.99;
  double nu;

  double kmin=1e-10;
  double kmax=1e10;
  long i,Ntest = 1000;
  double dlnk = log(kmax/kmin)/Ntest;
  double k,p = 1e-6;

  fp = fopen("phtest_integ.txt","w");
  for(i=0;i<Ntest;++i)
    {
      k = (i*dlnk + log(kmin));
      fprintf(fp,"%.20e \t %.20e \n",exp(k),tophatradnorm_linear_powspec_exact_nonorm_lnk_integ_funct_I0(k,&p));
    }
  fclose(fp);

  fp = fopen("phtest.txt","w");
  r = rmax;
  while(r >= rmin)
    {
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",r,sigmaRtophat(r,a),sigmaRtophat_exact(r,a));
      r *= dr;
    }
  fclose(fp);

  fp = fopen("phtestrev.txt","w");
  nu = numax;
  while(nu >= numin)
    {
      fprintf(fp,"%.20e \t %.20e \n",inverse_nuRtophat(nu,a),DELTAC/nu);
      nu *= dnu;
    }
  fclose(fp);
}

void test_gf(void)
{
  FILE *fp;
  double a;
  double amin = 0.1;
  double amax = 1.0;
  double da = 0.0001;

  fp = fopen("gftest.txt","w");
  a = amax;
  while(a >= amin)
    {
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",a,growth_function(a),growth_function_exact(a));
      a -= da;
    }
  fclose(fp);
}

void test_linpk(void)
{
  double kmin=1e-8;
  double kmax=1e40;
  long i,Ntest = 10000;
  double dlnk = log(kmax/kmin)/Ntest;
  FILE *fp;
  double k;
  double a = 1.0;

  fp = fopen("pktest.txt","w");
  for(i=0;i<Ntest;++i)
    {
      k = exp(i*dlnk + log(kmin));
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",k,linear_powspec(k,a),linear_powspec_exact(k,a));
    }
  fclose(fp);

}

void test_transfunct(void)
{
  double kmin=1e-8;
  double kmax=1e30;
  long i,Ntest = 10000;
  double dlnk = log(kmax/kmin)/Ntest;
  FILE *fp;
  double k;

  fp = fopen("tftest.txt","w");
  for(i=0;i<Ntest;++i)
    {
      k = exp(i*dlnk + log(kmin));
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",k,transfer_function(k),transfunct_eh98(k));
    }
  fclose(fp);

  cosmoData.OmegaM = 0.2;
  cosmoData.h = 0.5;
  cosmoData.OmegaB = cosmoData.OmegaM*0.5;
  cosmoData.cosmoNum = 2;

  fp = fopen("tftest_highb.txt","w");
  for(i=0;i<Ntest;++i)
    {
      k = exp(i*dlnk + log(kmin));
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",k,transfer_function(k),transfunct_eh98(k));
    }
  fclose(fp);
}

void test_distances(void)
{
  FILE *fp;

  double amin = 0.1;
  double amax = 1.0;
  double da = 0.0001;
  double dmin = 0.0;
  double dmax = 3000.0;
  double dd = 1.0;
  double a,d;

  fp = fopen("disttest_a2d.txt","w");
  a = amax;
  while(a >= amin)
    {
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",a,comvdist(a),comvdist_exact(a));
      a -= da;
    }
  fclose(fp);

  fp = fopen("disttest_d2a.txt","w");
  d = dmin;
  while(d <= dmax)
    {
      a = acomvdist(d);
      fprintf(fp,"%.20e \t %.20e \t %.20e \n",d,a,comvdist_exact(a));
      d += dd;
    }
  fclose(fp);
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <gsl\gsl>

namespace Library
{
	class DepthStencilStates final
	{
	public:
		inline static winrt::com_ptr<ID3D11DepthStencilState> DefaultDepthCulling;
		inline static winrt::com_ptr<ID3D11DepthStencilState> NoDepthCulling;

		static void Initialize(gsl::not_null<ID3D11Device*> direct3DDevice);
		static void Shutdown();

		DepthStencilStates() = delete;
		DepthStencilStates(const DepthStencilStates&) = delete;
		DepthStencilStates& operator=(const DepthStencilStates&) = delete;
		DepthStencilStates(DepthStencilStates&&) = delete;
		DepthStencilStates& operator=(DepthStencilStates&&) = delete;
		~DepthStencilStates() = default;
	};
}

ext: h


text:
/* permutation/gsl_permutation.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_PERMUTATION_H__
#define __GSL_PERMUTATION_H__

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

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

__BEGIN_DECLS

struct gsl_permutation_struct
{
  size_t size;
  size_t *data;
};

typedef struct gsl_permutation_struct gsl_permutation;

GSL_EXPORT gsl_permutation *gsl_permutation_alloc (const size_t n);
GSL_EXPORT gsl_permutation *gsl_permutation_calloc (const size_t n);
GSL_EXPORT void gsl_permutation_init (gsl_permutation * p);
GSL_EXPORT void gsl_permutation_free (gsl_permutation * p);
GSL_EXPORT int gsl_permutation_memcpy (gsl_permutation * dest, const gsl_permutation * src);

GSL_EXPORT int gsl_permutation_fread (FILE * stream, gsl_permutation * p);
GSL_EXPORT int gsl_permutation_fwrite (FILE * stream, const gsl_permutation * p);
GSL_EXPORT int gsl_permutation_fscanf (FILE * stream, gsl_permutation * p);
GSL_EXPORT int gsl_permutation_fprintf (FILE * stream, const gsl_permutation * p, const char *format);

GSL_EXPORT size_t gsl_permutation_size (const gsl_permutation * p);
GSL_EXPORT size_t * gsl_permutation_data (const gsl_permutation * p);

GSL_EXPORT size_t gsl_permutation_get (const gsl_permutation * p, const size_t i);
GSL_EXPORT int gsl_permutation_swap (gsl_permutation * p, const size_t i, const size_t j);

GSL_EXPORT int gsl_permutation_valid (gsl_permutation * p);
GSL_EXPORT void gsl_permutation_reverse (gsl_permutation * p);
GSL_EXPORT int gsl_permutation_inverse (gsl_permutation * inv, const gsl_permutation * p);
GSL_EXPORT int gsl_permutation_next (gsl_permutation * p);
GSL_EXPORT int gsl_permutation_prev (gsl_permutation * p);
GSL_EXPORT int gsl_permutation_mul (gsl_permutation * p, const gsl_permutation * pa, const gsl_permutation * pb);

GSL_EXPORT int gsl_permutation_linear_to_canonical (gsl_permutation * q, const gsl_permutation * p);
GSL_EXPORT int gsl_permutation_canonical_to_linear (gsl_permutation * p, const gsl_permutation * q);

GSL_EXPORT size_t gsl_permutation_inversions (const gsl_permutation * p);
GSL_EXPORT size_t gsl_permutation_linear_cycles (const gsl_permutation * p);
GSL_EXPORT size_t gsl_permutation_canonical_cycles (const gsl_permutation * q);

#ifdef HAVE_INLINE

extern inline
size_t
gsl_permutation_get (const gsl_permutation * p, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= p->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return p->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_PERMUTATION_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/* randist/fdist.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The F distribution has the form

   p(x) dx = (nu1^(nu1/2) nu2^(nu2/2) Gamma((nu1 + nu2)/2) /
   Gamma(nu1/2) Gamma(nu2/2)) *
   x^(nu1/2 - 1) (nu2 + nu1 * x)^(-nu1/2 -nu2/2) dx

   The method used here is the one described in Knuth */

double
gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2)
{

  double Y1 =  gsl_ran_gamma (r, nu1 / 2, 2.0);
  double Y2 =  gsl_ran_gamma (r, nu2 / 2, 2.0);

  double f = (Y1 * nu2) / (Y2 * nu1);

  return f;
}

double
gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2)
{
  if (x < 0)
    {
      return 0 ;
    }
  else
    {
      double p;
      double lglg = (nu1 / 2) * log (nu1) + (nu2 / 2) * log (nu2) ;

      double lg12 = gsl_sf_lngamma ((nu1 + nu2) / 2);
      double lg1 = gsl_sf_lngamma (nu1 / 2);
      double lg2 = gsl_sf_lngamma (nu2 / 2);

      p = exp (lglg + lg12 - lg1 - lg2)
	* pow (x, nu1 / 2 - 1) * pow (nu2 + nu1 * x, -nu1 / 2 - nu2 / 2);

      return p;
    }
}

ext: c


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

//
// Created by cleve on 2/1/2022.
//
#pragma once

#include "types.h"

#include <utility>
#include <cstdint>
#include <string_view>
#include <bit>

#include <gsl/gsl>

namespace jpeg_lite::decoder::transform
{
[[nodiscard]] FORCE_INLINE gsl::index coordinate_to_zigzag(gsl::index x, gsl::index y);

[[nodiscard]] FORCE_INLINE std::pair<gsl::index, gsl::index> zigzag_to_coordinate(gsl::index zigzag);

[[nodiscard, maybe_unused]] int16_t binary_string_to_int16(std::string_view sv);

[[nodiscard]] static inline FORCE_INLINE constexpr int16_t value_category(int16_t val)
{
	if (val == 0x0000)
		return 0;
	return std::bit_width(static_cast<uint16_t>(val > 0 ? val : -val));
}
}

max_stars_repo_name: SmartPolarBear/jpeg-lite

max_issues_repo_name: SmartPolarBear/jpeg-lite

max_forks_repo_name: SmartPolarBear/jpeg-lite


text:
//============================================================================
// Name        : main.c
// Author      : Cesare De Cal
// Version     :
// Copyright   : Cesare De Cal
// Description : Exercise 7 - System of Linear Equations
//============================================================================

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

// Useful: https://gist.github.com/bjd2385/7f4685e703f7437e513608f41c65bbd7 (LU decomposition)
// Useful: https://github.com/rwinlib/gsl/blob/master/include/gsl/gsl_linalg.h#L173 (GEPP)

gsl_matrix *createMatrix(int size) {
  gsl_matrix *matrix = gsl_matrix_alloc(size, size);
  gsl_matrix_set_zero(matrix);

  for (int i = 1; i <= size; i++) {
    for (int j = 1; j <= size; j++) {
      // Fill diagonal
      if (i == j) {
        if ((i+1) % 3 == 0) {
          double value = (2.0*(i+1))/3;
          gsl_matrix_set(matrix, i-1, j-1, value);
        } else {
          gsl_matrix_set(matrix, i-1, j-1, 1);
        }
        // Fill adjacent lower diagonal
      } else if (i == j+1) {
        gsl_matrix_set(matrix, i-1, j-1, 1);
        // Fill adjacent upper diagonal
      } else if (i+1 == j) {
        gsl_matrix_set(matrix, i-1, j-1, -1);
      }
    }
  }
  return matrix;
}

// Helper methods

void printVectorContents(gsl_vector *vector) {
  for (int i = 0; i < vector->size; i++) {
    printf("%.9e\\\\\n", gsl_vector_get(vector, i));
  }
  printf("\n");
}

void printMatrixContents(gsl_matrix *matrix) {
  for (int i = 0; i < matrix->size1; ++i) {
    for (int j = 0; j < matrix->size2; ++j) {
      printf("%.9e ", gsl_matrix_get(matrix, i, j));

      if (j == matrix->size2-1) {
        printf("\\\\");
      } else {
        printf("& ");
      }
    }
    printf("\n");
  }
  printf("\n");
}

void computeForOrder(int size, int shouldPrintVerboseOutput) {
  gsl_matrix *matrixA = createMatrix(size);
  if (shouldPrintVerboseOutput) {
    printf("Creating matrix A:\n");
    printMatrixContents(matrixA);
  }

  // Right-hand side vector y
  gsl_vector *bVector = gsl_vector_alloc(matrixA->size2);
  gsl_vector_set_zero(bVector);
  // Fill element at index 0 with 1
  gsl_vector_set(bVector, 0, 1);

  if (shouldPrintVerboseOutput) {
    printf("Creating vector y (b):\n");
    printVectorContents(bVector);
  }

  // Solve system using LU decomposition
  gsl_matrix *luMatrix = gsl_matrix_alloc(matrixA->size1, matrixA->size2);
  gsl_permutation *permutation = gsl_permutation_alloc(matrixA->size1);
  int signum; // Sign of the permutation
  // Copy A over newly created matrixA
  gsl_matrix_memcpy(luMatrix, matrixA);
  // Compute LU decomposition
  gsl_linalg_LU_decomp(luMatrix, permutation, &signum);

  if (shouldPrintVerboseOutput) {
    printf("LU decomposition of A:\n");
    printMatrixContents(luMatrix);
  }

  // Result vector v
  gsl_vector *xVector = gsl_vector_alloc(matrixA->size2);
  gsl_vector_set_zero(xVector);

  gsl_linalg_LU_solve(luMatrix, permutation, bVector, xVector);

  if (shouldPrintVerboseOutput) {
    printf("Solutions x vector:\n");
    printVectorContents(xVector);
  }

  // Calculate condition number
  // There is no function in GSL to directly compute this number
  // the condition number is given by taking the absolute value of the ratio of the largest singular value and the smallest singular value
  // cond(A) = abs( largest_sing_val / smallest_sing_val )
  gsl_matrix *matrixV = gsl_matrix_alloc(matrixA->size1, matrixA->size2);
  gsl_vector *vectorS = gsl_vector_alloc(matrixA->size1);
  gsl_vector *vectorWorkspace = gsl_vector_alloc(matrixA->size1);
  gsl_linalg_SV_decomp(matrixA, matrixV, vectorS, vectorWorkspace);

  if (shouldPrintVerboseOutput) {
    printf("Singular diagonal vector:\n");
    printVectorContents(vectorS);
  }

  double minSingularValue, maxSingularValue;
  gsl_vector_minmax(vectorS, &minSingularValue, &maxSingularValue);
  if (shouldPrintVerboseOutput) {
    printf("Max: %.9e, min: %.9e\n", maxSingularValue, minSingularValue);
  }

  double conditionNumber = fabs(maxSingularValue/minSingularValue);
  double x1 = gsl_vector_get(xVector, 0);
  double error = (M_E - 2) - x1;

  printf("%d & %.9e & %.9e & %.9e \\\\", size, x1, error, conditionNumber);

  // Free memory space
  gsl_matrix_free(matrixA);
  gsl_vector_free(bVector);
  gsl_matrix_free(luMatrix);
  gsl_permutation_free(permutation);
  gsl_vector_free(xVector);
  gsl_matrix_free(matrixV);
  gsl_vector_free(vectorS);
  gsl_vector_free(vectorWorkspace);
}

int main() {
  printf("This program is going to generate a LaTeX friendly output to be used in the report table.\n");
  printf("------------------------------\n");
  printf("n, x1, error, condition number\n");
  for (int n = 1; n <= 50; n++) {
    computeForOrder(n, 0);
    printf("\n");
  }
  printf("------------------------------\n");
  return 0;
}

max_stars_repo_name: csr/MATLAB-Scientific-Programming

max_issues_repo_name: csr/MATLAB-Scientific-Programming

max_forks_repo_name: csr/MATLAB-Scientific-Programming


text:
/**
* Copyright 2016 BitTorrent Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*    http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#pragma once

#include <scraps/config.h>

#include <scraps/Byte.h>

#if SCRAPS_APPLE
    #include <scraps/apple/SHA256.h>
    namespace scraps { using SHA256 = apple::SHA256; }
#else
    #include <scraps/sodium/SHA256.h>
    namespace scraps { using SHA256 = sodium::SHA256; }
#endif

#include <gsl.h>

namespace scraps {

template <typename BaseByteType>
struct SHA256ByteTag {};

template <typename BaseByteType>
using SHA256Byte = StrongByte<SHA256ByteTag<BaseByteType>>;

template <typename ByteT, std::ptrdiff_t N>
std::array<SHA256Byte<std::remove_const_t<ByteT>>, SHA256::kHashSize>
GetSHA256(gsl::span<ByteT, N> data) {
    std::array<SHA256Byte<std::remove_const_t<ByteT>>, SHA256::kHashSize> ret;

    SHA256 sha256;
    sha256.update(data.data(), data.size());
    sha256.finish(ret.data());

    return ret;
}

} // namespace scraps

max_stars_repo_name: carlbrown/scraps

max_issues_repo_name: carlbrown/scraps

max_forks_repo_name: carlbrown/scraps


text:
/*
 *  Copyright 2008-2009 NVIDIA Corporation
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

/*! \file defs.h
 *  \brief CBLAS utility definitions for interface routines
 */

#pragma once

#include <cblas.h>

namespace cusp
{
namespace system
{
namespace cpp
{
namespace detail
{
namespace cblas
{

struct cblas_format {};

struct upper   : public cblas_format {};
struct lower   : public cblas_format {};
struct unit    : public cblas_format {};
struct nonunit : public cblas_format {};

struct cblas_row_major { const static CBLAS_ORDER order = CblasRowMajor; };
struct cblas_col_major { const static CBLAS_ORDER order = CblasColMajor; };

template< typename LayoutFormat >
struct Orientation : thrust::detail::eval_if<
                        thrust::detail::is_same<LayoutFormat, cusp::row_major>::value,
                        thrust::detail::identity_<cblas_row_major>,
                        thrust::detail::identity_<cblas_col_major>
                     >
{};

} // end namespace cblas
} // end namespace detail
} // end namespace cpp
} // end namespace system
} // end namespace cusp

max_stars_repo_name: Raman-sh/cusplibrary


text:
/*

cp ../../../../native_ref/<target binary> libnetlib.so

gcc-mp-4.8 -O3 dgemmtest.c common.c -o dgemmtest -L. -lnetlib -I../../../../netlib/CBLAS
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-CBLAS.csv

gcc-mp-4.8 -O3 dgemmtest.c common.c -o dgemmtest -I/System/Library/Frameworks/vecLib.framework/Headers -framework veclib
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-veclib.csv

gcc-mp-4.8 -O3 dgemmtest.c common.c -o dgemmtest -I/opt/local/include /opt/local/lib/libatlas.a /opt/local/lib/libcblas.a /opt/local/lib/liblapack.a /opt/local/lib/libf77blas.a -lgfortran
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-atlas.csv

gcc-mp-4.8 -O3 dgemmtest.c common.c -o dgemmtest -I../../../../netlib/CBLAS -L/opt/intel/composerxe/mkl/lib -lmkl_rt
export DYLD_LIBRARY_PATH=/opt/intel/composerxe/mkl/lib:/opt/intel/composerxe/lib/
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-mkl.csv

gcc-mp-4.8 -O3 dgemmtest.c cudawrapper.c common.c -o dgemmtest -I../../../../netlib/CBLAS -I/usr/local/cuda/include/ -L/usr/local/cuda/lib -lcublas
export DYLD_LIBRARY_PATH=/usr/local/cuda/lib
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-cuda.csv

CLB=/Users/samuel/Documents/Projects/clBLAS/src
gcc-mp-4.8 -O3 dgemmtest.c clwrapper.c common.c -o dgemmtest -I../../../../netlib/CBLAS -I$CLB -I$CLB/include -L. -lclBLAS  -framework OpenCL
./dgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-clblas.csv

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <cblas.h>
#include "common.h"

// does AB == C ? If not, complain on stderr
void test(int m, double* a, double *b, double *c) {
	int i, j, k, exact = 0, wrong = 0;
	double diff;
	double* d = calloc(m * m, sizeof(double));
	for (i = 0 ; i < m ; i++) {
		for (j = 0 ; j < m ; j++) {
			for (k = 0 ; k < m ; k++) {
				d[i + j * m] += a[i + k * m] * b[j * m + k];
			}
		}
	}
	for (i = 0 ; i < m ; i++) {
		for (j = 0 ; j < m ; j++) {
			diff = c[i * m + j] - d[i * m + j];
			if (diff != 0.0) {
				exact++;
			}
			if (abs(diff) > 0.000001) {
				wrong++;
			}
		}
	}
	free(d);
	if (wrong > 0) {
		fprintf(stderr, "not exact = %d, wrong = %d\n", exact, wrong);
	}
}

long benchmark(int size) {
    int m = sqrt(size);
	long requestStart, requestEnd;

    double* a = random_array(m * m);
    double* b = random_array(m * m);
    double* c = calloc(m * m, sizeof(double));

	requestStart = currentTimeNanos();

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, m, m, 1, a, m, b, m, 0, c, m);

	requestEnd = currentTimeNanos();

#ifdef __TEST__
	test(m, a, b, c);
#endif

    free(a);
    free(b);
    free(c);

    return (requestEnd - requestStart);
  }

main()
{
	srand(time(NULL));

    double factor = 6.0 / 100.0;
    int i, j;
    for (i = 0 ; i < 10 ; i++) {
        for (j = 1 ; j <= 100 ; j++) {
            int size = (int) pow(10.0, factor * j);
            if (size < 10) continue;
            long took = benchmark(size);
            printf("\"%d\",\"%lu\"\n", size, took);
			fflush(stdout);
        }
    }
}


text:
/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_multimin.h>
#include <pthread.h>

/*User includes*/
#include "../Lib/FileUtils.h"
#include "../Lib/MatrixUtils.h"
#include "MetroLogStudent.h"

static char *usage[] = {"MetroLogStudent - Samples the Poisson Log-Student distn to SADs\n",
                        "Required parameters:\n",
			"   -out  filestub      output file stub\n",
			"   -in   filename      parameter file  \n",
                        "Optional:\n",
			"   -s    integer       generate integer mcmc samples\n",
			"   -seed   long        seed random number generator\n",
			"   -sigmaM float       std. dev. of mean prop. distn\n",
			"   -sigmaV float       \n",
			"   -sigmaN float       \n",
			"   -sigmaS float       \n",
			"   -v                  verbose\n"};

static int  nLines   = 11;

static int  verbose  = FALSE;

int main(int argc, char* argv[])
{
  int  i = 0, nNA     = 0;
  t_Params tParams;
  t_Data   tData;
  gsl_vector* ptX = gsl_vector_alloc(4); /*parameter estimates*/
  t_MetroInit atMetroInit[3];

  gsl_rng_env_setup();

  gsl_set_error_handler_off();

  /*get command line params*/
  getCommandLineParams(&tParams, argc, argv);

  /*read in abundance distribution*/
  readAbundanceData(tParams.szInputFile, &tData);

  /*set initial estimates for parameters*/
  gsl_vector_set(ptX, 0, INIT_M);
  gsl_vector_set(ptX, 1, INIT_V);
  gsl_vector_set(ptX, 2, INIT_N);
  gsl_vector_set(ptX, 3, tData.nL*2);

  printf("D = %d L = %d Chao = %f\n",tData.nL, tData.nJ, chao(&tData));

  minimiseSimplex(ptX, 4, (void*) &tData, &nLogLikelihood);

  outputResults(ptX, &tData);

  if(tParams.nIter > 0){
    mcmc(&tParams, &tData, ptX);
  }

  /*free up allocated memory*/
  gsl_vector_free(ptX);

  freeAbundanceData(&tData);

  exit(EXIT_SUCCESS);
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

  /*get out file stub*/
  ptParams->szOutFileStub  = extractParameter(argc,argv,OUT_FILE_STUB,ALWAYS);
  if(ptParams->szOutFileStub == NULL)
    goto error;

  /*get out file stub*/
  szTemp  = extractParameter(argc,argv,SEED,OPTION);
  if(szTemp != NULL){
    ptParams->lSeed = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->lSeed = 0;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  szTemp  = extractParameter(argc,argv,SIGMA_M,OPTION);
  if(szTemp != NULL){
    ptParams->dSigmaM = strtod(szTemp,&cError);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->dSigmaM = DEF_SIGMA;
  }

  szTemp  = extractParameter(argc,argv,SIGMA_V,OPTION);
  if(szTemp != NULL){
    ptParams->dSigmaV = strtod(szTemp,&cError);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->dSigmaV = DEF_SIGMA;
  }

  szTemp  = extractParameter(argc,argv,SIGMA_N,OPTION);
  if(szTemp != NULL){
    ptParams->dSigmaN = strtod(szTemp,&cError);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->dSigmaN = DEF_SIGMA;
  }

  szTemp  = extractParameter(argc,argv,SIGMA_S,OPTION);
  if(szTemp != NULL){
    ptParams->dSigmaS = strtod(szTemp,&cError);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->dSigmaS = DEF_SIGMA_S;
  }

  szTemp  = extractParameter(argc,argv,SAMPLE,OPTION);
  if(szTemp != NULL){
    ptParams->nIter = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nIter = 0;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readAbundanceData(const char *szFile, t_Data *ptData)
{
  int **aanAbund = NULL;
  int  i = 0, nNA = 0, nA = 0, nC = 0;
  int  nL = 0, nJ = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    szTok = strtok(szLine, DELIM);

    nNA = strtol(szTok,&pcError,10);
    if(*pcError != '\0'){
      goto formatError;
    }

    aanAbund = (int **) malloc(nNA*sizeof(int*));

    for(i = 0; i < nNA; i++){
      aanAbund[i] = (int *) malloc(sizeof(int)*2);

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM);

      nA = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      szTok = strtok(NULL, DELIM);

      nC = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      nL += nC;
      nJ += nC*nA;

      aanAbund[i][0]  = nA;
      aanAbund[i][1]  = nC;
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  ptData->nJ          = nJ;
  ptData->nL          = nL;
  ptData->aanAbund    = aanAbund;
  ptData->nNA         = nNA;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int compare_doubles(const void* a, const void* b)
{
  double* arg1 = (double *) a;
  double* arg2 = (double *) b;
  if( *arg1 < *arg2 ) return -1;
  else if( *arg1 == *arg2 ) return 0;
  else return 1;
}

double chao(t_Data *ptData)
{
  double n1 = 0.0, n2 = 0.0;
  int **aanAbund = ptData->aanAbund;

  if(aanAbund[0][0] == 1 && aanAbund[1][0] == 2){
    n1 = (double) aanAbund[0][1]; n2 = (double) aanAbund[1][1];

    return ((double) ptData->nL) + 0.5*((n1*n1)/n2);
  }
  else{
    return -1.0;
  }
}

double f1(double x, void *pvParams)
{
  t_LSParams *ptLSParams = (t_LSParams *) pvParams;
  double dMDash = ptLSParams->dMDash, dV = ptLSParams->dV, dNu = ptLSParams->dNu;
  int n = ptLSParams->n;
  double t = ((x - dMDash)*(x - dMDash))/dV;
  double dExp  = x*((double) n) - exp(x);
  double dF  = pow(1.0 + t/dNu, -0.5*(dNu + 1.0));

  return exp(dExp)*dF;
}

double f1Log(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV;
  int n = ptLNParams->n;
  double dTemp = (x - dMDash);
  double dExp  = x*((double) n) - exp(x) - 0.5*((dTemp*dTemp)/dV);
  double dRet  = exp(dExp);

  return dRet;
}

double derivExponent(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
  double dTemp = (x - dMDash)/dV, dRet = 0.0;

  dRet = ((double) n) - exp(x) - dTemp;

  return dRet;
}

double logStirlingsGamma(double dZ)
{
  return 0.5*log(2.0*M_PI) + (dZ - 0.5)*log(dZ) - dZ;
}

double logLikelihoodQuad(int n, double dMDash, double dV, double dNu)
{
  gsl_integration_workspace *ptGSLWS =
    gsl_integration_workspace_alloc(1000);
  double dLogFac1   = 0.0, dLogFacN  = 0.0;
  double dN = (double) n, dResult = 0.0, dError = 0.0, dPrecision = 0.0;
  gsl_function tGSLF;
  t_LSParams tLSParams;
  double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;

  tLSParams.n = n; tLSParams.dMDash = dMDash; tLSParams.dV = dV; tLSParams.dNu = dNu;

  tGSLF.function = &f1;
  tGSLF.params   = (void *) &tLSParams;

  if(dNu < MAX_MU_GAMMA){
    dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
  }
  else{
    dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
  }

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else if(n < 100){
    dLogFacN = gsl_sf_lngamma(dN + 1.0);
  }
  else{
    dLogFacN = logStirlingsGamma(dN + 1.0);
  }

  dA = -100.0; dB = 100.0;

  if(n < 10){
    dPrecision = HI_PRECISION;
  }
  else{
    dPrecision = LO_PRECISION;
  }

  gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);

  //printf("%f %f\n", dResult, dError);

  gsl_integration_workspace_free(ptGSLWS);

  return log(dResult) - dLogFacN + dLogFac1 - 0.5*log(dV);
}

int solveF(double x_lo, double x_hi, double (*f)(double, void*),
	   void* params, double tol, double *xsolve)
{
  int status, iter = 0, max_iter = 100;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double r = 0;
  gsl_function F;
  t_LNParams *ptLNParams = (t_LNParams *) params;
  F.function = f;
  F.params = params;

  //printf("%f %f %d %f %f\n",ptLNParams->dMDash, ptLNParams->dV, ptLNParams->n, x_lo, x_hi);
  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);

  do{
    iter++;
    status = gsl_root_fsolver_iterate (s);
    r = gsl_root_fsolver_root (s);
    x_lo = gsl_root_fsolver_x_lower (s);
    x_hi = gsl_root_fsolver_x_upper (s);

    status = gsl_root_test_interval (x_lo, x_hi, 0, tol);
  }
  while (status == GSL_CONTINUE && iter < max_iter);

  (*xsolve) = gsl_root_fsolver_root (s);
  gsl_root_fsolver_free (s);

  return status;
}

double logLikelihoodLNQuad(int n, double dMDash, double dV)
{
  gsl_integration_workspace *ptGSLWS =
    gsl_integration_workspace_alloc(1000);
  double dLogFac1   = 0.0, dLogFacN  = 0.0;
  double dResult = 0.0, dError = 0.0, dPrecision = 0.0;
  gsl_function tGSLF;
  t_LNParams tLNParams;
  double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;

  tLNParams.n = n; tLNParams.dMDash = dMDash; tLNParams.dV = dV;

  tGSLF.function = &f1Log;
  tGSLF.params   = (void *) &tLNParams;

  dLogFac1 = log(2.0*M_PI*dV);

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else{
    dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
  }

  if(dEst > dV){
    double dMax = 0.0;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    if(fabs(dUpper) > 1.0e-7){
      solveF(0.0, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);
    }

    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }
  else{
    double dMax = 0.0;
    double dLower = dEst - dV;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    if(fabs(dUpper - dLower) > 1.0e-7){
      solveF(dLower, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);
    }
    else{
      dMax = 0.5*(dLower + dUpper);
    }
    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }

  if(n < 10){
    dPrecision = HI_PRECISION;
  }
  else{
    dPrecision = LO_PRECISION;
  }

  gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);

  gsl_integration_workspace_free(ptGSLWS);

  return log(dResult) - dLogFacN -0.5*dLogFac1;
}

double logLikelihoodLNRampal(int n, double dMDash, double dV)
{
  double dN = (double) n;
  double dLogLik = 0.0, dTemp = gsl_pow_int(log(dN) - dMDash,2), dTemp3 = gsl_pow_int(log(dN) - dMDash,3);

  dLogLik = -0.5*log(2.0*M_PI*dV) - log(dN) - (dTemp/(2.0*dV));

  dLogLik += log(1.0 + 1.0/(2.0*dN*dV)*(dTemp/dV + log(dN) - dMDash - 1.0)
		 + 1.0/(6.0*dN*dN*dV*dV*dV)*(3.0*dV*dV - (3.0*dV - 2.0*dV*dV)*(dMDash - log(dN))
		 - 3.0*dV*dTemp + dTemp3));

  return dLogLik;
}

double logLikelihoodRampal(int n, double dMDash, double dV, double dNu)
{
  double dGamma = 0.5*(dNu + 1.0), dN = (double) n, dRN = 1.0/dN, dRSV = 1.0/(sqrt(dV)*sqrt(dNu));
  double dZ = (log(dN) - dMDash)*dRSV;
  double dDZDX = dRN*dRSV, dDZDX2 = -dRN*dRN*dRSV;
  double dF = (1.0 + dZ*dZ);
  double dA = 0.0, dB = 0.0, dTemp = 0.0;
  double dLogFac1 = 0.0;

  if(dNu < MAX_MU_GAMMA){
    dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
  }
  else{
    dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
  }

  dA = 4.0*dZ*dZ*dDZDX*dDZDX*dGamma*(dGamma + 1.0);
  dA /= dF*dF;

  dB = -2.0*dGamma*(dDZDX*dDZDX + dZ*dDZDX2);
  dB /= dF;

  dTemp = dRN + dA + dB;

  return -dGamma*log(dF) + log(dTemp) + dLogFac1 - 0.5*log(dV);
}

double nLogLikelihood(const gsl_vector * x, void * params)
{
  double dMDash  = gsl_vector_get(x,0), dV = gsl_vector_get(x,1);
  double dNu  = gsl_vector_get(x,2);
  int    nS = (int) floor(gsl_vector_get(x, 3));
  t_Data *ptData = (t_Data *) params;
  int    i       = 0;
  double dLogNot0 = 0.0, dLogL   = 0.0;
  double dLog0 = 0.0, dLog1 = 0.0, dLog2 = 0.0, dLog3 = 0.0;

  if(dV <= 0.0 || dNu < 1.0){
    return PENALTY;
  }

  for(i = 0; i < ptData->nNA; i++){
    double dLogP = 0.0;
    int    nA    = ptData->aanAbund[i][0];

    if(nA < MAX_QUAD){
      dLogP = logLikelihoodQuad(nA, dMDash, dV, dNu);
    }
    else{
      dLogP = logLikelihoodRampal(nA, dMDash, dV, dNu);
    }

    dLogL += ((double) ptData->aanAbund[i][1])*dLogP;

    dLogL -= gsl_sf_lnfact(ptData->aanAbund[i][1]);

  }

  dLog0 = logLikelihoodQuad(0, dMDash, dV, dNu);

  dLog1 = (nS - ptData->nL)*dLog0;

  dLog2 = - gsl_sf_lnfact(nS - ptData->nL);

  dLog3 = gsl_sf_lnfact(nS);

  dLogL += dLog1 + dLog2 + dLog3;

  /*return*/
  return -dLogL;
}

double negLogLikelihood(double dMDash, double dV, double dNu, int nS, void * params)
{
  t_Data *ptData = (t_Data *) params;
  int    i       = 0;
  double dLogNot0 = 0.0, dLogL   = 0.0;
  double dLog0 = 0.0, dLog1 = 0.0, dLog2 = 0.0, dLog3 = 0.0;

  if(dV <= 0.0 || dNu < 1.0){
    return PENALTY;
  }

  for(i = 0; i < ptData->nNA; i++){
    double dLogP = 0.0;
    int    nA    = ptData->aanAbund[i][0];

    if(nA < MAX_QUAD){
      dLogP = logLikelihoodQuad(nA, dMDash, dV, dNu);
    }
    else{
      dLogP = logLikelihoodRampal(nA, dMDash, dV, dNu);
    }

    dLogL += ((double) ptData->aanAbund[i][1])*dLogP;

    dLogL -= gsl_sf_lnfact(ptData->aanAbund[i][1]);

  }

  dLog0 = logLikelihoodQuad(0, dMDash, dV, dNu);

  dLog1 = (nS - ptData->nL)*dLog0;

  dLog2 = - gsl_sf_lnfact(nS - ptData->nL);

  dLog3 = gsl_sf_lnfact(nS);

  dLogL += dLog1 + dLog2 + dLog3;

  /*return*/
  return -dLogL;
}

int minimiseSimplex(gsl_vector* ptX, size_t nP, void* pvData, double (*f)(const gsl_vector*, void* params))
{
  const gsl_multimin_fminimizer_type *T =
    gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *ss;
  gsl_multimin_function minex_func;
  size_t iter = 0;
  int i = 0, status;
  double size;

  /* Initial vertex size vector */
  ss = gsl_vector_alloc (nP);

  /* Set all step sizes to default constant */
  for(i = 0; i < nP; i++){
    gsl_vector_set(ss, i,INIT_SIMPLEX_SIZE*fabs(gsl_vector_get(ptX,i)));
  }

  /* Initialize method and iterate */
  minex_func.f = f;
  minex_func.n = nP;
  minex_func.params = pvData;

  s = gsl_multimin_fminimizer_alloc (T, nP);
  gsl_multimin_fminimizer_set(s, &minex_func, ptX, ss);

  do{
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);

    if(status)
      break;

    size = gsl_multimin_fminimizer_size(s);
    status = gsl_multimin_test_size(size, MIN_SIMPLEX_SIZE);

    if(status == GSL_SUCCESS){
      for(i = 0; i < nP; i++){
	gsl_vector_set(ptX, i, gsl_vector_get(s->x, i));
      }

      if(verbose) printf("converged to minimum at\n");
    }

    if(verbose){
      printf ("%5d ", iter);

      for (i = 0; i < nP; i++) printf("%10.3e ", gsl_vector_get(s->x, i));

      printf("f() = %7.3f size = %.3f\n", s->fval, size);
    }
  }
  while(status == GSL_CONTINUE && iter < MAX_SIMPLEX_ITER);

  for(i = 0; i < nP; i++){
    gsl_vector_set(ptX, i, gsl_vector_get(s->x, i));
  }

  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  return status;
}

void freeAbundanceData(t_Data *ptData)
{
  int i = 0;

  for(i = 0; i < ptData->nNA; i++){
    free(ptData->aanAbund[i]);
  }
  free(ptData->aanAbund);
}

void getProposal(gsl_rng *ptGSLRNG, gsl_vector *ptXDash, gsl_vector *ptX, int* pnSDash, int nS, t_Params *ptParams)
{
  double dDeltaS =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaS);
  double dDeltaM =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaM);
  double dDeltaV =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaV);
  double dDeltaN =  gsl_ran_gaussian(ptGSLRNG, ptParams->dSigmaN);
  int    nSDash = 0;

  gsl_vector_set(ptXDash, 0, gsl_vector_get(ptX,0) + dDeltaM);
  gsl_vector_set(ptXDash, 1, gsl_vector_get(ptX,1) + dDeltaV);
  gsl_vector_set(ptXDash, 2, gsl_vector_get(ptX,2) + dDeltaN);

  //printf("%e %e %e\n",dDeltaA,dDeltaB,dDeltaG);

  nSDash = nS + (int) floor(dDeltaS);
  if(nSDash < 1){
    nSDash = 1;
  }
  (*pnSDash) = nSDash;
}

void outputResults(gsl_vector *ptX, t_Data *ptData)
{
  double dAlpha = 0.0, dBeta = 0.0, dGamma = 0.0, dS = 0.0, dL = 0.0;

  dAlpha = gsl_vector_get(ptX, 0);

  dBeta  = gsl_vector_get(ptX, 1);

  dGamma = gsl_vector_get(ptX, 2);

  dS = gsl_vector_get(ptX, 3);

  dL = nLogLikelihood(ptX, ptData);

  printf("\nML simplex: M = %.2f V = %.2f Nu = %.2f S = %.2f NLL = %.2f\n",dAlpha, dBeta, dGamma, dS, dL);
}

void* metropolis (void * pvInitMetro)
{
  t_MetroInit *ptMetroInit  = (t_MetroInit *) pvInitMetro;
  gsl_vector  *ptX          = ptMetroInit->ptX;
  t_Data      *ptData       = ptMetroInit->ptData;
  t_Params    *ptParams     = ptMetroInit->ptParams;
  gsl_vector  *ptXDash      = gsl_vector_alloc(4); /*proposal*/
  char *szSampleFile = (char *) malloc(MAX_LINE_LENGTH*sizeof(char));
  const gsl_rng_type *T;
  gsl_rng            *ptGSLRNG;
  FILE    *sfp = NULL;
  int nS = 0, nSDash = 0, nIter = 0;
  double dRand = 0.0, dNLL = 0.0;
  void   *pvRet = NULL;

  /*set up random number generator*/
  T        = gsl_rng_default;
  ptGSLRNG = gsl_rng_alloc (T);

  nS = (int) floor(gsl_vector_get(ptX,3));

  dNLL = negLogLikelihood(gsl_vector_get(ptX,0), gsl_vector_get(ptX,1), gsl_vector_get(ptX,2), nS,(void*) ptData);

  sprintf(szSampleFile,"%s_%d%s", ptParams->szOutFileStub, ptMetroInit->nThread, SAMPLE_FILE_SUFFIX);

  sfp = fopen(szSampleFile, "w");
  if(!sfp){
    exit(EXIT_FAILURE);
  }

  /*seed random number generator*/
  gsl_rng_set(ptGSLRNG, ptMetroInit->lSeed);

  /*now perform simple Metropolis algorithm*/
  while(nIter < ptParams->nIter){
    double dA = 0.0, dNLLDash = 0.0;

    getProposal(ptGSLRNG, ptXDash, ptX, &nSDash, nS, ptParams);

    dNLLDash = negLogLikelihood(gsl_vector_get(ptXDash,0), gsl_vector_get(ptXDash,1), gsl_vector_get(ptXDash,2), nSDash, (void*) ptData);
    //printf("X' %e %e %e %d %f\n", gsl_vector_get(ptXDash,0), gsl_vector_get(ptXDash,1), gsl_vector_get(ptXDash,2), nSDash, dNLLDash);
    //printf("X %e %e %e %d %f\n", gsl_vector_get(ptX,0), gsl_vector_get(ptX,1), gsl_vector_get(ptX,2), nS, dNLL);
    dA = exp(dNLL - dNLLDash);
    if(dA > 1.0){
      dA = 1.0;
    }

    dRand = gsl_rng_uniform(ptGSLRNG);

    if(dRand < dA){
      gsl_vector_memcpy(ptX, ptXDash);
      nS = nSDash;
      dNLL = dNLLDash;
      ptMetroInit->nAccepted++;
    }

    if(nIter % SLICE == 0){
      fprintf(sfp, "%d,%.10e,%.10e,%.10e,%d,%f\n",nIter,gsl_vector_get(ptX, 0), gsl_vector_get(ptX, 1), gsl_vector_get(ptX, 2), nS, dNLL);
      fflush(sfp);
    }

    nIter++;
  }

  fclose(sfp);

  /*free up allocated memory*/
  gsl_vector_free(ptXDash);
  free(szSampleFile);
  gsl_rng_free(ptGSLRNG);

  return pvRet;
}

void writeThread(t_MetroInit *ptMetroInit)
{
  gsl_vector *ptX = ptMetroInit->ptX;
    printf("%d: a = %.2f b = %.2f g = %.2f S = %.2f\n", ptMetroInit->nThread,
	   gsl_vector_get(ptX, 0),
	   gsl_vector_get(ptX, 1),
	   gsl_vector_get(ptX, 2),
	   gsl_vector_get(ptX, 3));
}

void mcmc(t_Params *ptParams, t_Data *ptData, gsl_vector* ptX)
{
  pthread_t thread1, thread2, thread3;
  int       iret1  , iret2  , iret3;
  gsl_vector *ptX1 = gsl_vector_alloc(4),
             *ptX2 = gsl_vector_alloc(4),
             *ptX3 = gsl_vector_alloc(4);
  t_MetroInit atMetroInit[3];

  printf("\nMCMC iter = %d sigmaM = %.2f sigmaV = %.2f sigmaN = %.2f sigmaS = %.2f\n",
	   ptParams->nIter, ptParams->dSigmaM, ptParams->dSigmaV, ptParams->dSigmaN, ptParams->dSigmaS);

  gsl_vector_memcpy(ptX1, ptX);

  gsl_vector_set(ptX2, 0, gsl_vector_get(ptX,0) + 2.0*ptParams->dSigmaM);
  gsl_vector_set(ptX2, 1, gsl_vector_get(ptX,1) + 2.0*ptParams->dSigmaV);
  gsl_vector_set(ptX2, 2, gsl_vector_get(ptX,2) + 2.0*ptParams->dSigmaN);
  gsl_vector_set(ptX2, 3, gsl_vector_get(ptX,3) + 2.0*ptParams->dSigmaS);

  gsl_vector_set(ptX3, 0, gsl_vector_get(ptX,0) - 2.0*ptParams->dSigmaM);
  gsl_vector_set(ptX3, 1, gsl_vector_get(ptX,1) - 2.0*ptParams->dSigmaV);
  gsl_vector_set(ptX3, 2, gsl_vector_get(ptX,2) - 2.0*ptParams->dSigmaN);

  if(gsl_vector_get(ptX,3) - 2.0*ptParams->dSigmaS > (double) ptData->nL){
	gsl_vector_set(ptX3, 3, gsl_vector_get(ptX,3) - 2.0*ptParams->dSigmaS);   }
  else{
	gsl_vector_set(ptX3, 3, (double) ptData->nL);   }

  atMetroInit[0].ptParams = ptParams;
  atMetroInit[0].ptData   = ptData;
  atMetroInit[0].ptX      = ptX1;
  atMetroInit[0].nThread  = 0;
  atMetroInit[0].lSeed    = ptParams->lSeed;
  atMetroInit[0].nAccepted = 0;

  atMetroInit[1].ptParams = ptParams;
  atMetroInit[1].ptData   = ptData;
  atMetroInit[1].ptX      = ptX2;
  atMetroInit[1].nThread  = 1;
  atMetroInit[1].lSeed    = ptParams->lSeed + 1;
  atMetroInit[1].nAccepted = 0;

  atMetroInit[2].ptParams = ptParams;
  atMetroInit[2].ptData   = ptData;
  atMetroInit[2].ptX      = ptX3;
  atMetroInit[2].nThread  = 2;
  atMetroInit[2].lSeed    = ptParams->lSeed + 2;
  atMetroInit[2].nAccepted = 0;

  writeThread(&atMetroInit[0]);
  writeThread(&atMetroInit[1]);
  writeThread(&atMetroInit[2]);

  iret1 = pthread_create(&thread1, NULL, metropolis, (void*) &atMetroInit[0]);
  iret2 = pthread_create(&thread2, NULL, metropolis, (void*) &atMetroInit[1]);
  iret3 = pthread_create(&thread3, NULL, metropolis, (void*) &atMetroInit[2]);
  pthread_join(thread1, NULL);
  pthread_join(thread2, NULL);
  pthread_join(thread3, NULL);

  printf("%d: accept. ratio %d/%d = %f\n", atMetroInit[0].nThread,
	 atMetroInit[0].nAccepted, ptParams->nIter,((double) atMetroInit[0].nAccepted)/((double) ptParams->nIter));

  printf("%d: accept. ratio %d/%d = %f\n", atMetroInit[1].nThread,
	 atMetroInit[1].nAccepted, ptParams->nIter,((double) atMetroInit[1].nAccepted)/((double) ptParams->nIter));

  printf("%d: accept. ratio %d/%d = %f\n", atMetroInit[2].nThread,
	 atMetroInit[2].nAccepted, ptParams->nIter, ((double) atMetroInit[2].nAccepted)/((double) ptParams->nIter));

  gsl_vector_free(ptX1); gsl_vector_free(ptX2); gsl_vector_free(ptX3);
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


text:
/* multifit_nlinear/gsl_multifit_nlinear.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 * Copyright (C) 2015, 2016 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MULTIFIT_NLINEAR_H__
#define __GSL_MULTIFIT_NLINEAR_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>

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

__BEGIN_DECLS

typedef enum
{
  GSL_MULTIFIT_NLINEAR_FWDIFF,
  GSL_MULTIFIT_NLINEAR_CTRDIFF
} gsl_multifit_nlinear_fdtype;

/* Definition of vector-valued functions and gradient with parameters
   based on gsl_vector */

typedef struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
  int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
               gsl_vector * fvv);
  size_t n;        /* number of functions */
  size_t p;        /* number of independent variables */
  void * params;   /* user parameters */
  size_t nevalf;   /* number of function evaluations */
  size_t nevaldf;  /* number of Jacobian evaluations */
  size_t nevalfvv; /* number of fvv evaluations */
} gsl_multifit_nlinear_fdf;

/* trust region subproblem method */
typedef struct
{
  const char *name;
  void * (*alloc) (const void * params, const size_t n, const size_t p);
  int (*init) (const void * vtrust_state, void * vstate);
  int (*preloop) (const void * vtrust_state, void * vstate);
  int (*step) (const void * vtrust_state, const double delta,
               gsl_vector * dx, void * vstate);
  int (*preduction) (const void * vtrust_state, const gsl_vector * dx,
                     double * pred, void * vstate);
  void (*free) (void * vstate);
} gsl_multifit_nlinear_trs;

/* scaling matrix specification */
typedef struct
{
  const char *name;
  int (*init) (const gsl_matrix * J, gsl_vector * diag);
  int (*update) (const gsl_matrix * J, gsl_vector * diag);
} gsl_multifit_nlinear_scale;

/*
 * linear least squares solvers - there are three steps to
 * solving a least squares problem using a trust region
 * method:
 *
 * 1. init: called once per iteration when a new Jacobian matrix
 *          is computed; perform factorization of Jacobian (qr,svd)
 *          or form normal equations matrix (cholesky)
 * 2. presolve: called each time a new LM parameter value mu is available;
 *              used for cholesky method in order to factor
 *              the (J^T J + mu D^T D) matrix
 * 3. solve: solve the least square system for a given rhs
 */
typedef struct
{
  const char *name;
  void * (*alloc) (const size_t n, const size_t p);
  int (*init) (const void * vtrust_state, void * vstate);
  int (*presolve) (const double mu, const void * vtrust_state, void * vstate);
  int (*solve) (const gsl_vector * f, gsl_vector * x,
                const void * vtrust_state, void * vstate);
  int (*rcond) (double * rcond, void * vstate);
  void (*free) (void * vstate);
} gsl_multifit_nlinear_solver;

/* tunable parameters */
typedef struct
{
  const gsl_multifit_nlinear_trs *trs;        /* trust region subproblem method */
  const gsl_multifit_nlinear_scale *scale;    /* scaling method */
  const gsl_multifit_nlinear_solver *solver;  /* solver method */
  gsl_multifit_nlinear_fdtype fdtype;         /* finite difference method */
  double factor_up;                           /* factor for increasing trust radius */
  double factor_down;                         /* factor for decreasing trust radius */
  double avmax;                               /* max allowed |a|/|v| */
  double h_df;                                /* step size for finite difference Jacobian */
  double h_fvv;                               /* step size for finite difference fvv */
} gsl_multifit_nlinear_parameters;

typedef struct
{
  const char *name;
  void * (*alloc) (const gsl_multifit_nlinear_parameters * params,
                   const size_t n, const size_t p);
  int (*init) (void * state, const gsl_vector * wts,
               gsl_multifit_nlinear_fdf * fdf, const gsl_vector * x,
               gsl_vector * f, gsl_matrix * J, gsl_vector * g);
  int (*iterate) (void * state, const gsl_vector * wts,
                  gsl_multifit_nlinear_fdf * fdf, gsl_vector * x,
                  gsl_vector * f, gsl_matrix * J, gsl_vector * g,
                  gsl_vector * dx);
  int (*rcond) (double * rcond, void * state);
  double (*avratio) (void * state);
  void (*free) (void * state);
} gsl_multifit_nlinear_type;

/* current state passed to low-level trust region algorithms */
typedef struct
{
  const gsl_vector * x;             /* parameter values x */
  const gsl_vector * f;             /* residual vector f(x) */
  const gsl_vector * g;             /* gradient J^T f */
  const gsl_matrix * J;             /* Jacobian J(x) */
  const gsl_vector * diag;          /* scaling matrix D */
  const gsl_vector * sqrt_wts;      /* sqrt(diag(W)) or NULL for unweighted */
  const double *mu;                 /* LM parameter */
  const gsl_multifit_nlinear_parameters * params;
  void *solver_state;               /* workspace for linear least squares solver */
  gsl_multifit_nlinear_fdf * fdf;
  double *avratio;                  /* |a| / |v| */
} gsl_multifit_nlinear_trust_state;

typedef struct
{
  const gsl_multifit_nlinear_type * type;
  gsl_multifit_nlinear_fdf * fdf ;
  gsl_vector * x;             /* parameter values x */
  gsl_vector * f;             /* residual vector f(x) */
  gsl_vector * dx;            /* step dx */
  gsl_vector * g;             /* gradient J^T f */
  gsl_matrix * J;             /* Jacobian J(x) */
  gsl_vector * sqrt_wts_work; /* sqrt(W) */
  gsl_vector * sqrt_wts;      /* ptr to sqrt_wts_work, or NULL if not using weights */
  size_t niter;               /* number of iterations performed */
  gsl_multifit_nlinear_parameters params;
  void *state;
} gsl_multifit_nlinear_workspace;

gsl_multifit_nlinear_workspace *
gsl_multifit_nlinear_alloc (const gsl_multifit_nlinear_type * T,
                            const gsl_multifit_nlinear_parameters * params,
                            size_t n, size_t p);

void gsl_multifit_nlinear_free (gsl_multifit_nlinear_workspace * w);

gsl_multifit_nlinear_parameters gsl_multifit_nlinear_default_parameters(void);

int
gsl_multifit_nlinear_init (const gsl_vector * x,
                           gsl_multifit_nlinear_fdf * fdf,
                           gsl_multifit_nlinear_workspace * w);

int gsl_multifit_nlinear_winit (const gsl_vector * x,
                                const gsl_vector * wts,
                                gsl_multifit_nlinear_fdf * fdf,
                                gsl_multifit_nlinear_workspace * w);

int
gsl_multifit_nlinear_iterate (gsl_multifit_nlinear_workspace * w);

double
gsl_multifit_nlinear_avratio (const gsl_multifit_nlinear_workspace * w);

int
gsl_multifit_nlinear_driver (const size_t maxiter,
                             const double xtol,
                             const double gtol,
                             const double ftol,
                             void (*callback)(const size_t iter, void *params,
                                              const gsl_multifit_nlinear_workspace *w),
                             void *callback_params,
                             int *info,
                             gsl_multifit_nlinear_workspace * w);

gsl_matrix *
gsl_multifit_nlinear_jac (const gsl_multifit_nlinear_workspace * w);

const char *
gsl_multifit_nlinear_name (const gsl_multifit_nlinear_workspace * w);

gsl_vector *
gsl_multifit_nlinear_position (const gsl_multifit_nlinear_workspace * w);

gsl_vector *
gsl_multifit_nlinear_residual (const gsl_multifit_nlinear_workspace * w);

size_t
gsl_multifit_nlinear_niter (const gsl_multifit_nlinear_workspace * w);

int
gsl_multifit_nlinear_rcond (double *rcond, const gsl_multifit_nlinear_workspace * w);

const char *
gsl_multifit_nlinear_trs_name (const gsl_multifit_nlinear_workspace * w);

int gsl_multifit_nlinear_eval_f(gsl_multifit_nlinear_fdf *fdf,
                                const gsl_vector *x,
                                const gsl_vector *swts,
                                gsl_vector *y);

int gsl_multifit_nlinear_eval_df(const gsl_vector *x,
                                 const gsl_vector *f,
                                 const gsl_vector *swts,
                                 const double h,
                                 const gsl_multifit_nlinear_fdtype fdtype,
                                 gsl_multifit_nlinear_fdf *fdf,
                                 gsl_matrix *df, gsl_vector *work);

int
gsl_multifit_nlinear_eval_fvv(const double h,
                              const gsl_vector *x,
                              const gsl_vector *v,
                              const gsl_vector *f,
                              const gsl_matrix *J,
                              const gsl_vector *swts,
                              gsl_multifit_nlinear_fdf *fdf,
                              gsl_vector *yvv, gsl_vector *work);

/* covar.c */
int
gsl_multifit_nlinear_covar (const gsl_matrix * J, const double epsrel,
                            gsl_matrix * covar);

/* convergence.c */
int
gsl_multifit_nlinear_test (const double xtol, const double gtol,
                           const double ftol, int *info,
                           const gsl_multifit_nlinear_workspace * w);

/* fdjac.c */
int
gsl_multifit_nlinear_df(const double h, const gsl_multifit_nlinear_fdtype fdtype,
                        const gsl_vector *x, const gsl_vector *wts,
                        gsl_multifit_nlinear_fdf *fdf,
                        const gsl_vector *f, gsl_matrix *J, gsl_vector *work);

/* fdfvv.c */
int
gsl_multifit_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
                           const gsl_vector *f, const gsl_matrix *J,
                           const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
                           gsl_vector *fvv, gsl_vector *work);

/* top-level algorithms */
GSL_VAR const gsl_multifit_nlinear_type * gsl_multifit_nlinear_trust;

/* trust region subproblem methods */
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lm;
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lmaccel;
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_dogleg;
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_ddogleg;
GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_subspace2D;

/* scaling matrix strategies */
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_levenberg;
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_marquardt;
GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_more;

/* linear solvers */
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_cholesky;
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_qr;
GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_svd;

__END_DECLS

#endif /* __GSL_MULTIFIT_NLINEAR_H__ */

max_stars_repo_name: Tjoppen/fmigo

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * 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.
 *
 */
/*
 *
 *
 *       Filename:  initial.h
 *
 *    Description:  Initial values
 *
 *        Version:  1.0
 *        Created:  16/05/2014 00:44:34
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <gsl/gsl_const_mksa.h>

#define	BOLTZ	GSL_CONST_MKSA_BOLTZMANN	/* Boltzmann constant  */
#define HBAR	GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR /* hbar */

const double omega_c = 1000 ;			/* critical ohmic frequency */
const double alpha = 5e-3 ;			/* coupling strength */
const double Delta = 8.0e+9 ;			/* pumping amplitude (GHz) */

double T = .1 ;				/*
						 *  The real temperature is given by
						 *  temp = T*HBAR*Delta/BOLTZ
						 *
						 *  HBAR*Delta/BOLTZ = 0.061
						 *
						 */

const double D = 1 ;                            /* normalized delta */
double OMEGA = 2 ;				/* normalized pumping frequency */

const double gamma0 = 0.05 ;			/* energy hopping between sites */

const double t_end = 200 ;			/* time end */
const double STEP = .01 ;			/* time step */

const double R[] = { 1, 0, 0.5, -0.4 } ;   /* initial state: |z,-> */

/* const double R[] = { 1, 0, 0.5, -0.4 } ; */ /* initial state with neg. e.p. */
/* const double r[] = { 1, 0, 1, 0 } ;  	 initial state with pos. t.d. */


text:
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_sf_gamma.h>

#include "allvars.h"
#include "proto.h"
#ifdef COSMIC_RAYS
#include "cosmic_rays.h"
#endif
#ifdef CS_MODEL
#include "cs_metals.h"
#endif

#ifdef MACHNUM
#ifdef COSMIC_RAYS
#define h  All.HubbleParam
#define cm (h/All.UnitLength_in_cm)
#define s  (h/All.UnitTime_in_s)
#define LightSpeed (2.9979e10*cm/s)
#define c2   ( LightSpeed * LightSpeed )
#endif
#endif

/*! \file init.c
 *  \brief code for initialisation of a simulation from initial conditions
 */

/*! This function reads the initial conditions, and allocates storage for the
 *  tree(s). Various variables of the particle data are initialised and An
 *  intial domain decomposition is performed. If SPH particles are present,
 *  the inial SPH smoothing lengths are determined.
 */
void init(void)
{
  int i, j;
  double a3, atime;
#ifdef COSMIC_RAYS
  int CRpop;
#endif

#if defined(COSMIC_RAYS) && defined(MACHNUM)
  double Pth1, PCR1[NUMCRPOP], rBeta[NUMCRPOP], C_phys[NUMCRPOP], q_phys[NUMCRPOP];
#endif
#ifdef CR_INITPRESSURE
  double cr_pressure, q_phys, C_phys[NUMCRPOP];
#endif
#ifdef BLACK_HOLES
  int count_holes = 0;
#endif
#ifdef CHEMISTRY
  int ifunc;
  double min_t_cool, max_t_cool;
  double min_t_elec, max_t_elec;
  double a_start, a_end;
#endif

#ifdef EULERPOTENTIALS
  double a0, a1, a2;
  double b0, b1, b2;
#endif

#ifdef START_WITH_EXTRA_NGBDEV
  double MaxNumNgbDeviationMerk;
#endif

#ifdef DISTORTIONTENSORPS
  int i1, i2;
#endif

  All.Time = All.TimeBegin;

  if(RestartFlag == 3 && RestartSnapNum < 0)
    {
      if(ThisTask == 0)
	printf("Need to give the snapshot number if FOF/SUBFIND is selected for output\n");
      endrun(0);
    }

  if(RestartFlag == 4 && RestartSnapNum < 0)
    {
      if(ThisTask == 0)
	printf("Need to give the snapshot number if snapshot should be converted\n");
      endrun(0);
    }

  switch (All.ICFormat)
    {
    case 1:
    case 2:
    case 3:
    case 4:
      if(RestartFlag >= 2 && RestartSnapNum >= 0)
	{
	  char fname[1000];

	  if(All.NumFilesPerSnapshot > 1)
	    sprintf(fname, "%s/snapdir_%03d/%s_%03d", All.OutputDir, RestartSnapNum, All.SnapshotFileBase,
		    RestartSnapNum);
	  else
	    sprintf(fname, "%s%s_%03d", All.OutputDir, All.SnapshotFileBase, RestartSnapNum);
	  read_ic(fname);
	  if(RestartFlag == 4)
	    {
	      sprintf(All.SnapshotFileBase, "%s_converted", All.SnapshotFileBase);
	      if(ThisTask == 0)
		printf("Start writing file %s\n", All.SnapshotFileBase);
	      printf("RestartSnapNum %d\n", RestartSnapNum);
	      savepositions(RestartSnapNum);
	      endrun(0);
	    }

	}
      else
	{
	  if(All.ICFormat == 4)
	    read_ic_cluster(All.InitCondFile);
	  else
	    read_ic(All.InitCondFile);
	}
      break;
    default:
      if(ThisTask == 0)
	printf("ICFormat=%d not supported.\n", All.ICFormat);
      endrun(0);
    }

  All.Time = All.TimeBegin;

#ifdef COOLING
#ifdef CS_MODEL
  if(RestartFlag == 0)
    XH = HYDROGEN_MASSFRAC;
#endif
  IonizeParams();
#endif

#ifdef CHEMISTRY
  InitChem();
#endif

  if(All.ComovingIntegrationOn)
    {
      All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
      All.Ti_Current = 0;
      a3 = All.Time * All.Time * All.Time;
      atime = All.Time;
    }
  else
    {
      All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
      All.Ti_Current = 0;
      a3 = 1;
      atime = 1;
    }

#ifdef RADTRANSFER
  All.Radiation_Ti_begstep = 0;
#endif

  set_softenings();

  All.NumCurrentTiStep = 0;	/* setup some counters */
  All.SnapshotFileCount = 0;
  if(RestartFlag == 2)
    {
      if(RestartSnapNum < 0)
	All.SnapshotFileCount = atoi(All.InitCondFile + strlen(All.InitCondFile) - 3) + 1;
      else
	All.SnapshotFileCount = RestartSnapNum + 1;
    }

#ifdef OUTPUTLINEOFSIGHT
  All.Ti_nextlineofsight = (int) (log(All.TimeFirstLineOfSight / All.TimeBegin) / All.Timebase_interval);
  if(RestartFlag == 2)
    endrun(78787);
#endif

  All.TotNumOfForces = 0;
  All.NumForcesSinceLastDomainDecomp = 0;
#if defined(MAGNETIC) && defined(BSMOOTH)
#ifdef SETMAINTIMESTEPCOUNT
  All.MainTimestepCounts = All.MainTimestepCountIni;
#else
  All.MainTimestepCounts = 0;
#endif
#endif

  All.TopNodeAllocFactor = 0.008;
  All.TreeAllocFactor = 0.7;

  All.Cadj_Cost = 1.0e-30;
  All.Cadj_Cpu = 1.0e-3;

  if(All.ComovingIntegrationOn)
    if(All.PeriodicBoundariesOn == 1)
      check_omega();

  All.TimeLastStatistics = All.TimeBegin - All.TimeBetStatistics;
#ifdef BLACK_HOLES
  All.TimeNextBlackHoleCheck = All.TimeBegin;
#endif

#ifdef BUBBLES
  if(All.ComovingIntegrationOn)
    All.TimeOfNextBubble = 1. / (1. + All.FirstBubbleRedshift);
  else
    All.TimeOfNextBubble = All.TimeBegin + All.BubbleTimeInterval / All.UnitTime_in_Megayears;
  if(ThisTask == 0)
    printf("Initial time: %g and first bubble time %g \n", All.TimeBegin, All.TimeOfNextBubble);

  if(RestartFlag == 2 && All.TimeBegin > All.TimeOfNextBubble)
    {
      printf("Restarting from the snapshot file with the wrong FirstBubbleRedshift! \n");
      endrun(0);
    }
#endif

#ifdef MULTI_BUBBLES
  if(All.ComovingIntegrationOn)
    All.TimeOfNextBubble = 1. / (1. + All.FirstBubbleRedshift);
  else
    All.TimeOfNextBubble = All.TimeBegin + All.BubbleTimeInterval / All.UnitTime_in_Megayears;
  if(ThisTask == 0)
    printf("Initial time: %g and time of the first bubbles %g \n", All.TimeBegin, All.TimeOfNextBubble);
  if(RestartFlag == 2 && All.TimeBegin > All.TimeOfNextBubble)
    {
      printf("Restarting from the snapshot file with the wrong FirstBubbleRedshift! \n");
      endrun(0);
    }
#endif

  if(All.ComovingIntegrationOn)	/*  change to new velocity variable */
    {
      for(i = 0; i < NumPart; i++)
	for(j = 0; j < 3; j++)
	  P[i].Vel[j] *= sqrt(All.Time) * All.Time;
    }

#ifdef REINIT_AT_TURNAROUND
  double turnaround_radius_local=0.0, turnaround_radius_global=0.0, v_part, r_part;

   /* get local turnaroundradius */
  for (i = 0; i < NumPart; i++)
    {
     r_part = sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]);
     v_part = (P[i].Pos[0]*P[i].Vel[0]+P[i].Pos[1]*P[i].Vel[1]+P[i].Pos[2]*P[i].Vel[2]);
     if ((v_part < 0.0) && (r_part > turnaround_radius_local)) turnaround_radius_local = r_part;
    }

  /* find global turnaround radius by taking maximum of all CPUs */
  MPI_Allreduce(&turnaround_radius_local, &turnaround_radius_global, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

  All.CurrentTurnaroundRadius = turnaround_radius_global;

  if (ThisTask == 0)
   {
    printf("REINIT_AT_TURNAROUND: current (initial) turnaround radius = %g\n", All.CurrentTurnaroundRadius);
    fflush(stdout);
   }
#endif

  for(i = 0; i < NumPart; i++)	/*  start-up initialization */
    {
      for(j = 0; j < 3; j++)
	P[i].g.GravAccel[j] = 0;

#ifdef DISTORTIONTENSORPS
      /*init tidal tensor for first output (this is not used for any calculation) */
      for(i1 = 0; i1 < 3; i1++)
	for(i2 = 0; i2 < 3; i2++)
	  P[i].tidal_tensorps[i1][i2] = 0.0;

      /* find caustics on the fly by sign analysis of configuration space distortion */
      P[i].last_stream_determinant = 1.0;
#ifdef DISTORTION_READALL
      MyDouble product_matrix[3][3];
      product_matrix[0][0] = P[i].distortion_tensorps[0][0] +
                             P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][0] +
                             P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][0] +
                             P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][0];
      product_matrix[0][1] = P[i].distortion_tensorps[0][1] +
                             P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][1] +
                             P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][1] +
                             P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][1];
      product_matrix[0][2] = P[i].distortion_tensorps[0][2] +
                             P[i].distortion_tensorps[0][3] * P[i].V_matrix[0][2] +
                             P[i].distortion_tensorps[0][4] * P[i].V_matrix[1][2] +
                             P[i].distortion_tensorps[0][5] * P[i].V_matrix[2][2];
      product_matrix[1][0] = P[i].distortion_tensorps[1][0] +
                             P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][0] +
                             P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][0] +
                             P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][0];
      product_matrix[1][1] = P[i].distortion_tensorps[1][1] +
                             P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][1] +
                             P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][1] +
                             P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][1];
      product_matrix[1][2] = P[i].distortion_tensorps[1][2] +
                             P[i].distortion_tensorps[1][3] * P[i].V_matrix[0][2] +
                             P[i].distortion_tensorps[1][4] * P[i].V_matrix[1][2] +
                             P[i].distortion_tensorps[1][5] * P[i].V_matrix[2][2];
      product_matrix[2][0] = P[i].distortion_tensorps[2][0] +
                             P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][0] +
                             P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][0] +
                             P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][0];
      product_matrix[2][1] = P[i].distortion_tensorps[2][1] +
                             P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][1] +
                             P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][1] +
                             P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][1];
      product_matrix[2][2] = P[i].distortion_tensorps[2][2] +
                             P[i].distortion_tensorps[2][3] * P[i].V_matrix[0][2] +
                             P[i].distortion_tensorps[2][4] * P[i].V_matrix[1][2] +
                             P[i].distortion_tensorps[2][5] * P[i].V_matrix[2][2];

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

#endif
#ifdef ANNIHILATION_RADIATION
      /* integrated annihilation rate */
      P[i].annihilation = 0.0;

      P[i].rho_normed_cutoff_current = 1.0;
      P[i].rho_normed_cutoff_last = 1.0;

      P[i].stream_density = P[i].init_density;

      P[i].analytic_caustics = 0.0;
#endif

#ifdef REINIT_AT_TURNAROUND
      All.SIM_epsilon = REINIT_AT_TURNAROUND;
      /* no phase-space analysis for particles that are initially within the turnaround radius */
      if (sqrt(P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2]) < All.CurrentTurnaroundRadius)
#ifdef DISTORTION_READALL
       P[i].turnaround_flag = 1;
#else
       P[i].turnaround_flag = -1;
#endif
      else
       P[i].turnaround_flag = 0;
#endif

#ifdef OUTPUT_LAST_CAUSTIC
      /* all entries zero -> no caustic yet */
      P[i].lc_Time       = 0.0;
      P[i].lc_Pos[0]     = 0.0;
      P[i].lc_Pos[1]     = 0.0;
      P[i].lc_Pos[2]     = 0.0;
      P[i].lc_Vel[0]     = 0.0;
      P[i].lc_Vel[1]     = 0.0;
      P[i].lc_Vel[2]     = 0.0;
      P[i].lc_rho_normed_cutoff = 0.0;

      P[i].lc_Dir_x[0]   = 0.0;
      P[i].lc_Dir_x[1]   = 0.0;
      P[i].lc_Dir_x[2]   = 0.0;
      P[i].lc_Dir_y[0]   = 0.0;
      P[i].lc_Dir_y[1]   = 0.0;
      P[i].lc_Dir_y[2]   = 0.0;
      P[i].lc_Dir_z[0]   = 0.0;
      P[i].lc_Dir_z[1]   = 0.0;
      P[i].lc_Dir_z[2]   = 0.0;

      P[i].lc_smear_x    = 0.0;
      P[i].lc_smear_y    = 0.0;
      P[i].lc_smear_z    = 0.0;
#endif

#ifdef PMGRID
      /* long range tidal field init */
      P[i].tidal_tensorpsPM[0][0] = 0;
      P[i].tidal_tensorpsPM[0][1] = 0;
      P[i].tidal_tensorpsPM[0][2] = 0;
      P[i].tidal_tensorpsPM[1][0] = 0;
      P[i].tidal_tensorpsPM[1][1] = 0;
      P[i].tidal_tensorpsPM[1][2] = 0;
      P[i].tidal_tensorpsPM[2][0] = 0;
      P[i].tidal_tensorpsPM[2][1] = 0;
      P[i].tidal_tensorpsPM[2][2] = 0;
#endif

#ifndef DISTORTION_READALL
      /*6D phase space distortion tensor equals unity in the beginning */

      /*
       * Distortion tensor D' (comoving) equals unity in the beginning
       * So the initial stream density is given by |det(D')|*mean density
       */

      /* cold dark matter => velocity dispersion=0 => d/dtD' = 0 => Z'=0 (comoving) */
      for(i1 = 0; i1 < 6; i1++)
	for(i2 = 0; i2 < 6; i2++)
	  {
	    if((i1 == i2))
	      {
		P[i].distortion_tensorps[i1][i2] = 1.0;
	      }
	    else
	      {
		P[i].distortion_tensorps[i1][i2] = 0.0;
	      }
	  }
#endif
#endif

#ifdef KEEP_DM_HSML_AS_GUESS
      if(RestartFlag != 1)
	P[i].DM_Hsml = -1;
#endif

#ifdef PMGRID
      for(j = 0; j < 3; j++)
	P[i].GravPM[j] = 0;
#endif
      P[i].Ti_begstep = 0;
      P[i].Ti_current = 0;
      P[i].TimeBin = 0;

      if(header.flag_ic_info != FLAG_SECOND_ORDER_ICS)
	P[i].OldAcc = 0;	/* Do not zero in 2lpt case as masses are stored here */

      P[i].GravCost = 1;
#if defined(EVALPOTENTIAL) || defined(COMPUTE_POTENTIAL_ENERGY)
      P[i].p.Potential = 0;
#endif
#ifdef STELLARAGE
      if(RestartFlag == 0)
	P[i].StellarAge = 0;
#endif

#ifdef METALS
#ifndef CS_MODEL
      if(RestartFlag == 0)
	P[i].Metallicity = 0;
#else
      if(RestartFlag == 0)
	{
	  for(j = 0; j < 12; j++)
	    {
	      P[i].Zm[j] = 0;
	      P[i].ZmReservoir[j] = 0;
	    }

	  P[i].Zm[6] = HYDROGEN_MASSFRAC * (P[i].Mass);
	  P[i].Zm[0] = (1 - HYDROGEN_MASSFRAC) * (P[i].Mass);
	  /*     Order of chemical elements:   He, Carbon,Mg,O,Fe,Si,H,N,Ne,S,Ca,Zn */

	  PPP[i].Hsml = 0;
#ifdef CS_FEEDBACK
	  PPP[i].EnergySN = 0;
	  PPP[i].EnergySNCold = 0;
#endif
	}
#endif

#endif

#ifdef BLACK_HOLES
      if(P[i].Type == 5)
	{
	  count_holes++;

	  if(RestartFlag == 0)
	    P[i].BH_Mass = All.SeedBlackHoleMass;
#ifdef BH_BUBBLES
	  if(RestartFlag == 0)
	    {
	      P[i].BH_Mass_bubbles = All.SeedBlackHoleMass;
	      P[i].BH_Mass_ini = All.SeedBlackHoleMass;
#ifdef UNIFIED_FEEDBACK
	      P[i].BH_Mass_radio = All.SeedBlackHoleMass;
#endif
	    }
#endif
	}
#endif
    }

#ifdef BLACK_HOLES
  MPI_Allreduce(&count_holes, &All.TotBHs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#endif

  for(i = 0; i < TIMEBINS; i++)
    TimeBinActive[i] = 1;

  reconstruct_timebins();

#ifdef PMGRID
  All.PM_Ti_endstep = All.PM_Ti_begstep = 0;
#endif

#ifdef CONDUCTION
  All.Conduction_Ti_endstep = All.Conduction_Ti_begstep = 0;
#endif
#ifdef CR_DIFFUSION
  All.CR_Diffusion_Ti_endstep = All.CR_Diffusion_Ti_begstep = 0;
#endif

  for(i = 0; i < N_gas; i++)	/* initialize sph_properties */
    {
      for(j = 0; j < 3; j++)
	{
	  SphP[i].VelPred[j] = P[i].Vel[j];
	  SphP[i].a.HydroAccel[j] = 0;
	}

      SphP[i].e.DtEntropy = 0;

#ifdef CHEMISTRY
      SphP[i].Gamma = GAMMA;	/* set universal value */
      SphP[i].t_cool = 0;
      SphP[i].t_elec = 0;
#endif

#ifdef CS_MODEL
      SphP[i].DensityOld = 0;
#endif
#ifdef CS_FEEDBACK
      SphP[i].da.DensityAvg = 0;
      SphP[i].ea.EntropyAvg = 0;
      SphP[i].DensPromotion = 0;
      SphP[i].TempPromotion = 0;
      SphP[i].HotHsml = 0;
#endif

      if(RestartFlag == 0)
	{
#ifndef READ_HSML
	  PPP[i].Hsml = 0;
#endif
	  SphP[i].d.Density = -1;
#ifdef VOLUME_CORRECTION
	  SphP[i].DensityOld = 1;
#endif
#ifdef COOLING
	  SphP[i].Ne = 1.0;
#endif
	  SphP[i].v.DivVel = 0;
	}
#ifdef WINDS
      SphP[i].DelayTime = 0;
#endif
#ifdef SFR
      SphP[i].Sfr = 0;
#endif
#ifdef MAGNETIC
#ifdef BINISET
      SphP[i].BPred[0] = All.BiniX;
      SphP[i].BPred[1] = All.BiniY;
      SphP[i].BPred[2] = All.BiniZ;
#ifdef EULERPOTENTIALS
      if(All.BiniY == 0 && All.BiniZ == 0)
	{
	  a0 = 0;
	  a1 = All.BiniX;
	  a2 = 0;
	  b0 = 0;
	  b1 = 1;
	  b2 = 1;
	}
      else
	{
	  if(All.BiniX != 0 && (All.BiniY != 0 || All.BiniZ != 0))
	    {
	      b0 = -(All.BiniZ + All.BiniY) / All.BiniX;
	      b1 = 1;
	      b2 = 1;
	      a0 = 0;
	      a1 = -All.BiniZ / b0;
	      a2 = All.BiniY / b0;
	    }
	  else
	    {
	      a0 = a1 = a2 = b0 = b1 = b2 = 0;
	      if(ThisTask == 0)
		{
		  printf("Can not reconstruct Euler potentials from Bini values !\n");
		  endrun(6723);
		}
	    }
	}
      SphP[i].EulerA = (a0 * P[i].Pos[0] + a1 * P[i].Pos[1] + a2 * P[i].Pos[2]) * atime * All.HubbleParam;
      SphP[i].EulerB = (b0 * P[i].Pos[0] + b1 * P[i].Pos[1] + b2 * P[i].Pos[2]) * atime * All.HubbleParam;
#endif
#endif
#ifndef EULERPOTENTIALS
      for(j = 0; j < 3; j++)
	{
	  SphP[i].DtB[j] = 0;
	  SphP[i].B[j] = SphP[i].BPred[j];
	}
#endif
#ifdef TIME_DEP_MAGN_DISP
#ifdef HIGH_MAGN_DISP_START
      SphP[i].Balpha = All.ArtMagDispConst;
#else
      SphP[i].Balpha = All.ArtMagDispMin;
#endif
      SphP[i].DtBalpha = 0.0;
#endif
#ifdef DIVBCLEANING_DEDNER
      SphP[i].Phi = SphP[i].PhiPred = SphP[i].DtPhi = 0;
#endif
#endif

#ifdef TIME_DEP_ART_VISC
#ifdef HIGH_ART_VISC_START
      if(HIGH_ART_VISC_START == 0)
	SphP[i].alpha = All.ArtBulkViscConst;
      if(HIGH_ART_VISC_START > 0)
	if(P[i].Pos[0] > HIGH_ART_VISC_START)
	  SphP[i].alpha = All.ArtBulkViscConst;
	else
	  SphP[i].alpha = All.AlphaMin;
      if(HIGH_ART_VISC_START < 0)
	if(P[i].Pos[0] < -HIGH_ART_VISC_START)
	  SphP[i].alpha = All.ArtBulkViscConst;
	else
	  SphP[i].alpha = All.AlphaMin;
#else
      SphP[i].alpha = All.AlphaMin;
#endif
      SphP[i].Dtalpha = 0.0;
#endif

#if defined(BH_THERMALFEEDBACK) || defined(BH_KINETICFEEDBACK)
      SphP[i].i.Injected_BH_Energy = 0;
#endif
    }

#ifdef TWODIMS
  for(i = 0; i < NumPart; i++)
    {
      P[i].Pos[2] = 0;
      P[i].Vel[2] = 0;
    }
#endif

#ifdef ASSIGN_NEW_IDS
  assign_unique_ids();
#endif

#ifndef NOTEST_FOR_IDUNIQUENESS
  test_id_uniqueness();
#endif

  All.NumForcesSinceLastDomainDecomp = (long long) (1 + All.TotNumPart * All.TreeDomainUpdateFrequency);

  Flag_FullStep = 1;		/* to ensure that Peano-Hilber order is done */

  TreeReconstructFlag = 1;

  domain_Decomposition();	/* do initial domain decomposition (gives equal numbers of particles) */

  set_softenings();

  /* will build tree */
  ngb_treebuild();

  All.Ti_Current = 0;

#ifdef START_WITH_EXTRA_NGBDEV
  MaxNumNgbDeviationMerk = All.MaxNumNgbDeviation;
  All.MaxNumNgbDeviation = All.MaxNumNgbDeviationStart;
#endif

  if(RestartFlag != 3)
  setup_smoothinglengths();

#ifdef VORONOI
  for(i = 0; i < N_gas; i++)
    SphP[i].MaxDelaunayRadius = SphP[i].Hsml;

  voronoi_mesh();

#ifdef MESHRELAX_DENSITY_IN_INPUT
  if(RestartFlag == 0)
    for(i = 0; i < N_gas; i++)
      P[i].Mass *= SphP[i].Volume;
#endif

#ifdef VORONOI_MESHRELAX
  double mass, masstot, egy, egytot;
  for(i = 0, mass = 0, egy = 0; i < N_gas; i++)
    {
      mass += P[i].Mass;
      egy += P[i].Mass * SphP[i].Entropy;
      SphP[i].Pressure = GAMMA_MINUS1 * SphP[i].Entropy * P[i].Mass / SphP[i].Volume;
   }
  MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&egy, &egytot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  All.MeanMass = masstot / All.TotN_gas;
#ifdef TWODIMS
  All.MeanPressure = GAMMA_MINUS1 * egytot / (boxSize_X * boxSize_Y);
#else
  All.MeanPressure = GAMMA_MINUS1 * egytot / (boxSize_X * boxSize_Y * boxSize_Z);
#endif
#endif

  voronoi_density();

  myfree(List_P);
  myfree(ListExports);

  myfree(DT);
  myfree(DP - 5);
  myfree(VF);
#endif

#ifdef START_WITH_EXTRA_NGBDEV
  All.MaxNumNgbDeviation = MaxNumNgbDeviationMerk;
#endif

#ifdef HEALPIX
  //this should be readed in the parameterfile
  All.Nside=8;
  //
  if(ThisTask==0)
  printf(" First calculation of Healpix %i with %i \n",All.Nside,NSIDE2NPIX(All.Nside));
  // initialize the healpix array (just in case)
  All.healpixmap=(float *) malloc(NSIDE2NPIX(All.Nside) * sizeof(float));
  for(i=0;i<NSIDE2NPIX(All.Nside);i++)All.healpixmap[i]=0;

  healpix_halo(All.healpixmap);

#endif
  /* at this point, the entropy variable actually contains the
   * internal energy, read in from the initial conditions file.
   * Once the density has been computed, we can convert to entropy.
   */

  for(i = 0; i < N_gas; i++)	/* initialize sph_properties */
    {
      if(header.flag_entropy_instead_u == 0)
	{
	  if(ThisTask == 0 && i == 0)
	    printf("Converting u -> entropy !\n");
#ifndef EOS_DEGENERATE
#ifndef VORONOI_MESHRELAX
	  SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].d.Density / a3, GAMMA_MINUS1);
#endif
#else
	  SphP[i].Entropy *= All.UnitEnergy_in_cgs;
	  /* call eos with physical units, energy and entropy are always stored in physical units */
#ifdef EOS_ENERGY
	  SphP[i].temp = 1.0;
	  eos_calc_egiven2(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc, SphP[i].Entropy,
			   &SphP[i].temp, &SphP[i].Pressure);
	  SphP[i].Pressure /= All.UnitPressure_in_cgs;
#else
	  SphP[i].u = SphP[i].Entropy;
	  eos_calc_egiven(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc, SphP[i].Entropy,
			  &SphP[i].temp, &SphP[i].Entropy);

	  /* get pressure */
	  eos_calc_sgiven(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc, SphP[i].Entropy,
			  &SphP[i].temp, &SphP[i].Pressure, &SphP[i].u);
	  SphP[i].Pressure /= All.UnitPressure_in_cgs;
#endif
#endif
	}
      SphP[i].e.DtEntropy = 0;

      SphP[i].v.DivVel = 0;

#ifdef MACHNUM
      SphP[i].Shock_MachNumber = 1.0;
#ifdef COSMIC_RAYS
      Pth1 = SphP[i].Entropy * pow(SphP[i].d.Density / a3, GAMMA);

#ifdef CR_IC_PHYSICAL
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	{
	  C_phys[CRpop] = SphP[i].CR_C0[CRpop];
	  q_phys[CRpop] = SphP[i].CR_q0[CRpop];
	}
#else
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	{
	  C_phys[CRpop] = SphP[i].CR_C0[CRpop] * pow(SphP[i].d.Density, (All.CR_Alpha[CRpop] - 1.0) / 3.0);
	  q_phys[CRpop] = SphP[i].CR_q0[CRpop] * pow(SphP[i].d.Density, 1.0 / 3.0);
	}
#endif
      SphP[i].PreShock_XCR = 0.0;
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	{
	  rBeta[CRpop] = gsl_sf_beta((All.CR_Alpha[CRpop] - 2.0) * 0.5, (3.0 - All.CR_Alpha[CRpop]) * 0.5) *
	    gsl_sf_beta_inc((All.CR_Alpha[CRpop] - 2.0) * 0.5, (3.0 - All.CR_Alpha[CRpop]) * 0.5,
			    1.0 / (1.0 + q_phys[CRpop] * q_phys[CRpop]));

	  PCR1[CRpop] = C_phys[CRpop] * c2 * SphP[i].d.Density * rBeta[CRpop] / 6.0;
	  PCR1[CRpop] *= pow(atime, -3.0 * GAMMA);
	  SphP[i].PreShock_XCR += PCR1[CRpop] / Pth1;
	}

      SphP[i].PreShock_PhysicalDensity = SphP[i].d.Density / a3;
      SphP[i].PreShock_PhysicalEnergy =
	SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].d.Density / a3, GAMMA_MINUS1);

      SphP[i].Shock_DensityJump = 1.0001;
      SphP[i].Shock_EnergyJump = 1.0;
#endif /* COSMIC_RAYS */
#ifdef OUTPUT_PRESHOCK_CSND
      Pth1 = SphP[i].Entropy * pow(SphP[i].d.Density / a3, GAMMA);
      SphP[i].PreShock_PhysicalSoundSpeed =
	sqrt(GAMMA * Pth1 / SphP[i].d.Density) * pow(atime, -3. / 2. * GAMMA_MINUS1);
      SphP[i].PreShock_PhysicalDensity = SphP[i].d.Density / a3;
#endif /* OUTPUT_PRESHOCK_CSND */
#endif /* MACHNUM */

#ifdef REIONIZATION
      All.not_yet_reionized = 1;
#endif

#ifdef CR_IC_PHYSICAL
      /* Scale CR variables so that values from IC file are now the
       * physical values, not the adiabatic invariants
       */
      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	{
	  SphP[i].CR_C0[CRpop] *= pow(SphP[i].d.Density, (1.0 - All.CR_Alpha[CRpop]) / 3.0);
	  SphP[i].CR_q0[CRpop] *= pow(SphP[i].d.Density, -1.0 / 3.0);
	}
#endif

#ifdef CR_INITPRESSURE

      cr_pressure = CR_INITPRESSURE * SphP[i].Entropy * pow(SphP[i].d.Density / a3, GAMMA);
      SphP[i].Entropy *= (1 - CR_INITPRESSURE);
      q_phys = 1.685;

      for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
	{
	  C_phys[CRpop] =
	    cr_pressure / (SphP[i].d.Density / a3 * CR_Tab_Beta(q_phys, CRpop) * (C / All.UnitVelocity_in_cm_per_s) *
		       (C / All.UnitVelocity_in_cm_per_s) / 6.0);

	  SphP[i].CR_C0[CRpop] = C_phys[CRpop] * pow(SphP[i].d.Density, (1.0 - All.CR_Alpha[CRpop]) / 3.0);
	  SphP[i].CR_q0[CRpop] = q_phys * pow(SphP[i].d.Density, -1.0 / 3.0);
	}
#endif
    }

  if(RestartFlag == 3)
    {
#ifdef FOF
      fof_fof(RestartSnapNum);
#endif
      endrun(0);
    }

#ifdef CHEMISTRY
  if(ThisTask == 0)
    {
      printf("Initial abundances: \n");
      printf("HI=%g, HII=%g, HeI=%g, HeII=%g, HeIII=%g \n",
	     SphP[1].HI, SphP[1].HII, SphP[1].HeI, SphP[1].HeII, SphP[1].HeIII);

      printf("HM=%g, H2I=%g, H2II=%g, elec=%g, %d\n",
	     SphP[1].HM, SphP[1].H2I, SphP[1].H2II, SphP[1].elec, P[1].ID);

      printf("x=%g, y=%g, z=%g, vx=%g, vy=%g, vz=%g, density=%g, entropy=%g\n",
	     P[N_gas - 1].Pos[0], P[N_gas - 1].Pos[1], P[N_gas - 1].Pos[2], P[N_gas - 1].Vel[0],
	     P[N_gas - 1].Vel[1], P[N_gas - 1].Vel[2], SphP[N_gas - 1].Density, SphP[N_gas - 1].Entropy);
    }

  /* need predict the cooling time and elec_dot here */
  min_t_cool = min_t_elec = 1.0e30;
  max_t_cool = max_t_elec = -1.0e30;

  for(i = 0; i < N_gas; i++)
    {
      a_start = All.Time;
      a_end = All.Time + 0.001;	/* 0.001 as an arbitrary value */

      ifunc = compute_abundances(0, i, a_start, a_end);

      if(fabs(SphP[i].t_cool) < min_t_cool)
	min_t_cool = fabs(SphP[i].t_cool);
      if(fabs(SphP[i].t_cool) > max_t_cool)
	max_t_cool = fabs(SphP[i].t_cool);

      if(fabs(SphP[i].t_elec) < min_t_elec)
	min_t_elec = fabs(SphP[i].t_elec);
      if(fabs(SphP[i].t_elec) > max_t_elec)
	max_t_elec = fabs(SphP[i].t_elec);

    }

  fprintf(stdout, "PE %d t_cool min= %g, max= %g in yrs \n", ThisTask, min_t_cool, max_t_cool);
  fflush(stdout);
  fprintf(stdout, "PE %d t_elec min= %g, max= %g in yrs \n", ThisTask, min_t_elec, max_t_elec);
  fflush(stdout);

#endif

}

/*! This routine computes the mass content of the box and compares it to the
 * specified value of Omega-matter.  If discrepant, the run is terminated.
 */
void check_omega(void)
{
  double mass = 0, masstot, omega;
  int i;

  for(i = 0; i < NumPart; i++)
    mass += P[i].Mass;

  MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  omega =
    masstot / (All.BoxSize * All.BoxSize * All.BoxSize) / (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
#ifdef TIMEDEPGRAV
  omega *= All.Gini / All.G;
#endif
  if(fabs(omega - All.Omega0) > 1.0e-3)
    {
      if(ThisTask == 0)
	{
	  printf("\n\nI've found something odd!\n");
	  printf
	    ("The mass content accounts only for Omega=%g,\nbut you specified Omega=%g in the parameterfile.\n",
	     omega, All.Omega0);
	  printf("\nI better stop.\n");

	  fflush(stdout);
	}
      endrun(1);
    }
}

/*! This function is used to find an initial smoothing length for each SPH
 *  particle. It guarantees that the number of neighbours will be between
 *  desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
 *  of the smoothing length is provided to the function density(), which will
 *  then iterate if needed to find the right smoothing length.
 */
void setup_smoothinglengths(void)
{
  int i, no, p;

  if(RestartFlag == 0)
    {
      for(i = 0; i < N_gas; i++)
	{
	  no = Father[i];

	  while(10 * All.DesNumNgb * P[i].Mass > Nodes[no].u.d.mass)
	    {
	      p = Nodes[no].u.d.father;

	      if(p < 0)
		break;

	      no = p;
	    }

#ifndef READ_HSML
#ifndef TWODIMS
	  PPP[i].Hsml =
	    pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
#else
	  PPP[i].Hsml =
	    pow(1.0 / (M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 2) * Nodes[no].len;
#endif
	  if(PPP[i].Hsml > 200.0 * All.SofteningTable[0])
	    PPP[i].Hsml = All.SofteningTable[0];
#endif
	}
    }

#ifdef BLACK_HOLES
  if(RestartFlag == 0 || RestartFlag == 2)
    {
      for(i = 0; i < NumPart; i++)
	if(P[i].Type == 5)
	  PPP[i].Hsml = All.SofteningTable[5];
    }
#endif

#ifdef RADTRANSFER
  if(RestartFlag == 0 || RestartFlag == 2)
    {
      for(i = 0; i < NumPart; i++)
	if(P[i].Type == 4)
	  PPP[i].Hsml = All.SofteningTable[4];
    }
#endif

  density();

#if defined(MAGNETIC) && defined(BFROMROTA)
  if(RestartFlag == 0)
    {
      if(ThisTask == 0)
	printf("Converting: Vector Potential -> Bfield\n");
      rot_a();
    }
#endif

}

void assign_unique_ids(void)
{
  int i, *numpartlist;
  MyIDType idfirst;

  numpartlist = mymalloc(NTask * sizeof(int));

  MPI_Allgather(&NumPart, 1, MPI_INT, numpartlist, 1, MPI_INT, MPI_COMM_WORLD);

  idfirst = 1;

  for(i = 0; i < ThisTask; i++)
    idfirst += numpartlist[i];

  for(i = 0; i < NumPart; i++)
    {
      P[i].ID = idfirst;
      idfirst++;
    }

  myfree(numpartlist);
}

void test_id_uniqueness(void)
{
  int i;
  double t0, t1;
  MyIDType *ids, *ids_first;

  if(ThisTask == 0)
    {
      printf("Testing ID uniqueness...\n");
      fflush(stdout);
    }

  if(NumPart == 0)
    {
      printf("need at least one particle per cpu\n");
      endrun(8);
    }

  t0 = second();

#ifndef SPH_BND_PARTICLES
  ids = (MyIDType *) mymalloc(NumPart * sizeof(MyIDType));
  ids_first = (MyIDType *) mymalloc(NTask * sizeof(MyIDType));

  for(i = 0; i < NumPart; i++)
    ids[i] = P[i].ID;

  parallel_sort(ids, NumPart, sizeof(MyIDType), compare_IDs);

  for(i = 1; i < NumPart; i++)
    if(ids[i] == ids[i - 1])
      {
#ifdef LONGIDS
	printf("non-unique ID=%d%09d found on task=%d (i=%d NumPart=%d)\n",
	       (int) (ids[i] / 1000000000), (int) (ids[i] % 1000000000), ThisTask, i, NumPart);

#else
	printf("non-unique ID=%d found on task=%d   (i=%d NumPart=%d)\n", (int) ids[i], ThisTask, i, NumPart);
#endif
	endrun(12);
      }

  MPI_Allgather(&ids[0], sizeof(MyIDType), MPI_BYTE, ids_first, sizeof(MyIDType), MPI_BYTE, MPI_COMM_WORLD);

  if(ThisTask < NTask - 1)
    if(ids[NumPart - 1] == ids_first[ThisTask + 1])
      {
	printf("non-unique ID=%d found on task=%d\n", (int) ids[NumPart - 1], ThisTask);
	endrun(13);
      }

  myfree(ids_first);
  myfree(ids);
#endif

  t1 = second();

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

int compare_IDs(const void *a, const void *b)
{
  if(*((MyIDType *) a) < *((MyIDType *) b))
    return -1;

  if(*((MyIDType *) a) > *((MyIDType *) b))
    return +1;

  return 0;
}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp