text:
/* eigen/gsl_eigen.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Brian Gough, 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_EIGEN_H__
#define __GSL_EIGEN_H__

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

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

__BEGIN_DECLS

typedef struct {
  size_t size;
  double * d;
  double * sd;
} gsl_eigen_symm_workspace;

gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * gc;
  double * gs;
} gsl_eigen_symmv_workspace;

gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
} gsl_eigen_herm_workspace;

gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
                         gsl_eigen_herm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
  double * gc;
  double * gs;
} gsl_eigen_hermv_workspace;

gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
                           gsl_matrix_complex * evec,
                           gsl_eigen_hermv_workspace * w);

typedef struct {
  size_t size;           /* matrix size */
  size_t max_iterations; /* max iterations since last eigenvalue found */
  size_t n_iter;         /* number of iterations since last eigenvalue found */
  size_t n_evals;        /* number of eigenvalues found so far */

  int compute_t;         /* compute Schur form T = Z^t A Z */

  gsl_matrix *H;         /* pointer to Hessenberg matrix */
  gsl_matrix *Z;         /* pointer to Schur vector matrix */
} gsl_eigen_francis_workspace;

gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
void gsl_eigen_francis_T (const int compute_t,
                          gsl_eigen_francis_workspace * w);
int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
                       gsl_eigen_francis_workspace * w);
int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
                         gsl_matrix * Z,
                         gsl_eigen_francis_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *diag;            /* diagonal matrix elements from balancing */
  gsl_vector *tau;             /* Householder coefficients */
  gsl_matrix *Z;               /* pointer to Z matrix */
  int do_balance;              /* perform balancing transformation? */
  size_t n_evals;              /* number of eigenvalues found */

  gsl_eigen_francis_workspace *francis_workspace_p;
} gsl_eigen_nonsymm_workspace;

gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
                               gsl_eigen_nonsymm_workspace *w);
int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
                       gsl_eigen_nonsymm_workspace * w);
int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
                         gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *work;            /* scratch workspace */
  gsl_vector *work2;           /* scratch workspace */
  gsl_vector *work3;           /* scratch workspace */

  gsl_matrix *Z;               /* pointer to Schur vectors */

  gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
} gsl_eigen_nonsymmv_workspace;

gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
void gsl_eigen_nonsymmv_params (const int balance,
                                gsl_eigen_nonsymmv_workspace *w);
int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                        gsl_matrix_complex * evec,
                        gsl_eigen_nonsymmv_workspace * w);
int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                          gsl_matrix_complex * evec, gsl_matrix * Z,
                          gsl_eigen_nonsymmv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_symm_workspace *symm_workspace_p;
} gsl_eigen_gensymm_workspace;

gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc (const size_t n);
void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace * w);
int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
                       gsl_vector * eval, gsl_eigen_gensymm_workspace * w);
int gsl_eigen_gensymm_standardize (gsl_matrix * A, const gsl_matrix * B);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_symmv_workspace *symmv_workspace_p;
} gsl_eigen_gensymmv_workspace;

gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc (const size_t n);
void gsl_eigen_gensymmv_free (gsl_eigen_gensymmv_workspace * w);
int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
                        gsl_vector * eval, gsl_matrix * evec,
                        gsl_eigen_gensymmv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_herm_workspace *herm_workspace_p;
} gsl_eigen_genherm_workspace;

gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc (const size_t n);
void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace * w);
int gsl_eigen_genherm (gsl_matrix_complex * A, gsl_matrix_complex * B,
                       gsl_vector * eval, gsl_eigen_genherm_workspace * w);
int gsl_eigen_genherm_standardize (gsl_matrix_complex * A,
                                   const gsl_matrix_complex * B);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_eigen_hermv_workspace *hermv_workspace_p;
} gsl_eigen_genhermv_workspace;

gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc (const size_t n);
void gsl_eigen_genhermv_free (gsl_eigen_genhermv_workspace * w);
int gsl_eigen_genhermv (gsl_matrix_complex * A, gsl_matrix_complex * B,
                        gsl_vector * eval, gsl_matrix_complex * evec,
                        gsl_eigen_genhermv_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */
  gsl_vector *work;       /* scratch workspace */

  size_t n_evals;         /* number of eigenvalues found */
  size_t max_iterations;  /* maximum QZ iterations allowed */
  size_t n_iter;          /* number of iterations since last eigenvalue found */
  double eshift;          /* exceptional shift counter */

  int needtop;            /* need to compute top index? */

  double atol;            /* tolerance for splitting A matrix */
  double btol;            /* tolerance for splitting B matrix */

  double ascale;          /* scaling factor for shifts */
  double bscale;          /* scaling factor for shifts */

  gsl_matrix *H;          /* pointer to hessenberg matrix */
  gsl_matrix *R;          /* pointer to upper triangular matrix */

  int compute_s;          /* compute generalized Schur form S */
  int compute_t;          /* compute generalized Schur form T */

  gsl_matrix *Q;          /* pointer to left Schur vectors */
  gsl_matrix *Z;          /* pointer to right Schur vectors */
} gsl_eigen_gen_workspace;

gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const size_t n);
void gsl_eigen_gen_free (gsl_eigen_gen_workspace * w);
void gsl_eigen_gen_params (const int compute_s, const int compute_t,
                           const int balance, gsl_eigen_gen_workspace * w);
int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
                   gsl_vector_complex * alpha, gsl_vector * beta,
                   gsl_eigen_gen_workspace * w);
int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
                      gsl_vector_complex * alpha, gsl_vector * beta,
                      gsl_matrix * Q, gsl_matrix * Z,
                      gsl_eigen_gen_workspace * w);

typedef struct {
  size_t size;            /* size of matrices */

  gsl_vector *work1;      /* 1-norm of columns of A */
  gsl_vector *work2;      /* 1-norm of columns of B */
  gsl_vector *work3;      /* real part of eigenvector */
  gsl_vector *work4;      /* imag part of eigenvector */
  gsl_vector *work5;      /* real part of back-transformed eigenvector */
  gsl_vector *work6;      /* imag part of back-transformed eigenvector */

  gsl_matrix *Q;          /* pointer to left Schur vectors */
  gsl_matrix *Z;          /* pointer to right Schur vectors */

  gsl_eigen_gen_workspace *gen_workspace_p;
} gsl_eigen_genv_workspace;

gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const size_t n);
void gsl_eigen_genv_free (gsl_eigen_genv_workspace * w);
int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
                    gsl_vector_complex * alpha, gsl_vector * beta,
                    gsl_matrix_complex * evec,
                    gsl_eigen_genv_workspace * w);
int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
                       gsl_vector_complex * alpha, gsl_vector * beta,
                       gsl_matrix_complex * evec,
                       gsl_matrix * Q, gsl_matrix * Z,
                       gsl_eigen_genv_workspace * w);

typedef enum {
  GSL_EIGEN_SORT_VAL_ASC,
  GSL_EIGEN_SORT_VAL_DESC,
  GSL_EIGEN_SORT_ABS_ASC,
  GSL_EIGEN_SORT_ABS_DESC
}
gsl_eigen_sort_t;

/* Sort eigensystem results based on eigenvalues.
 * Sorts in order of increasing value or increasing
 * absolute value.
 *
 * exceptions: GSL_EBADLEN
 */

int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
                         gsl_eigen_sort_t sort_type);

int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
                            gsl_matrix_complex * evec,
                            gsl_eigen_sort_t sort_type);

int gsl_eigen_gensymmv_sort (gsl_vector * eval, gsl_matrix * evec,
                             gsl_eigen_sort_t sort_type);

int gsl_eigen_genhermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
                             gsl_eigen_sort_t sort_type);

int gsl_eigen_genv_sort (gsl_vector_complex * alpha, gsl_vector * beta,
                         gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

/* Prototypes for the schur module */

int gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B,
                          double *wr1, double *wr2, double *wi,
                          double *scale1, double *scale2);

int gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
                             double d1, double d2, const gsl_vector *b,
                             gsl_vector *x, double *s, double *xnorm,
                             double smin);

int gsl_schur_solve_equation_z(double ca, const gsl_matrix *A,
                               gsl_complex *z, double d1, double d2,
                               const gsl_vector_complex *b,
                               gsl_vector_complex *x, double *s,
                               double *xnorm, double smin);

/* The following functions are obsolete: */

/* Eigensolve by Jacobi Method
 *
 * The data in the matrix input is destroyed.
 *
 * exceptions:
 */
int
gsl_eigen_jacobi(gsl_matrix * matrix,
                      gsl_vector * eval,
                      gsl_matrix * evec,
                      unsigned int max_rot,
                      unsigned int * nrot);

/* Invert by Jacobi Method
 *
 * exceptions:
 */
int
gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
                             gsl_matrix * ainv,
                             unsigned int max_rot);

__END_DECLS

#endif /* __GSL_EIGEN_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


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

#include "allvars.h"
#include "proto.h"

#ifdef SUBFIND

#include "fof.h"
#include "subfind.h"

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct contamdata_in
{
  MyDouble Pos[3];
  MyFloat R200;
  int NodeList[NODELISTLENGTH];
}
 *ContamIn, *ContamGet;

static struct contamdata_out
{
  double ContaminationMass;
  int ContaminationLen;
}
 *ContamResult, *ContamOut;

void subfind_contamination(void)
{
  int i, j, ndone, ndone_flag, dummy, count;
  int ngrp, sendTask, recvTask, place, nexport, nimport;
  struct unbind_data *d;

  d = (struct unbind_data *) mymalloc(NumPart * sizeof(struct unbind_data));

  for(i = 0, count = 0; i < NumPart; i++)
#ifdef DENSITY_SPLIT_BY_TYPE
    if(!((1 << P[i].Type) & (DENSITY_SPLIT_BY_TYPE)))
#else
    if(!((1 << P[i].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
      d[count++].index = i;

  force_treebuild(count, d);	/* construct tree only with boundary particles */

  myfree(d);

  /* allocate buffers to arrange communication */

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct contamdata_in) + sizeof(struct contamdata_out) +
					     sizemax(sizeof(struct contamdata_in),
						     sizeof(struct contamdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  /* we will repeat the whole thing for those groups where we didn't converge to a SO radius yet */

  i = 0;			/* begin with this index */

  do
    {
      for(j = 0; j < NTask; j++)
	{
	  Send_count[j] = 0;
	  Exportflag[j] = -1;
	}

      /* do local particles and prepare export list */

      for(nexport = 0; i < Ngroups; i++)
	{
	  if(Group[i].R_Mean200 > 0)
	    {
	      if(subfind_contamination_evaluate(i, 0, &nexport, Send_count) < 0)
		break;
	    }
	  else
	    {
	      Group[i].ContaminationLen = 0;
	      Group[i].ContaminationMass = 0;
	    }
	}

      qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);

      MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

      for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	{
	  Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	  nimport += Recv_count[j];

	  if(j > 0)
	    {
	      Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
	      Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
	    }
	}

      ContamGet = (struct contamdata_in *) mymalloc(nimport * sizeof(struct contamdata_in));
      ContamIn = (struct contamdata_in *) mymalloc(nexport * sizeof(struct contamdata_in));

      /* prepare particle data for export */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  ContamIn[j].Pos[0] = Group[place].Pos[0];
	  ContamIn[j].Pos[1] = Group[place].Pos[1];
	  ContamIn[j].Pos[2] = Group[place].Pos[2];
	  ContamIn[j].R200 = Group[place].R_Mean200;

	  memcpy(ContamIn[j].NodeList,
		 DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	}

      /* exchange data */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;

	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* get the data */
		  MPI_Sendrecv(&ContamIn[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct contamdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A,
			       &ContamGet[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct contamdata_in), MPI_BYTE,
			       recvTask, TAG_DENS_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      myfree(ContamIn);
      ContamResult = (struct contamdata_out *) mymalloc(nimport * sizeof(struct contamdata_out));
      ContamOut = (struct contamdata_out *) mymalloc(nexport * sizeof(struct contamdata_out));

      /* now do the locations that were sent to us */
      for(j = 0; j < nimport; j++)
	subfind_contamination_evaluate(j, 1, &dummy, &dummy);

      if(i >= Ngroups)
	ndone_flag = 1;
      else
	ndone_flag = 0;

      MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

      /* get the result */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;
	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* send the results */
		  MPI_Sendrecv(&ContamResult[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct contamdata_out),
			       MPI_BYTE, recvTask, TAG_DENS_B,
			       &ContamOut[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct contamdata_out),
			       MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      /* add the result to the local particles */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;
	  Group[place].ContaminationLen += ContamOut[j].ContaminationLen;
	  Group[place].ContaminationMass += ContamOut[j].ContaminationMass;
	}

      myfree(ContamOut);
      myfree(ContamResult);
      myfree(ContamGet);
    }
  while(ndone < NTask);

  myfree(DataNodeList);
  myfree(DataIndexTable);
}

/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
int subfind_contamination_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int startnode, listindex = 0;
  double h, mass, masssum;
  int count, countsum;
  MyDouble *pos;

  masssum = 0;
  countsum = 0;

  if(mode == 0)
    {
      pos = Group[target].Pos;
      h = Group[target].R_Mean200;
    }
  else
    {
      pos = ContamGet[target].Pos;
      h = ContamGet[target].R200;
    }

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = ContamGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  count =
	    subfind_contamination_treefind(pos, h, target, &startnode, mode, nexport, nsend_local, &mass);

	  if(count < 0)
	    return -1;

	  masssum += mass;
	  countsum += count;
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = ContamGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
      Group[target].ContaminationMass = masssum;
      Group[target].ContaminationLen = countsum;
    }
  else
    {
      ContamResult[target].ContaminationMass = masssum;
      ContamResult[target].ContaminationLen = countsum;
    }

  return 0;
}

int subfind_contamination_treefind(MyDouble searchcenter[3], MyFloat hsml, int target, int *startnode,
				   int mode, int *nexport, int *nsend_local, double *Mass)
{
  int no, p, task, nexport_save;
  struct NODE *current;
  double mass;
  int count;
  MyDouble dx, dy, dz, dist, r2;

#ifdef PERIODIC
  MyDouble xtmp;
#endif
  nexport_save = *nexport;

  mass = 0;
  count = 0;

  no = *startnode;

  while(no >= 0)
    {
      if(no < All.MaxPart)	/* single particle */
	{
	  p = no;
	  no = Nextnode[no];

	  dist = hsml;
	  dx = NGB_PERIODIC_LONG_X(P[p].Pos[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(P[p].Pos[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(P[p].Pos[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  if(dx * dx + dy * dy + dz * dz > dist * dist)
	    continue;

	  mass += P[p].Mass;
	  count++;
	}
      else
	{
	  if(no >= All.MaxPart + MaxNodes)	/* pseudo particle */
	    {
	      if(mode == 1)
		endrun(12312);

	      if(mode == 0)
		{
		  if(Exportflag[task = DomainTask[no - (All.MaxPart + MaxNodes)]] != target)
		    {
		      Exportflag[task] = target;
		      Exportnodecount[task] = NODELISTLENGTH;
		    }

		  if(Exportnodecount[task] == NODELISTLENGTH)
		    {
		      if(*nexport >= All.BunchSize)
			{
			  *nexport = nexport_save;
			  if(nexport_save == 0)
			    endrun(13005);	/* in this case, the buffer is too small to process even a single particle */
			  for(task = 0; task < NTask; task++)
			    nsend_local[task] = 0;
			  for(no = 0; no < nexport_save; no++)
			    nsend_local[DataIndexTable[no].Task]++;
			  return -1;
			}
		      Exportnodecount[task] = 0;
		      Exportindex[task] = *nexport;
		      DataIndexTable[*nexport].Task = task;
		      DataIndexTable[*nexport].Index = target;
		      DataIndexTable[*nexport].IndexGet = *nexport;
		      *nexport = *nexport + 1;
		      nsend_local[task]++;
		    }

		  DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]++] =
		    DomainNodeIndex[no - (All.MaxPart + MaxNodes)];

		  if(Exportnodecount[task] < NODELISTLENGTH)
		    DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]] = -1;
		}

	      no = Nextnode[no - MaxNodes];
	      continue;
	    }

	  current = &Nodes[no];

	  if(mode == 1)
	    {
	      if(current->u.d.bitflags & (1 << BITFLAG_TOPLEVEL))	/* we reached a top-level node again, which means that we are done with the branch */
		{
		  *startnode = -1;
		  *Mass = mass;
		  return count;
		}
	    }

	  no = current->u.d.sibling;	/* in case the node can be discarded */

	  dist = hsml + 0.5 * current->len;;
	  dx = NGB_PERIODIC_LONG_X(current->center[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(current->center[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(current->center[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  /* now test against the minimal sphere enclosing everything */
	  dist += FACT1 * current->len;
	  if((r2 = (dx * dx + dy * dy + dz * dz)) > dist * dist)
	    continue;

	  no = current->u.d.nextnode;	/* ok, we need to open the node */
	}
    }

  *startnode = -1;

  *Mass = mass;
  return count;
}

#endif

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#ifndef _GENMAP_QUALITY_H_
#define _GENMAP_QUALITY_H_

#include <gslib.h>

void printPartStat(long long *vtx, int nel, int nv, comm_ext ce);

#endif

ext: h

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>  /* GSL Integration Library */
#include "funs.c"                 /* Integration functions   */
#include "timing_funs.h"          /* Timing functions        */

#define NMAX 100

int main(void)
{
	/* Declaration and initialization */
	int i;
	double vals1[NMAX+1], vals2[NMAX+1];

	/*Uses functions declared in "funs.c" to compute I(n) using the recurrence and the 		generic integrator methods, and stores the results in their respective arrays.*/
	integral_recur(0, 100, vals1);
	integral_gen(0, 100, vals2);

	/*Prints "n" and the result I(n) for each method above.  The absolute error 		multiplied by a factor of 1.e9 is reported in the last column.*/

	printf("n    Recurrence    Generic       Scaled Abs Err\n");

	for(i = 0; i <=NMAX; i++)
	{
		printf("%d    %.8f    %.8f    %f\n", i, vals1[i], vals2[i],
		(fabs(vals1[i]-vals2[i]))*1.e9);
	}

	/* Define the intial number of counts, tmin, tmax*/
	double time, trr, tgsl;
	int tmin = 1;
	int tmax = 2;
	int count = 1000;

	/* Computes the time per function call of the recurrence integration method.  Calls 		the function integral_recur() "count" times and adjusts count until the total time 		is between 1 and 2 seconds.*/
	printf("\n Timing for Recurrence Relation:\n");
	do
	{
		timer_start ();
		for(i = 0; i < count; i++)
		{
			integral_recur(0, 100, vals1);
		}
		time = timer_stop ();
		trr = time /count;
		printf (" %10.2f usec     %10.6f sec    %10d\n",trr * 1.e6, time, count);
		count = adjust_rep_count (count, time, tmin, tmax);
	}
	while ((time >tmax) || (time <tmin));

	/* Computes the time per function call of the GSL integrator method.  Calls 		the function integral_gen() "count" times and adjusts count until the total time is 		between 1 and 2 seconds. */
	printf("\n Timing for GSL Integrator:\n");
	do
	{
		timer_start ();
		for(i = 0; i < count; i++)
		{
			integral_gen(0, 100, vals1);
		}
		time = timer_stop ();
		tgsl = time /count;
		printf (" %10.2f usec     %10.6f sec    %10d\n",tgsl * 1.e6, time, count);
		count = adjust_rep_count (count, time, tmin, tmax);
	}
	while ((time >tmax) || (time <tmin));

	/*Prints time per function call for leibniz and bbp series, and then takes
	the ratio between the times to compare the relative speed*/
	printf("\nTime per function call for Recur:      %10.6f usec", trr*1.e6);
	printf("\nTime per function call for GSL:        %10.6f usec", tgsl*1.e6);
	int ratio = round(tgsl/trr);
	printf("\nRatio of function times tgsl/trr:        %d\n",ratio);

}


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

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

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

/* THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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. */

#ifndef UKF_PARAMETER_SCALAR_H
#define UKF_PARAMETER_SCALAR_H

#include <gsl/gsl_linalg.h> // For the Cholesky decomposition
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>

#include "ukf_types.h"
#include "ukf_math.h"

namespace ukf
{
    /**
      * @short UKF for parameter estimation.
      * The notations follow "Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models",p93, PhD, van Der Merwe
      */
    namespace parameter
    {
        /**
         * @short Allocation of the vectors/matrices and initialization
         *
         */
        inline void ukf_scalar_init(ukf_param &p, ukf_scalar_state &s)
        {
            // Init the lambda
            p.lambda = p.alpha * p.alpha * (p.n + p.kpa) - p.n;
            p.gamma = sqrt(p.n + p.lambda);
            p.nbSamples = 2 * p.n + 1;

            // Init the matrices used to iterate
            s.Kk = gsl_vector_alloc(p.n); // Kalman gain
            gsl_vector_set_zero(s.Kk);

            s.Kk_mat = gsl_matrix_alloc(p.n,1);
            gsl_matrix_set_zero(s.Kk_mat);

            s.Kk_mat_T = gsl_matrix_alloc(1,p.n);
            gsl_matrix_set_zero(s.Kk_mat_T);

            s.Pwdk = gsl_vector_alloc(p.n); // Covariance of the parameters and output
            gsl_vector_set_zero(s.Pwdk);

            // Whatever the type of evolution noise, its covariance is set to evolution_noise
            s.Prrk = gsl_matrix_alloc(p.n,p.n);
            p.evolution_noise->init(p,s);

            s.Peek = p.observation_noise; // Covariance of the observation noise

            s.Pddk = 0.0; // Covariance of the output

            s.w = gsl_vector_alloc(p.n); // Parameter vector
            gsl_vector_set_zero(s.w);

            s.wk = gsl_vector_alloc(p.n); // Vector holding one sigma point
            gsl_vector_set_zero(s.wk);

            s.Pk = gsl_matrix_alloc(p.n,p.n); // Covariance matrix
            gsl_matrix_set_identity(s.Pk);
            gsl_matrix_scale(s.Pk,p.prior_pi);

            s.Sk = gsl_matrix_alloc(p.n,p.n); // Matrix holding the cholesky decomposition of Pk
            // Initialize Sk to the cholesky decomposition of Pk
            gsl_matrix_memcpy(s.Sk, s.Pk);
            gsl_linalg_cholesky_decomp(s.Sk);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; k++)
                for(int j = 0 ; j < k ; j++)
                    gsl_matrix_set(s.Sk,j,k,0.0);

            s.cSk = gsl_vector_alloc(p.n); // Vector holding one column of Lpi
            gsl_vector_set_zero(s.cSk);

            s.wm  = gsl_vector_alloc(p.nbSamples); // Weights used to compute the mean of the sigma points images
            s.wc = gsl_vector_alloc(p.nbSamples); // Weights used to update the covariance matrices

            // Set the weights
            gsl_vector_set(s.wm, 0, p.lambda / (p.n + p.lambda));
            gsl_vector_set(s.wc, 0, p.lambda / (p.n + p.lambda) + (1.0 - p.alpha*p.alpha + p.beta));
            for(int j = 1 ; j < p.nbSamples; j ++)
            {
                gsl_vector_set(s.wm, j, 1.0 / (2.0 * (p.n + p.lambda)));
                gsl_vector_set(s.wc, j, 1.0 / (2.0 * (p.n + p.lambda)));
            }

            s.dk = gsl_vector_alloc(p.nbSamples); // Holds the image of the sigma points
            gsl_vector_set_zero(s.dk);

            s.d_mean = 0; // Holds the mean of the sigma points images

            s.sigmaPoints = gsl_matrix_alloc(p.n,p.nbSamples); // Holds the sigma points in the columns
            gsl_matrix_set_zero(s.sigmaPoints);

            s.temp_n = gsl_vector_alloc(p.n);
            gsl_vector_set_zero(s.temp_n);

            s.temp_n_n = gsl_matrix_alloc(p.n,p.n);
            gsl_matrix_set_zero(s.temp_n_n);
        }

        /**
         * @short Free of memory allocation
         *
         */
        inline void ukf_scalar_free(ukf_param &p, ukf_scalar_state &s)
        {
            gsl_vector_free(s.Kk);
            gsl_matrix_free(s.Kk_mat);
            gsl_matrix_free(s.Kk_mat_T);
            gsl_vector_free(s.Pwdk);
            gsl_matrix_free(s.Prrk);

            gsl_vector_free(s.w);
            gsl_vector_free(s.wk);

            gsl_matrix_free(s.Pk);
            gsl_matrix_free(s.Sk);
            gsl_vector_free(s.cSk);

            gsl_vector_free(s.wm);
            gsl_vector_free(s.wc);

            gsl_vector_free(s.dk);

            gsl_matrix_free(s.sigmaPoints);

            gsl_vector_free(s.temp_n);
            gsl_matrix_free(s.temp_n_n);
        }

        /**
         * @short Iteration for UKF for parameter estimation, in case of a scalar output
         *
         */
	template<typename GFUNC>
        void ukf_scalar_iterate(ukf_param &p, ukf_scalar_state &s,
				       GFUNC g,
				       gsl_vector * xk, double dk)
        {
            // Here, we implement the UKF for parameter estimation in the scalar case
            // The notations follow p93 of the PhD thesis of Van Der Merwe, "Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models"

            // ************************************************** //
            // ************ Time update equations    ************ //
            // ************************************************** //
            // Add the evolution noise to the parameter covariance Eq 3.137
            gsl_matrix_add(s.Pk, s.Prrk);

            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //
            // Equations 3.138
            // w_k^j = w_(k-1)  <-- this is here denoted s.w
            // w_k^j = w_(k-1) + gamma Sk_j for 1 <= j <= n
            // w_k^j = w_(k-1) - gamma Sk_j for n+1 <= j <= 2n

            // Perform a cholesky decomposition of Pk
            gsl_matrix_memcpy(s.Sk, s.Pk);
            gsl_linalg_cholesky_decomp(s.Sk);
            // Set all the elements of Lpi strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; k++)
                for(int j = 0 ; j < k ; j++)
                    gsl_matrix_set(s.Sk,j,k,0.0);

            gsl_matrix_set_col(s.sigmaPoints,0, s.w);
            for(int j = 1 ; j < p.n+1 ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) + p.gamma * gsl_matrix_get(s.Sk,i,j-1));

            for(int j = p.n+1 ; j < p.nbSamples ; ++j)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) - p.gamma * gsl_matrix_get(s.Sk,i,j-(p.n+1)));

            // ************************************************** //
            // ***** Compute the images of the sigma points ***** //
            // ************************************************** //

            // Compute the images of the sigma points
            // and the mean of the yj
            s.d_mean = 0.0;
            gsl_vector_set_zero(s.dk);
            for(int j = 0 ; j < p.nbSamples ; j++)
            {
                // Equation 3.139
                gsl_matrix_get_col(s.wk, s.sigmaPoints,j);
                s.dk->data[j] = g(s.wk,xk);

                // Update the mean : y_mean = sum_[j=0..2n] wm_j y_j
                // Equation 3.140
                s.d_mean += s.wm->data[j] * s.dk->data[j];
            }

            // ************************************************** //
            // ************** Update the statistics ************* //
            // ************************************************** //

            gsl_vector_set_zero(s.Pwdk);
            // The covariance of the output is initialized with the observation noise covariance
            // Eq 3.142
            s.Pddk = s.Peek;
            for(int j = 0 ; j < p.nbSamples ; j++)
            {
                // Eq 3.142
                s.Pddk += s.wc->data[j] * gsl_pow_2(s.dk->data[j] - s.d_mean);
                // Eq 3.143
                for(int i = 0 ; i < p.n ; ++i)
                {
                    s.Pwdk->data[i] = s.Pwdk->data[i] + s.wc->data[j] * (gsl_matrix_get(s.sigmaPoints,i,j) - s.w->data[i]) * (s.dk->data[j] - s.d_mean) ;
                }
            }

            // ************************************************** //
            // ******* Kalman gain and parameters update ******** //
            // ************************************************** //

            //if(s.Pddk == 0.0)
            //    printf("[Error] Output covariance is null !");
            // May not occur as soon as the observation noise covariance is set != 0.0

            // Eq. 3.144
            for(int i = 0 ; i < p.n ; ++i)
                s.Kk->data[i] = s.Pwdk->data[i] / s.Pddk;

            // Eq 3.145
            // wk = w_(k-1) + Kk * (dk - d_mean)
            s.ino_dk = dk - s.d_mean;
            for(int i = 0 ; i < p.n ; ++i)
                s.w->data[i] = s.w->data[i] + s.Kk->data[i] * s.ino_dk;

            // Eq. 3.146
            // Pk = P_(k-1) - Pddk . Kk Kk^T
            for(int i = 0 ; i < p.n ; ++i)
            {
                gsl_matrix_set(s.Kk_mat, i, 0, s.Kk->data[i]);
                gsl_matrix_set(s.Kk_mat_T, 0, i, s.Kk->data[i]);
            }
            gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,-s.Pddk, s.Kk_mat, s.Kk_mat_T,1.0,s.Pk);

            // Update of the evolution noise
            p.evolution_noise->updateEvolutionNoise(p,s);
        }

        /**
        * @short Evaluation of the output from the sigma points
        *
        */
        void ukf_scalar_evaluate(ukf_param &p, ukf_scalar_state &s,
				 double(*g)(gsl_vector*, gsl_vector*),
				 gsl_vector * xk, double &dk)
        {
            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //

            // 1- Compute the cholesky decomposition of Pk
            gsl_matrix_memcpy(s.temp_n_n, s.Pk);
            gsl_linalg_cholesky_decomp(s.temp_n_n);

            // 2 - Set all the elements of Sk_temp strictly above the diagonal to zero
            for(int k = 0 ; k < p.n ; k++)
                for(int j = 0 ; j < k ; j++)
                    gsl_matrix_set(s.temp_n_n,j,k,0.0);
            // Now Sk_temp is a lower triangular matrix containing the cholesky decomposition of Pk

            // 3- Compute the sigma points
            // Equations 3.138
            // w_k^j = w_(k-1)  <-- this is here denoted s.w
            // w_k^j = w_(k-1) + gamma Sk_j for 1 <= j <= n
            // w_k^j = w_(k-1) - gamma Sk_j for n+1 <= j <= 2n
            gsl_matrix_set_col(s.sigmaPoints,0, s.w);
            for(int j = 1 ; j < p.n+1 ; j++)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) + p.gamma * gsl_matrix_get(s.temp_n_n,i,j-1));

            for(int j = p.n+1 ; j < p.nbSamples ; j++)
                for(int i = 0 ; i < p.n ; ++i)
                    gsl_matrix_set(s.sigmaPoints, i, j, gsl_vector_get(s.w, i) - p.gamma * gsl_matrix_get(s.temp_n_n,i,j-(p.n+1)));

            // ************************************************** //
            // ***** Compute the images of the sigma points ***** //
            // ************************************************** //

            // Compute the images of the sigma points
            // and their mean
            dk = 0.0;
            for(int j = 0 ; j < p.nbSamples ; j++)
            {
                gsl_matrix_get_col(s.wk, s.sigmaPoints,j);
                // Update the mean : d_mean = sum_[j=0..2n] w_j y_j
                dk += gsl_vector_get(s.wm,j) * g(s.wk,xk);
            }
        }

    } // parameter
} // ukf

#endif // UKF_PARAMETER_SCALAR_H


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

#pragma once

#include <gsl/gsl-lite.hpp>
#include <units/bits/constexpr_math.h>
#include <units/bits/math_concepts.h>
#include <units/bits/pow.h>
#include <units/bits/ratio_maths.h>
#include <cmath>

namespace units::detail {

template<std::intmax_t N, typename F>
  requires gt_zero<N>
[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept
{
  if constexpr (N == 1) {
    return v;
  } else {
    gsl_Expects(v >= 0);
    if (v == 0) {
      return 0;
    }

    constexpr double exponent = 1.0 / static_cast<double>(N);
    const auto root = static_cast<std::intmax_t>(pow_function(static_cast<double>(v), exponent));

    // integer roots may be truncated down by 1 or in even rarer cases up by 1 due to finite precision of pow and
    // exponent, check both cases
    if (v == pow_impl<N>(root + 1)) {
      return root + 1;
    }
    if (v < pow_impl<N>(root)) {
      return root - 1;
    }
    return root;
  }
}

// maximum v is std::numeric_limits<std::intmax_t>::max() which is the worst case for exp convergence
// ExpOrder = 12 and Factor = 64 give a precision of about O(1e-15) for a wide range of 1 / N exponents
// Factor = 32 needs quite a few more terms to converge
// https://godbolt.org/z/odWq1o
template<std::intmax_t N, std::size_t ExpOrder = 12, std::intmax_t Factor = 64>
  requires gt_zero<N>
[[nodiscard]] constexpr std::intmax_t iroot_compile(std::intmax_t v) noexcept
{
  return iroot_impl<N>(v, [](double x, double exponent) { return constexpr_pow<ExpOrder, Factor>(x, exponent); });
}

template<std::intmax_t N>
  requires gt_zero<N>
[[nodiscard]] std::intmax_t iroot_runtime(std::intmax_t v) noexcept
{
  return iroot_impl<N>(v, [](double x, [[maybe_unused]] double exponent) {
    if constexpr (N == 2) {
      return std::sqrt(x);
    } else if constexpr (N == 3) {
      return std::cbrt(x);
    } else {
      return std::pow(x, exponent);
    }
  });
}

template<std::intmax_t N>
  requires gt_zero<N>
[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept
{
  // compile time version is much slower, use faster version at runtime
  if (std::is_constant_evaluated()) {
    return iroot_compile<N>(v);
  }

  return iroot_runtime<N>(v);
}

}  // namespace units::detail

max_stars_repo_name: HazardyKnusperkeks/units

max_issues_repo_name: HazardyKnusperkeks/units

max_forks_repo_name: HazardyKnusperkeks/units


text:
/*
   LAPACKE Example : Calling DGELS using col-major layout
   =====================================================

   The program computes the solution to the system of linear
   equations with a square matrix A and multiple
   right-hand sides B, where A is the coefficient matrix
   and b is the right-hand side matrix:

   Description
   ===========

   In this example, we wish solve the least squares problem min_x || B - Ax ||
   for two right-hand sides using the LAPACK routine DGELS. For input we will
   use the 5-by-3 matrix

         ( 1  1  1 )
         ( 2  3  4 )
     A = ( 3  5  2 )
         ( 4  2  5 )
         ( 5  4  3 )
    and the 5-by-2 matrix

         ( -10 -3 )
         (  12 14 )
     B = (  14 12 )
         (  16 16 )
         (  18 16 )
    We will first store the input matrix as a static C two-dimensional array,
    which is stored in col-major layout, and let LAPACKE handle the work space
    array allocation. The LAPACK base name for this function is gels, and we
    will use double precision (d), so the LAPACKE function name is LAPACKE_dgels.

    lda=5 and ldb=5. The output for each right hand side is stored in b as
    consecutive vectors of length 3. The correct answer for this problem is
    the 3-by-2 matrix

         ( 2 1 )
         ( 1 1 )
         ( 1 2 )

    A complete C program for this example is given below. Note that when the arrays
     are passed to the LAPACK routine, they must be dereferenced, since LAPACK is
      expecting arrays of type double *, not double **.

   LAPACKE Interface
   =================

   LAPACKE_dgels (col-major, high-level) Example Program Results

  -- LAPACKE Example routine (version 3.7.0) --
  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
     December 2016

*/
/* Calling DGELS using col-major layout */

/* Includes */
#include <stdio.h>
#include <lapacke.h>

 int main() {
   double A[5][3] = {{1,2,3},{4,5,1},{3,5,2},{4,1,4},{2,5,3}};
   double b[5][2] = {{-10,12},{14,16},{18,-3},{14,12},{16,16}};

   const lapack_int m = 5;
   const lapack_int n = 3;
   const lapack_int lda = 5;
   const lapack_int ldb = 5;
   const lapack_int nrhs = 2;

   printf("LAPACKE_dgels (col-major, high-level) Example Program Results\n" );
   const lapack_int info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m, n, nrhs, *A, lda, *b, ldb);
   printf( "Solution: %d %d %lf %d\n", n, nrhs, b[0][0], ldb);

   return info;
}

max_stars_repo_name: nabbelbabbel/conan-lapack

max_issues_repo_name: nabbelbabbel/conan-lapack

max_forks_repo_name: nabbelbabbel/conan-lapack


text:
#pragma once

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

namespace Library
{
	class RasterizerStates final
	{
	public:
		inline static winrt::com_ptr<ID3D11RasterizerState> BackCulling;
		inline static winrt::com_ptr<ID3D11RasterizerState> FrontCulling;
		inline static winrt::com_ptr<ID3D11RasterizerState> DisabledCulling;
		inline static winrt::com_ptr<ID3D11RasterizerState> Wireframe;

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

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


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <Python.h>

static char module_docstring[] = "GSL Statistics module";
static char t_test_docstring[] = "Perform an independent samples t-test";

double t_test(const double *data1, int n1, const double *data2, int n2)
{
    double mean1, mean2, var1, var2, num, sp, denom, t, dof, p;
    mean1 = gsl_stats_mean(data1, 1, n1);
    mean2 = gsl_stats_mean(data2, 1, n2);
    var1 = gsl_stats_variance(data1, 1, n1);
    var2 = gsl_stats_variance(data2, 1, n2);
    num = mean1 - mean2;
    sp = sqrt((((n1 - 1) * var1) + ((n2 - 1) * var2)) / (n1 + n2 - 2));
    denom = sp * sqrt((1 / (double) n1) + (1 / (double) n2));
    t = fabs(num / denom);
    dof = (double) n1 + n2 - 2;
    p = (1 - gsl_cdf_tdist_P(t, dof)) * 2;
    return p;
}

static double* PySequence_ToDoubleArray(PyObject* inp, int n)
{
    // Allocate array
    double *out;
    out = malloc(n * sizeof(double));
    for (int i = 0; i < n; i++)
    {
        PyObject *fitem;
        PyObject *item = PySequence_Fast_GET_ITEM(inp, i);
        if(!item) {
            Py_DECREF(inp);
            free(out);
            return NULL;
        }
        fitem = PyNumber_Float(item);
        if(!fitem) {
            Py_DECREF(inp);
            free(out);
            PyErr_SetString(PyExc_TypeError, "all items must be numbers");
            return NULL;
        }
        out[i] = PyFloat_AS_DOUBLE(fitem);
        Py_DECREF(fitem);
    }
    return out;
}

static PyObject *t_test_py(PyObject *self, PyObject *args)
{
    PyObject *inp1, *inp2;
    PyObject* data1;
    PyObject* data2;
    double *d1, *d2;
    int n1, n2;
    double res;

    if (!PyArg_ParseTuple(args, "OO", &inp1, &inp2))
        return 0;
    data1 = PySequence_Fast(inp1, "argument must be iterable");
    data2 = PySequence_Fast(inp2, "argument must be iterable");
    n1 = PySequence_Fast_GET_SIZE(data1);
    n2 = PySequence_Fast_GET_SIZE(data2);

    d1 = PySequence_ToDoubleArray(data1, n1);
    d2 = PySequence_ToDoubleArray(data2, n2);

    Py_DECREF(data1);
    Py_DECREF(data2);

    res = t_test(d1, n1, d2, n2);
    free(d1);
    free(d2);

    return Py_BuildValue("d", res);
}

static PyMethodDef module_methods[] = {
    {"t_test_py", t_test_py, METH_VARARGS, t_test_docstring},
    {NULL, NULL, 0, NULL}
};

static struct PyModuleDef gslstatsmodule = {
    PyModuleDef_HEAD_INIT,
    "gslstats",
    module_docstring,
    -1,
    module_methods,
};

PyMODINIT_FUNC PyInit_gslstats(void)
{
    Py_Initialize();
    return PyModule_Create(&gslstatsmodule);
}

max_stars_repo_name: kgoettler/python-c-ext

max_issues_repo_name: kgoettler/python-c-ext

max_forks_repo_name: kgoettler/python-c-ext


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <mpi.h>
#include "compearth.h"
#include "parmt_mtsearch.h"
#ifdef PARMT_USE_INTEL
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

/*
int main_temp()
{
    double *betas, *gammas, *kappas, *M0s, *mts, *sigmas, *thetas;
    int i, ierr, nmt;
    const int ng = 29;
    const int nb = 29;
    const int nk = 6;
    const int ns = 6;
    const int nt = 6;
    const int nm = 1;
    const int ldm = 8;
    const double oneDeg = M_PI/180.0;
    const double betaMin = 0.0 + oneDeg;
    const double betaMax = M_PI - oneDeg;
    const double gammaMin =-M_PI/6.0 + oneDeg;
    const double gammaMax = M_PI/6.0 - oneDeg;
    const double kappaMin = 0.0 + oneDeg;
    const double kappaMax = 2.0*M_PI - oneDeg;
    const double thetaMin = 0.0 + oneDeg;
    const double thetaMax = 0.5*M_PI - oneDeg;
    const double sigmaMin =-0.5*M_PI + oneDeg;
    const double sigmaMax = 0.5*M_PI - oneDeg;
    nmt = nm*nb*ng*nk*ns*nt;
    mts = memory_calloc64f(nmt*ldm); //(double *)calloc((size_t) (nmt*ldm), sizeof(double));
    M0s = memory_calloc64f(nm); //(double *)calloc((size_t) nm, sizeof(double));
    betas = memory_calloc64f(nb); //(double *)calloc((size_t) nb, sizeof(double));
    gammas = memory_calloc64f(ng); //(double *)calloc((size_t) ng, sizeof(double));
    kappas = memory_calloc64f(nk); //(double *)calloc((size_t) nk, sizeof(double));
    sigmas = memory_calloc64f(ns); //(double *)calloc((size_t) ns, sizeof(double));
    thetas = memory_calloc64f(nt); //(double *)calloc((size_t) nt, sizeof(double));

    ierr = array_linspace64f_work(1.0/sqrt(2.0), 1.0/sqrt(2.0), nm, M0s);
    ierr = array_linspace64f_work(betaMin+0.001,  betaMax-0.001,  nb, betas);
    ierr = array_linspace64f_work(gammaMin, gammaMax, ng, gammas);
    ierr = array_linspace64f_work(kappaMin, kappaMax, nk, kappas);
    ierr = array_linspace64f_work(thetaMin, thetaMax, nt, thetas);
    ierr = array_linspace64f_work(sigmaMin, sigmaMax, ns, sigmas);
    // avoid warnings on theta in tt2cmt
    ierr = gridSearch_discretizeMT64f(ng, gammas, nb, betas,  nm, M0s,
                                      nk, kappas, nt, thetas, ns, sigmas,
                                      ldm, nmt, mts);

double Mdc[6], lam[3], U[6], theta;
    compearth_tt2cmt(0.0, 0.0, 1./sqrt(2.),
                     0.0, 0.1, 0.0, Mdc, lam, U);

    FILE *ofl = fopen("mts.txt", "w");
    for (i=0; i<nmt; i++)
    {
        compearth_angleMT(1, Mdc, &mts[ldm*i], &theta);

        fprintf(ofl, "%f %f %f %f %f %f %f\n", //i+1,
                mts[ldm*i+0], mts[ldm*i+1], mts[ldm*i+2],
                mts[ldm*i+3], mts[ldm*i+4], mts[ldm*i+5], theta);
    }
    fclose(ofl);

    free(mts);
    free(M0s);
    free(betas);
    free(gammas);
    free(kappas);
    free(sigmas);
    free(thetas);
    return 0;
}
*/
//============================================================================//
/*!
 * @brief Provides a cell based discretization of the moment tensor space.
 *
 * @param[in] nb          Number of colatitudes.
 * @param[in] betaLower   Lower colatitude in grid search [0, betaUpper).
 * @param[in] betaUpper   Upper colatitude in grid search (betaLower, pi].
 * @param[in] ng          Number of longitudes.
 * @param[in] gammaLower  Lower longitude in grid search [-pi/6, gammaUpper).
 * @param[in] gammaUpper  Upper longitude in grid search (gammaLower, pi/6].
 * @param[in] nk          Number of strikes in grid search.
 * @param[in] kappaLower  Lower strike in grid search [0, kappaUpper).
 * @param[in] kappaUpper  Upper strike in grid search (kappaLower, 2*pi].
 * @param[in] ns          Number of slips in grid search.
 * @param[in] sigmaLower  Lower slip angle in grid search [-pi, sigmaUpper).
 * @param[in] sigmaUpper  Upper slip angle in grid search (sigmaUpper, pi].
 * @param[in] nt          Number of dips in grid search.
 * @param[in] thetaLower  Lower dip in grid search [0,dipUpper).
 * @param[in] thetaUpper  Upper dip in grid search (dipLower, pi/2].
 * @param[in] nm          Number of magnitudes.
 * @param[in] m0Lower     Lower scalar moment in grid search (Newton-meters).
 * @param[in] m0Upper     Upper scalar moment in grid search (Newton-meters).
 * @param[in] luseLog     If true then the the scalar moments will be
 *                        discretized on a log scale.
 *
 * @param[out] betas      Colatitudes (radians) at cell centers [nb].
 * @param[out] gammas     Longitudes (radians) at cell centers [ng].
 * @param[out] kappas     Strike angles (radians) at cell centers [nk].
 * @param[out] sigmas     Slip angles (radians) at cell centers [ns].
 * @param[out] thetas     Dip angles (radians) at cell centers [nt].
 * @param[out] M0s        Scalar moments (Newton-meters) in grid search [nm].
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_discretizeCells64f(
    const int nb, const double betaLower, const double betaUpper,
    const int ng, const double gammaLower, const double gammaUpper,
    const int nk, const double kappaLower, const double kappaUpper,
    const int ns, const double sigmaLower, const double sigmaUpper,
    const int nt, const double thetaLower, const double thetaUpper,
    const int nm, const double m0Lower, const double m0Upper,
    const bool luseLog,
    double **betas,  double **gammas, double **kappas,
    double **sigmas, double **thetas, double **M0s)
{
    const char *fcnm = "parmt_discretizeCells64f\0";
    double *mwl, *u, *v, *h;
    double du, du2, dv, dv2, dh, dh2, dk, dk2, ds, ds2,
           hLower, hUpper, mwll, mwul, uLower, uUpper, vLower, vUpper;
    int ierr;
    const double two = 2.0;
    // Initialize and do some basic error checks
    u = NULL;
    v = NULL;
    h = NULL;
    if (nb < 1 || ng < 1 || nk < 1 || ns < 1 || nt < 1 || nm < 1)
    {
        if (nb < 1){printf("%s: no betas\n", fcnm);}
        if (ng < 1){printf("%s: no gammas\n", fcnm);}
        if (nk < 1){printf("%s: no kappas\n", fcnm);}
        if (ns < 1){printf("%s: no sigmas\n", fcnm);}
        if (nt < 1){printf("%s: no thetas\n", fcnm);}
        if (nm < 1){printf("%s: no magnitudes\n", fcnm);}
        return -1;
    }
    // Discretize the moment tensor space in (u, v, h) space
    compearth_beta2u(1, &betaLower, &uLower);
    compearth_beta2u(1, &betaUpper, &uUpper);
    compearth_gamma2v(1, &gammaLower, &vLower);
    compearth_gamma2v(1, &gammaUpper, &vUpper);
    compearth_theta2h(1, &thetaLower, &hLower);
    compearth_theta2h(1, &thetaUpper, &hUpper);
    du = (uUpper - uLower)/(double) nb;
    dv = (vUpper - vLower)/(double) ng;
    dh =-(hUpper - hLower)/(double) nt; // negative sign makes +theta increasing
    dk = (kappaUpper - kappaLower)/(double) nk;
    ds = (sigmaUpper - sigmaLower)/(double) ns;
    du2 = du/two;
    dv2 = dv/two;
    dh2 = dh/two;
    u = array_linspace64f(uLower+du2, uUpper-du2, nb, &ierr);
    v = array_linspace64f(vLower+dv2, vUpper-dv2, ng, &ierr);
    h = array_linspace64f(hLower-dh2, hUpper+dh2, nt, &ierr);
    *betas  = memory_calloc64f(nb);
    *gammas = memory_calloc64f(ng);
    *thetas = memory_calloc64f(nt);
    //compearth_u2beta(nb, 20, 2, u, 1.e-6, *betas);
    compearth_u2beta(nb, u, *betas);
    compearth_v2gamma(ng, v, *gammas);
    compearth_h2theta(nt, h, *thetas);
    // discretize the magnitudes, kappa, and sigma
    dk2 = dk/two;
    ds2 = ds/two;
    if (!luseLog || nm == 1)
    {
        *M0s = array_linspace64f(m0Lower, m0Upper, nm, &ierr);
    }
    else
    {
        // TODO: I doubt this is correct because I don't understand how
        // lune volumes map to `spherical volumes'.  For simplicity I'll
        // use a shell but I think I should be looking at the `active
        // area' of the beta and gamma grid search.  In the interim I'll
        // use the volume of a spherical shell which is given by
        // \frac{4 \pi R^3}{3} where R is the radius (magnitude).
        // To get some form of uniform volume spacing I'll take the log
        // at the upper and lower limit and discretize evenly in logspace
        // then come back to M0.
        /*
        m0ll = 3.0*log(4.0*M_PI*m0Lower/3.0);
        m0ul = 3.0*log(4.0*M_PI*m0Upper/3.0);
        m0l = array_linspace64f(m0ll, m0ul, nm, &ierr);
        // Now solve for m0: 3/(4*pi)*exp(m0l/3)
        cblas_dscal(nm, 1.0/3.0, m0l, 1);         // m0l = m0l/3
        *M0s = array_exp64f(nm, m0l, &ierr);      // *M0 = exp(m0l/3)
        cblas_dscal(nm, 3.0/(4.0*M_PI), *M0s, 1); // *M0 = 3/(4*pi)exp(m0l/3)
        memory_free64f(&m0l);
        */
        // TODO: An alternative to doing math is just discretizing on a
        // the magnitude scale
        compearth_m02mw(1, CE_KANAMORI_1978, &m0Lower, &mwll);
        compearth_m02mw(1, CE_KANAMORI_1978, &m0Upper, &mwul);
        mwl = array_linspace64f(mwll, mwul, nm, &ierr);
        *M0s = memory_calloc64f(nm);
        ierr = compearth_mw2m0(nm, CE_KANAMORI_1978, mwl, *M0s);
        memory_free64f(&mwl);
    }
    *kappas = array_linspace64f(kappaLower+dk2, kappaUpper-dk2, nk, &ierr);
    *sigmas = array_linspace64f(sigmaLower+ds2, sigmaUpper-ds2, ns, &ierr);
    memory_free64f(&u);
    memory_free64f(&v);
    memory_free64f(&h);
    return ierr;
}

//============================================================================//
/*!
 * @brief Computes the moment tensors in the grid search.
 *        The grid search ordering is:
 *         Loop on magnitudes
 *          Loop on colatitudes
 *           Loop on longitudes
 *            Loop on strikes
 *             Loop on slips
 *              Loop on dips
 *
 * @param[in] comm    MPI communicator on which moment tensor will be split
 * @param[in] ng      number of longitudes
 * @param[in] gammas  longitudes (radians) on lune s.t.
 *                    \f$ \gamma \in [-\pi/6, \pi/6] \f$ [ng]
 * @param[in] nb      number of colatitudes
 * @param[in] betas   colatitudes (radians) on lune s.t.
 *                    \f$ \beta \in [0, \pi] \f$ [nb]
 * @param[in] nm      number of scalar moments
 * @param[in] M0s     scalar moments in Newton-meters [nm].
 *                    for a `unit' moment tensor choose M0 = 1/sqrt(2)
 * @param[in] nk      number of strike angles
 * @param[in] kappas  strike angles (radians) s.t.
 *                    \f$ \kappa \in [0, 2\pi] \f$ [nk]
 * @param[in] nt      number of dip angles
 * @param[in] thetas  dips angles (radians) s.t.
 *                    \f$ \theta \in [0, \pi/2] \f$ [nt]
 * @param[in] ns      number of slip angles
 * @param[in] sigmas  slip angles (radians) s.t.
 *                    \f$ \sigma \in [-\pi/2, \pi/2] \f$ [ns]
 * @param[in] ldm     leading dimension of mts.  must be >= 6.
 * @param[in] nmt     number of moment tensors (=nm*nb*ng*nk*ns*nt)
 *
 * @param[out] mts    moment tensors (Newton-meters) in north,east,down s.t.
 *                    \f$ \textbf{m}
 *                      = \f$ \{m_{xx}, m_{yy}, m_{zz},
 *                              m_{xy}, m_{xz}, m_{yz}\} \f$.
 *                    This has dimensions [nm*nb*ng*nk*ns*nt x ldm].  Observe
 *                    the ordering as defind above.
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI licensed under Apached 2
 *
 */
int parmt_discretizeMT64f_MPI(const MPI_Comm comm,
                              const int ng,
                              const double *__restrict__ gammas,
                              const int nb,
                              const double *__restrict__ betas,
                              const int nm,
                              const double *__restrict__ M0s,
                              const int nk,
                              const double *__restrict__ kappas,
                              const int nt,
                              const double *__restrict__ thetas,
                              const int ns,
                              const double *__restrict__ sigmas,
                              const int ldm,
                              struct localMT_struct *mts)
{
    const char *fcnm = "parmt_discretizeMT64f_MPI\0";
    double *mtWork;
    int64_t nmt64;
    int *displ, *nmtProc, *offset, *sendCounts;
    int dmt, i, ierr, imt1, imt2, myid, nmtall, nmt, nprocs;
    const int master = 0;
    //------------------------------------------------------------------------//
    ierr = 0;
    MPI_Comm_size(comm, &nprocs);
    MPI_Comm_rank(comm, &myid);
    memset(mts, 0, sizeof(struct localMT_struct));
    mts->comm = comm;
    // Check for integer overflow
    nmt64 = (int64_t) ng * (int64_t) nb * (int64_t) nm
           *(int64_t) nk * (int64_t) nt * (int64_t) ns;
    if (nmt64 > INT_MAX)
    {
        printf("%s: Error integer overflow - gridsearch too large\n", fcnm);
        return -1;
    }
    nmt = (int) nmt64;
    // Verify the inputs
    if (ldm < 6 || ng < 1 || nb < 1 || nm < 1 || nk < 1 || nt < 1 || ns < 1 ||
        gammas == NULL || betas == NULL || M0s == NULL || kappas == NULL ||
        thetas == NULL || sigmas == NULL)
    {
        if (ldm < 6){printf("%s: Invalid leading dimension\n", fcnm);}
        if (ng < 1){printf("%s: ng must be positive\n", fcnm);}
        if (nb < 1){printf("%s: nb must be positive\n", fcnm);}
        if (nm < 1){printf("%s: nm must be positive\n", fcnm);}
        if (nk < 1){printf("%s: nk must be positive\n", fcnm);}
        if (nt < 1){printf("%s: nt must be positive\n", fcnm);}
        if (ns < 1){printf("%s: ns must be positive\n", fcnm);}
        if (gammas == NULL){printf("%s: gammas is NULL\n", fcnm);}
        if (betas == NULL){printf("%s: betas is NULL\n", fcnm);}
        if (M0s == NULL){printf("%s: M0s is NULL\n", fcnm);}
        if (kappas == NULL){printf("%s: kappas is NULL\n", fcnm);}
        if (thetas == NULL){printf("%s: thetas is NULL\n", fcnm);}
        if (sigmas == NULL){printf("%s: sigmas is NULL\n", fcnm);}
        return -1;
    }
    // divide the moment tensor grid
    dmt = MAX(nmt/nprocs, 1);
    imt1 = myid*dmt;
    imt2 = (myid + 1)*dmt;
    if (myid == nprocs - 1){imt2 = nmt;}
    mts->nmt = imt2 - imt1;
    mts->ldm = ldm;
    MPI_Allreduce(&mts->nmt, &nmtall, 1, MPI_INTEGER, MPI_SUM, mts->comm);
    if (nmtall != nmt)
    {
        if (myid == master)
        {
            printf("%s: Failed to partition domain %d %d\n", fcnm, nmtall, nmt);
        }
        return -1;
    }
    // create the requisite information for an MPI_GATHER
    mts->nmtProc = memory_calloc32i(nprocs);
    mts->offset = memory_calloc32i(nprocs);
    nmtProc = memory_calloc32i(nprocs);
    offset = memory_calloc32i(nprocs);
    nmtProc[myid] = mts->nmt;
    offset[myid] = imt1;
    MPI_Allreduce(nmtProc, mts->nmtProc, nprocs, MPI_INTEGER,
                  MPI_SUM, mts->comm);
    MPI_Allreduce(offset, mts->offset,  nprocs, MPI_INTEGER,
                  MPI_SUM, mts->comm);
    mts->commSize = nprocs;
    mts->nmtAll = nmtall;
    mts->myid = myid;
    memory_free32i(&nmtProc);
    memory_free32i(&offset);
    // have master process compute all mts and scatter them
    // TODO this should be parallel
    if (mts->myid == master)
    {
        mtWork = memory_calloc64f(mts->ldm*mts->nmtAll);
        ierr = parmt_discretizeMT64f(ng, gammas,
                                     nb, betas,
                                     nm, M0s,
                                     nk, kappas,
                                     nt, thetas,
                                     ns, sigmas,
                                     mts->ldm, mts->nmtAll, mtWork);
    }
    else
    {
        mtWork = memory_calloc64f(1);
    }
    MPI_Bcast(&ierr, 1, MPI_INT, master, mts->comm);
    if (mts->commSize > 1)
    {
        // figure out who gets what part of the moment tensor space
        sendCounts = memory_calloc32i(mts->commSize);
        displ = memory_calloc32i(mts->commSize);
        for (i=0; i<mts->commSize; i++)
        {
            sendCounts[i] = mts->nmtProc[i]*mts->ldm;
            displ[i] = mts->ldm*mts->offset[i];
        }
        // set the memory and distribute it
        mts->mts = memory_calloc64f(mts->ldm*MAX(1,mts->nmt));
        MPI_Scatterv(mtWork, sendCounts, displ,
                     MPI_DOUBLE, mts->mts, mts->ldm*mts->nmt,
                     MPI_DOUBLE, master, mts->comm);
        memory_free64f(&mtWork);
    }
    else
    {
        mts->mts = mtWork;
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Computes the moment tensors in the grid search.
 *        The grid search ordering is:
 *         Loop on magnitudes
 *          Loop on colatitudes
 *           Loop on longitudes
 *            Loop on strikes
 *             Loop on slips
 *              Loop on dips
 *
 * @param[in] ng      number of longitudes
 * @param[in] gammas  longitudes (radians) on lune s.t.
 *                    \f$ \gamma \in [-\pi/6, \pi/6] \f$ [ng]
 * @param[in] nb      number of colatitudes
 * @param[in] betas   colatitudes (radians) on lune s.t.
 *                    \f$ \beta \in [0, \pi] \f$ [nb]
 * @param[in] nm      number of scalar moments
 * @param[in] M0s     scalar moments in Newton-meters [nm].
 *                    for a `unit' moment tensor choose M0 = 1/sqrt(2)
 * @param[in] nk      number of strike angles
 * @param[in] kappas  strike angles (radians) s.t.
 *                    \f$ \kappa \in [0, 2\pi] \f$ [nk]
 * @param[in] nt      number of dip angles
 * @param[in] thetas  dips angles (radians) s.t.
 *                    \f$ \theta \in [0, \pi/2] \f$ [nt]
 * @param[in] ns      number of slip angles
 * @param[in] sigmas  slip angles (radians) s.t.
 *                    \f$ \sigma \in [-\pi/2, \pi/2] \f$ [ns]
 * @param[in] ldm     leading dimension of mts.  must be >= 6.
 * @param[in] nmt     number of moment tensors (=nm*nb*ng*nk*ns*nt)
 *
 * @param[out] mts    moment tensors (Newton-meters) in north,east,down s.t.
 *                    \f$ \textbf{m}
 *                      = \f$ \{m_{xx}, m_{yy}, m_{zz},
 *                              m_{xy}, m_{xz}, m_{yz}\} \f$.
 *                    This has dimensions [nm*nb*ng*nk*ns*nt x ldm].  Observe
 *                    the ordering as defind above.
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI licensed under Apached 2
 *
 */
int parmt_discretizeMT64f(const int ng,
                          const double *__restrict__ gammas,
                          const int nb,
                          const double *__restrict__ betas,
                          const int nm,
                          const double *__restrict__ M0s,
                          const int nk,
                          const double *__restrict__ kappas,
                          const int nt,
                          const double *__restrict__ thetas,
                          const int ns,
                          const double *__restrict__ sigmas,
                          const int ldm,
                          const int nmt, double *__restrict__ mts)
{
    const char *fcnm = "parmt_discretizeMT64f\0";
    double lam[3], Muse[6], U[9] __attribute__ ((aligned (64)));
    double *mtWork, deltaDeg, gammaDeg, kappaDeg, sigmaDeg, thetaDeg;
    int i, ierr, ierr1, ib, ig, indx, im, imt, ik, is, it, nmtBase;
    const double M01[1] = {1.0};
    const double pi180i = 180.0/M_PI;
    const double betaMin = 0.0;
    const double betaMax = M_PI;
    const double gammaMin =-M_PI/6.0;
    const double gammaMax = M_PI/6.0;
    const double kappaMin = 0.0;
    const double kappaMax = 2.0*M_PI;
    const double thetaMin = 0.0;
    const double thetaMax = M_PI_2;
    const double sigmaMin =-M_PI_2;
    const double sigmaMax = M_PI_2;
    //const double sqrt2i = 1.0/sqrt(2.0); // unit magnitude MT
    //------------------------------------------------------------------------//
    ierr = 0;
    // Verify the inputs
    if (ldm < 6 || ng < 1 || nb < 1 || nm < 1 || nk < 1 || nt < 1 || ns < 1 ||
        nmt != ng*nb*nm*nk*nt*ns ||
        gammas == NULL || betas == NULL || M0s == NULL || kappas == NULL ||
        thetas == NULL || sigmas == NULL)
    {
        if (ldm < 6){printf("%s: Invalid leading dimension\n", fcnm);}
        if (ng < 1){printf("%s: ng must be positive\n", fcnm);}
        if (nb < 1){printf("%s: nb must be positive\n", fcnm);}
        if (nm < 1){printf("%s: nm must be positive\n", fcnm);}
        if (nk < 1){printf("%s: nk must be positive\n", fcnm);}
        if (nt < 1){printf("%s: nt must be positive\n", fcnm);}
        if (ns < 1){printf("%s: ns must be positive\n", fcnm);}
        if (nmt != ng*nb*nm*nk*nt*ns)
        {
            printf("%s: nmt != ng*nb*nm*nk*nt*ns %d %d\n",
                   fcnm, nmt, ng*nb*nm*nk*nt*ns);
        }
        if (gammas == NULL){printf("%s: gammas is NULL\n", fcnm);}
        if (betas == NULL){printf("%s: betas is NULL\n", fcnm);}
        if (M0s == NULL){printf("%s: M0s is NULL\n", fcnm);}
        if (kappas == NULL){printf("%s: kappas is NULL\n", fcnm);}
        if (thetas == NULL){printf("%s: thetas is NULL\n", fcnm);}
        if (sigmas == NULL){printf("%s: sigmas is NULL\n", fcnm);}
        return -1;
    }
    for (i=0; i<ng; i++)
    {
        if (gammas[i] < gammaMin || gammas[i] > gammaMax)
        {
            printf("%s: gammas %f out or range [%f,%f]\n",
                   fcnm, gammas[i], gammaMin, gammaMax);
            return -1;
        }
    }
    for (i=0; i<nb; i++)
    {
        if (betas[i] < betaMin || betas[i] > betaMax)
        {
            printf("%s: betas %f out or range [%f,%f]\n",
                   fcnm, betas[i], betaMin, betaMax);
            return -1;
        }
    }
    for (i=0; i<nk; i++)
    {
        if (kappas[i] < kappaMin || kappas[i] > kappaMax)
        {
            printf("%s: kappas %f out or range [%f,%f]\n",
                   fcnm, kappas[i], kappaMin, kappaMax);
            return -1;
        }
    }
    for (i=0; i<nt; i++)
    {
        if (thetas[i] < thetaMin || thetas[i] > thetaMax)
        {
            printf("%s: thetas %f out or range [%f,%f]\n",
                   fcnm, thetas[i], thetaMin, thetaMax);
            return -1;
        }
    }
    for (i=0; i<ns; i++)
    {
        if (sigmas[i] < sigmaMin || sigmas[i] > sigmaMax)
        {
            printf("%s: sigmas %f out or range [%f,%f]\n",
                   fcnm, sigmas[i], sigmaMin, sigmaMax);
            return -1;
        }
    }
    // The moment tensors are equivalent up to a scaling factor - so compute
    // the unit moment tensor
    nmtBase = nb*ng*nk*ns*nt;
    if (nm == 1)
    {
        mtWork = mts;
    }
    else
    {
        mtWork = memory_calloc64f(nmtBase*ldm);
    }
#ifdef _OPENMP
    #pragma omp parallel for collapse(5) \
     firstprivate (lam, M01, Muse, U) \
     private (deltaDeg, gammaDeg, kappaDeg, \
              sigmaDeg, thetaDeg, ierr1, imt, ib, ig, ik, is, it) \
     shared (fcnm, M0s, betas, gammas, kappas, sigmas, mtWork, thetas) \
     reduction (max:ierr) \
     default (none)
#endif
    // Loop on colatitudes
    for (ib=0; ib<nb; ib++)
    {
        // Loop on longitudes
        for (ig=0; ig<ng; ig++)
        {
            // Loop on strike
            for (ik=0; ik<nk; ik++)
            {
                // Loop on slip
                for (is=0; is<ns; is++)
                {
                    // Loop on dip
                    for (it=0; it<nt; it++)
                    {
                        // tape**2 space term
                        //M0 = sqrt2i;
                        deltaDeg = (M_PI/2.0 - betas[ib])*pi180i;
                        gammaDeg = gammas[ig]*pi180i;
                        kappaDeg = kappas[ik]*pi180i;
                        sigmaDeg = sigmas[is]*pi180i;
                        thetaDeg = thetas[it]*pi180i;
                        // index
                        imt = ib*ng*nk*ns*nt
                            + ig*nk*ns*nt
                            + ik*ns*nt
                            + is*nt
                            + it;
                        // Compute the corresponding moment tensor.
                        // Silver and Jordan compute M0 = 1/sqrt(2)norm(M)
                        // hence, if it is the scalar moment computation
                        // that handles the 1/sqrt(2) which means I should
                        // use unity here.  This has been verified b/c
                        // apply CMT2m0 on this moment tensor will yield
                        // the input M0.
                        ierr1 = compearth_TT2CMT(1, &gammaDeg, &deltaDeg, M01,
                                                 &kappaDeg, &thetaDeg,
                                                 &sigmaDeg,
                                                 Muse, lam, U);
/*
                        ierr1 = compearth_tt2cmt(gammaDeg, deltaDeg, 1.0, //sqrt2i,
                                                 kappaDeg, thetaDeg,
                                                 sigmaDeg,
                                                 Muse, lam, U);
*/
                        // Convert from USE to our NED estimation basis
                        ierr1 = compearth_convertMT(1, CE_USE, CE_NED, Muse,
                                                    &mtWork[imt*ldm]);
                        if (ierr1 != 0)
                        {
                            printf("%s: Error changing coords\n", fcnm);
                            ierr = ierr + 1;
                            continue;
                        }
                    } // loop on dip
                } // loop on slip
            } // loop on strike
        } // loop on longitudes
    } // loop on latitudes
    if (ierr != 0)
    {
        printf("%s: Errors during mt computation\n", fcnm);
        ierr = 1;
    }
    // Now copy all moment tensors over
    if (nm == 1)
    {
        cblas_dscal(ldm*nmtBase, M0s[0], mtWork, 1);
    }
    else
    {
        for (im=0; im<nm; im++)
        {
            indx = im*nmtBase*ldm; //nb*ng*nk*ns*nt;
            for (imt=0; imt<ldm*nmtBase; imt++)
            {
                mts[indx+imt] = M0s[im]*mtWork[imt];
            }
        }
        memory_free64f(&mtWork);
    }
    mtWork = NULL;
    return ierr;
}


text:
#ifndef __GSL_BLOCK_H__
#define __GSL_BLOCK_H__

#include <gsl/block/gsl_block_double.h>

#endif /* __GSL_BLOCK_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#include <stdio.h>
#include <math.h>
#include "wf.h"
#include "alphas.h"
#include "glasma.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

//workspace and global variables for integration
static double abserr = 1.e-200;
static double relerr = 1.e-20;

double phiKernel(double x, void *params);

static double result,error;
static size_t neval;
static gsl_function phiIntKernel;

double d2Nglasma0(double pT, double qT, double phipq, double yp, double yq, double rts)
{
    static struct glasma_params params;
    params.pT = pT;
    params.qT = qT;
    params.yp = yp;
    params.yq = yq;
    params.rts = rts;
    params.phipq = phipq;
    phiIntKernel.function = &phiKernel;
    (params.Kernel).function = &doubleKernel;

    phiIntKernel.params = &params;
    gsl_integration_qng(&phiIntKernel, -M_PI, M_PI, abserr, relerr,\
            &result, &error, &neval);

    double norm = (Nc*Nc)/gsl_pow_3(Nc*Nc-1.)/(4.*pow(M_PI,10.))/gsl_pow_2(pT*
        qT)*alpha(pT)*alpha(qT);

    return norm*result;
}

double d2Nglasma1(double pT, double yp, double yq, double rts)
{
    double sinres, cosres;
    static struct glasma_params params;
    params.pT = pT;
    params.yp = yp;
    params.yq = yq;
    params.rts = rts;

    phiIntKernel.function = &phiKernel;
    (params.Kernel).function = &cosKernel;
    phiIntKernel.params = &params;
    gsl_integration_qng(&phiIntKernel, -M_PI, M_PI, abserr, relerr, &cosres, &error, &neval);

    (params.Kernel).function = &sinKernel;
    phiIntKernel.params = &params;
    gsl_integration_qng(&phiIntKernel, -M_PI, M_PI, abserr, relerr, &sinres, &error, &neval);

    double norm =  (Nc*Nc)/gsl_pow_3(Nc*Nc-1.)/(4.*pow(M_PI,10.))/gsl_pow_2(pT*pT)*alpha(pT)*alpha(pT) ;

    return norm*(cosres*cosres + sinres*sinres);
}

double doubleKernel(double x, void *p)
{
   struct glasma_params params = *(struct glasma_params *)p;

   double rts = params.rts;
   double pT, qT, yp, yq;
   double x1p, x2p, x1q, x2q, x1max, x2max;

   double phi = params.phi;
   double phipq = params.phipq;
   double kT =  x;

   pT = params.pT;
   qT = params.qT;
   yp = params.yp;
   yq = params.yq;

   x1p = pT/rts*exp(+yp);
   x2p = pT/rts*exp(-yp);
   x1q = qT/rts*exp(+yq);
   x2q = qT/rts*exp(-yq);

   x1max = gsl_max(x1p,x1q);
   x2max = gsl_max(x2p,x2q);

   double qTmkT = sqrt( qT*qT + kT*kT - 2.0*qT*kT*cos(phipq-phi) );
   double qTpkT = sqrt( qT*qT + kT*kT + 2.0*qT*kT*cos(phipq-phi) );
   double pTmkT = sqrt( pT*pT + kT*kT - 2.0*pT*kT*cos(phi) );
   double pTpkT = sqrt( pT*pT + kT*kT + 2.0*pT*kT*cos(phi) );

   double pTmkTmqT = sqrt( pT*pT + kT*kT + qT*qT \
                          -2.0*pT*qT*cos(phipq) \
                          -2.0*pT*kT*cos(phi) \
                          +2.0*qT*kT*cos(phipq-phi)\
                        );

   double pTmkTpqT = sqrt( pT*pT + kT*kT + qT*qT \
                          +2.0*pT*qT*cos(phipq) \
                          -2.0*pT*kT*cos(phi) \
                          -2.0*qT*kT*cos(phipq-phi)\
                        );

   //diagram A & E
   //this factor of 0.5 is from
   //symmetrizing wrt pTmkT
   double AE = 0.5*kT*pow(wf(1,x1max,kT),2.0)\
          *( wf(2,x2p,pTmkT) + wf(2,x2p,pTpkT) )\
          *( wf(2,x2q,qTmkT) + wf(2,x2q,qTpkT) );

   //diagram B & F
   double BF = 0.5*kT*pow(wf(2,x2max,kT),2.0)\
          *( wf(1,x1p,pTmkT) + wf(1,x1p,pTpkT) )\
          *( wf(1,x1q,qTmkT) + wf(1,x1q,qTpkT) );

   double T1D, T2D, T1H, T2H, denD, denH;
   T1D = ( kT*pT*cos(phi)-kT*kT )\
        *( pT*pT - pT*qT*cos(phipq) - pT*kT*cos(phi) - pow(pTmkTmqT,2.0) ) \
        - (pT*kT*sin(phi))*(pT*qT*sin(phipq) + pT*kT*sin(phi) );

   T1H = ( kT*pT*cos(phi)-kT*kT )\
        *( pT*pT + pT*qT*cos(phipq) - pT*kT*cos(phi) - pow(pTmkTpqT,2.0) ) \
        - (pT*kT*sin(phi))*(-pT*qT*sin(phipq) + pT*kT*sin(phi) );

   T2D = ( kT*qT*cos(phipq-phi)+kT*kT )\
        *( -qT*qT + pT*qT*cos(phipq) - qT*kT*cos(phipq-phi) + pow(pTmkTmqT,2.0) ) \
        + (qT*kT*sin(phipq-phi))*(pT*qT*sin(phipq) - qT*kT*sin(phipq-phi) );

   T2H = ( -kT*qT*cos(phipq-phi)+kT*kT )\
        *( -qT*qT - pT*qT*cos(phipq) + qT*kT*cos(phipq-phi) + pow(pTmkTpqT,2.0) ) \
        + (qT*kT*sin(phipq-phi))*(pT*qT*sin(phipq) - qT*kT*sin(phipq-phi) );

   denD = pow(kT*pTmkTmqT*pTmkT*qTpkT,2.0);

   denH = pow(kT*pTmkTpqT*pTmkT*qTmkT,2.0);

   //this factor of 0.5 is from the different f^4\delta^4 structure
   double D = 0.5*kT*T1D*T2D/denD\
               *wf(1,x1max,kT)*wf(1,x1max,pTmkTmqT)*wf(2,x2max,pTmkT)*wf(2,x2max,qTpkT);

   double H = 0.5*kT*T1H*T2H/denH\
               *wf(1,x1max,kT)*wf(1,x1max,pTmkTpqT)*wf(2,x2max,pTmkT)*wf(2,x2max,qTmkT);

   return (AE + BF + D + H);
}

double cosKernel(double x, void *p)
{
   struct glasma_params params = *(struct glasma_params *)p;

   double rts = params.rts;
   double pT, yp, yq;
   double x1p, x2p, x1q, x2q, x1max, x2max;

   double phi = params.phi;
   double kT =  x;

   pT = params.pT;
   yp = params.yp;
   yq = params.yq;

   x1p = pT/rts*exp(+yp);
   x2p = pT/rts*exp(-yp);
   x1q = pT/rts*exp(+yq);
   x2q = pT/rts*exp(-yq);

   x1max = gsl_max(x1p,x1q);
   x2max = gsl_max(x2p,x2q);

   double pTmkT = sqrt( pT*pT + kT*kT - 2.0*pT*kT*cos(phi) );

   double num = pow( kT*pT*cos(phi)-kT*kT, 2.0);
   double den  = pow(kT*pTmkT,2.);

   return kT*wf(1,x1p,kT)*wf(2,x2p,pTmkT)*num/den;
}

double sinKernel(double x, void *p)
{
   struct glasma_params params = *(struct glasma_params *)p;

   double rts = params.rts;
   double pT, yp, yq;
   double x1p, x2p, x1q, x2q, x1max, x2max;

   double phi = params.phi;
   double kT =  x;

   pT = params.pT;
   yp = params.yp;
   yq = params.yq;

   x1p = pT/rts*exp(+yp);
   x2p = pT/rts*exp(-yp);
   x1q = pT/rts*exp(+yq);
   x2q = pT/rts*exp(-yq);

   x1max = gsl_max(x1p,x1q);
   x2max = gsl_max(x2p,x2q);

   double pTmkT = sqrt( pT*pT + kT*kT - 2.0*pT*kT*cos(phi) );

   double num = pow( kT*pT*sin(phi), 2.0);
   double den = pow( kT*pTmkT ,2.0);

   return kT*wf(1,x1max,kT)*wf(2,x2max,pTmkT)*num/den;
}

//does phi integration over d^2k_\perp
double phiKernel(double x, void *p)
{
   struct glasma_params params = *(struct glasma_params *)p;
   double result, error;
   params.phi = x;

   gsl_function F = params.Kernel;
   F.params = &params;

   gsl_integration_qng(&F, kT_min, 2.0*kT_max, abserr, relerr,\
           &result, &error, &neval);

return result;
}

void TabulateGlasma(FILE *out, double rts)
{
   //returns dN/d^2p_T d^2q_T dy_p dy_q / (S_\perp) [GeV^-2]
   //this is the delta function contribution
   //of the glasma graphs

   double yp, yq, rtpT;
   //double yrange = 6.0;
   double dy = 0.25;
   double yrange = dy;

   for (yp = -yrange; yp <= yrange; yp += dy)
   for (yq = -yrange; yq <= yrange; yq += dy)
   for (rtpT = .1; rtpT <= 5.11; rtpT += .2){
   fprintf(out,"%10.2f\t%10.2f\t%10.2f\t%10.5e\n", yp, yq, rtpT,\
   d2Nglasma1(pow(rtpT,2.),yp,yq,rts) );
   }
   fclose(out);

return;
}

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc


text:
/******************************************************************************
 * gsl_sprng.c: rewrite of gsl-sprng.h to use sprng streams                   *
 * Author: Yan Y. Liu <yanliu@illinois.edu>                                   *
 * Date: 2014/08/17                                                           *
 * Copyright and license of this source file are specified in LICENSE.TXT     *
 * under the root directory of this software package.                         *
 ******************************************************************************/
#ifndef GSL_SPRNG_C
#define GSL_SPRNG_C

#include <sys/time.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "gsl-sprng.h"
/***************************/
/* gsl definition of sprng */
/***************************/
void sprng_set(void * vstate, unsigned long int seed)
{
    gsl_sprng_state_t *mysprng = (gsl_sprng_state_t*) vstate;

    if (mysprng->stream != NULL) {
        free_sprng(mysprng->stream);
    }
    mysprng->stream = init_sprng(
                DEFAULT_RNG_TYPE,
                mysprng->streamnum,
                mysprng->nstreams,
                seed,
                SPRNG_DEFAULT);
}

unsigned long sprng_get(void * vstate)
{
    return( (long) isprng( ((gsl_sprng_state_t*) vstate)->stream ) );
}

double sprng_get_double(void * vstate)
{
    return( (double) sprng( ((gsl_sprng_state_t*) vstate)->stream ) );
}

/* global sprng 2.0 type for gsl initialization */
const gsl_rng_type gsl_rng_sprng20 = {
    "sprng20",        /* name */
    0x7fffffffUL,     /* RAND_MAX */
    0,                /* RAND_MIN */
    0,                /* size of state - not sure about this */
    &sprng_set,          /* initialisation */
    &sprng_get,          /* get integer RN */
    &sprng_get_double    /* get double RN */
};

/***************************/
/*      gsl_sprng API      */
/***************************/

// if stream=NULL, create a new stream
// otherwise, just use the stream
gsl_rng *gsl_sprng_alloc (int streamnum, int nstreams, int seed, int * stream)
{
    gsl_sprng_state_t *  gsl_sprng_state =
   	    (gsl_sprng_state_t*)  malloc( sizeof( gsl_sprng_state_t ) );

    if (gsl_sprng_state == NULL) {
        GSL_ERROR_VAL (
	   	  "failed to allocate space for gsl_sprng_state_t struct",
                     GSL_ENOMEM, 0);
        return NULL; // not reachable
    }
    gsl_sprng_state->streamnum = streamnum;
    gsl_sprng_state->nstreams = nstreams;
    // set seed
    if (seed == -1) {
        struct timeval tv; struct timezone tz;
        gettimeofday(&tv, &tz);
        seed = (int)(tv.tv_sec * streamnum);
    } else if (seed == 0) {
        seed = gsl_rng_default_seed; // zero, too, in fact
    }
    gsl_sprng_state->seed = seed;
    // sprng init
    if (stream == NULL) {
        stream = init_sprng(
                DEFAULT_RNG_TYPE,
                streamnum,
                nstreams,
                seed,
                SPRNG_DEFAULT);
    }
    gsl_sprng_state->stream = stream;
    // gsl_rng init
    gsl_rng * r = (gsl_rng *) malloc (sizeof (gsl_rng));
    if (r == NULL) {
        GSL_ERROR_VAL ("gsl_sprng: failed to allocate space for gsl_rng struct",
                        GSL_ENOMEM, 0);
        free(gsl_sprng_state); // not reachable
        return NULL;
    }
    r->type = &gsl_rng_sprng20;
    r->state = gsl_sprng_state;

    return r;
}

void gsl_sprng_free (gsl_rng * r)
{
    if (r==NULL) return;
    gsl_sprng_state_t *  gsl_sprng_state = (gsl_sprng_state_t *) r->state;
    if (gsl_sprng_state != NULL) {
        int numAvail = free_sprng(gsl_sprng_state->stream);
        fprintf(stdout, "gsl_sprng: free random stream. %d streams left\n", numAvail);
        free(r->state);
    }
    free(r);
}

void gsl_sprng_print (gsl_rng *r)
{
    if (r==NULL) return;
    gsl_sprng_state_t *  gsl_sprng_state = (gsl_sprng_state_t *) r->state;
    if (gsl_sprng_state == NULL) return;
    printf("gsl_sprng: streamnum=%d, nstreams=%d, seed=%d\n", gsl_sprng_state->streamnum, gsl_sprng_state->nstreams, gsl_sprng_state->seed);
    int * stream = gsl_sprng_state->stream;
    if (stream == NULL) return;
    print_sprng(stream);
}

int gsl_sprng_fread (char * fpath, gsl_rng * r)
{
    // delete old stream
    if (r!=NULL) {
        gsl_sprng_free(r);
        r = NULL;
    }
    // allocate memory to load data
    char buf[MAX_PACKED_LENGTH];
    // read stream from file to array
    FILE * f;
    int size;
    if ((f=fopen(fpath, "r")) == NULL) {
        fprintf(stderr, "gsl_sprng_fread(): failed to open %s for reading\n", fpath);
        return 0;
    }
    fread(&size,sizeof(int),1,f);
    fprintf(stdout, "sprng_fread[%s]: size=%d\n", fpath, size);
    int numRead = fread(buf, sizeof(char), size, f);
    if (numRead != size) {
        fprintf(stderr, "gsl_sprng_fread(): failed to read stream of size %d\n", size);
        fclose(f);
        return 0;
    }
    fclose(f);
    // load stream from array
    int * stream = unpack_sprng(buf);
    if (stream == NULL) {
        fprintf(stderr, "gsl_sprng_fread(): failed to read stream of size %d\n", size);
        return 0;
    }
    print_sprng(stream);
    // create new gsl_rng
    r = gsl_sprng_alloc(0, 0, 0, stream);
    if (r==NULL) return 0;
    return 1;
}

int gsl_sprng_fwrite (char * fpath, const gsl_rng * r)
{
    if (r==NULL) return 0;
    gsl_sprng_state_t *  gsl_sprng_state = (gsl_sprng_state_t*) r->state;
    if (gsl_sprng_state == NULL) return 0;
    int * stream = gsl_sprng_state->stream;
    if (stream==NULL) return 0;
    // write state to array
    char *buf;
    int size = pack_sprng(stream, &buf);
    // write array to file
    FILE * f;
    if ((f=fopen(fpath, "w")) == NULL) {
        fprintf(stderr, "gsl_sprng_fwrite(): failed to open %s for writing\n", fpath);
        return 0;
    }
    // our format: filesize, filecontent
    fwrite(&size,sizeof(int),1,f);
    int numWrote = fwrite(buf, sizeof(char), size, f);
    if (numWrote != size) {
        fprintf(stderr, "gsl_sprng_fwrite(): failed to write stream of size %d\n", size);
    }
    fclose(f);
    // sprng uses malloc() to allocate buf, we must free it properly
    free(buf);

    return 1;
}

#endif


text:
/**
 * @file batchf_zgemm_stride.c
 *
 *  @brief BBLAS batchf_zgemm_stride double _Complex routine.
 *
 *  BBLAS is a software package provided by Univ. of Manchester,
 *  Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-04-06
 *
 **/

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @precisions normal z -> c d s
 **/
#endif

#include <cblas.h>
#include "bblas.h"

#define COMPLEX

/**
    Purpose
    -------
    <b>batchf_zgemm_stride</b> is a batch version of zgemm for fixed size.
	It performs the matrix-matrix operations

	arrayC[i] = alpha[i]*op( arrayA[i] )*op( arrayB[i] ) + beta[i]*arrayC[i],

	where op( X ) is one of

	op( X ) = X      or
	op( X ) = X**T   or
	op( X ) = X**H,

	alpha[i] and beta[i] are scalars, and arrayA[i], arrayB[i] and C are matrices, with
	op( arrayA[i] ) an m by k matrix, op( arrayB[i] ) a k by n matrix and arrayC[i]
	an m by n matrix.

	Fixed Batch Operations
	-----------------------------------

	- all parameters that are arrays must have length at least batch_count.
	- all parameters that are arrays must have all values set.

	the values of transA, transB, M, N, K, alpha, beta, lda, ldb, and ldc
    are used for all computations.

	Parameters
	----------
	@param[in]
	transA   <tt>enum BBLAS_TRANS</tt>.
	         transA specifies the form of op( arrayA[i] ) to be used in
			 the matrix multiplication as follows:
      -     = BblasNoTrans:    op( arrayA[i] ) = arrayA[i].
      -     = BblasTrans:      op( arrayA[i] ) = arrayA[i]**T.
      -     = BblasConjTrans:  op( arrayA[i] ) = arrayA[i]**H.

	@param[in]
	transB   <tt>enum BBLAS_TRANS</tt>.
	         On entry, transB specifies the form of op( arrayB[i] ) to be used in
             the matrix multiplication as follows:
      -     = BblasNoTrans:    op( arrayB[i] ) = arrayB[i].
      -     = BblasTrans:      op( arrayB[i] ) = arrayB[i]**T.
      -     = BblasConjTrans:  op( arrayB[i] ) = arrayB[i]**H.

	@param[in]
	M       <tt>int</tt>.
            M specifies the number of rows of the matrix
            op( arrayA[i] ) and of the matrix arrayC[i].
            M must be greater than zero.

	@param[in]
	N       <tt>int</tt>.
            N specifies the number of columns of the matrix
            op( arrayB[i] ) and the number of columns of the matrix arrayC[i].
 		    N must be greater than zero.

    @param[in]
    K      <tt>int</tt>.
	       K specifies the number of columns of the matrix
           op( arrayA[i] ) and the number of rows of the matrix op( arrayB[i] ).
		   K must be greater than zero.

    @param[in]
    alpha   <tt>complex_16</tt>.

    @param[in]
    arrayA  Array of pointers.
            Each element arrayA[i] is a pointer to a COMPLEX_16 matrix of
 		    dimension lda by Ka, where Ka is K when transA = BblasNoTrans,
 		    and is M otherwise.
            When using transA = BblasNoTrans the leading M by K
            part of arrayA[i] must contain the matrix elements, otherwise
            the leading K by M part of arrayA[i] must contain the
            matrix elements.

    @param[in]
    lda     <tt>int</tt>.
            lda specifies the first dimension of arrayA[i] as declared
            in the calling (sub) program. When transA = BblasNoTrans then
            lda must be at least max( 1, M ), otherwise lda must be at
            least max( 1, K ).

    @param[in]
    strideA <tt>int</tt>.
            strideA specifies the storage spacing between A matrices.
                    strideA must be at least M by K.

    @param[in]
    arrayB  Array of pointers.
  		    Each element arrayB[i] is a pointer to a COMPLEX_16 matrix of
  		    dimension ldb by Kb, where Kb is N when transB = BblasNoTrans,
  		    and is K otherwise.
            When using transB = BblasNoTrans the leading K by N
            part of arrayB[i] must contain the matrix elements, otherwise
            the leading N by K part of arrayB[i] must contain the
            matrix elements.

    @param[in]
    ldb     <tt>int</tt>.
  		    ldb specifies the first dimension of arrayB[i] as declared
            in the calling (sub) program. When transB = BblasNoTrans then
            ldb must be at least max( 1, K ), otherwise ldb must be at
            least max( 1, N ).

    @param[in]
    strideB <tt>int</tt>.
            strideB specifies the storage spacing between A matrices.
            strideB must be at least K by N.

    @param[in]
    beta    <tt>complex_16</tt>.
            When beta is set to zero arrayC[i] need not be set on input.

    @param[in,out]
    arrayC  Array of pointers.
	        Each element arrayC[i] is a pointer to a COMPLEX_16 matrix of
      	    dimension ldc by N.
  	        Before entry, the leading M by N part of the arrayC[i] must
            contain a matrix C, except when beta is zero, in which
            case C need not be set on entry.
            On exit, the matrix arrayC[i] is overwritten by the M by N matrix
            ( alpha*op( arrayA[i] )*op( arrayB[i] ) + beta*arrayC[i] ).

    @param[in]
    ldc     <tt>int</tt>.
  	        Each element ldc specifies the first dimension of arrayC[i] as declared
  	        in the calling (sub) program. The value ldc must be at least
  	        max( 1, M )

    @param[in]
    strideC <tt>int</tt>.
            strideC specifies the storage spacing between A matrices.
            strideC must be at least M by N.

    @param[in]
    batch_count  <tt>int</tt>
                 The number of matrices to operate on.

    @param[in,out]
    info    <tt>int</tt>.
            info is the error return code of the batch.
  	        The error codes can be found in bblas_macros.h.
 **/
void batchf_zgemm_stride(
    const enum BBLAS_TRANS transA, const enum BBLAS_TRANS transB,
    const int M,  const int N, const int K,
    const BBLAS_Complex64_t alpha,
    const BBLAS_Complex64_t *arrayA, const int lda, const int strideA,
    const BBLAS_Complex64_t *arrayB, const int ldb, const int strideB,
    const BBLAS_Complex64_t beta,
    BBLAS_Complex64_t *arrayC, const int ldc, const int strideC,
    const int batch_count, int *info)
{
    /* Local variables */
    int batch_iter;
    int LDA,  LDB;
    char func_name[15] = "batchf_zgemm_op";

	/* Check input arguments */
    if (batch_count < 0)
    {
	    xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
    }
    if ((transA != BblasNoTrans) &&
        (transA != BblasTrans) &&
        (transA != BblasConjTrans))
    {
        xerbla_batch(func_name, BBLAS_ERR_TRANSA, -1);
        *info  = BBLAS_ERR_TRANSA;
        return;
    }
    if ((transB != BblasNoTrans) &&
        (transB != BblasTrans) &&
        (transB != BblasConjTrans))
    {
        xerbla_batch(func_name, BBLAS_ERR_TRANSB, -1);
        *info = BBLAS_ERR_TRANSB;
        return;
    }
    if ( transA == BblasNoTrans )
    {
        LDA = M;
    } else
    {
        LDA = K;
    }
    if ( transB == BblasNoTrans )
    {
        LDB = K;
    } else
    {
        LDB = N;
    }
    if (M < 0)
    {
        xerbla_batch(func_name, BBLAS_ERR_M, -1);
        *info = BBLAS_ERR_M;
        return;
    }
    if (N < 0) {
        xerbla_batch(func_name, BBLAS_ERR_N, -1);
        *info = BBLAS_ERR_N;
        return;
    }
    if (K < 0) {
        xerbla_batch(func_name, BBLAS_ERR_K, -1);
        *info = BBLAS_ERR_K;
        return;
    }
    if (lda < max(1, LDA)) {
        xerbla_batch(func_name, BBLAS_ERR_LDA, -1);
        *info =  BBLAS_ERR_LDA;
        return;
    }
    if (ldb < max(1, LDB)) {
        xerbla_batch(func_name, BBLAS_ERR_LDB, -1);
        *info = BBLAS_ERR_LDB;
        return;
    }
    if (ldc < max(1, M)) {
        xerbla_batch(func_name, BBLAS_ERR_LDC, -1);
        *info = BBLAS_ERR_LDC;
        return;
    }
    if (strideA < M*K) {
        xerbla_batch(func_name, BBLAS_ERR_STRIDEA, -1);
        *info =  BBLAS_ERR_STRIDEA;
        return;
    }
    if (strideB < K*N) {
        xerbla_batch(func_name, BBLAS_ERR_STRIDEB, -1);
        *info =  BBLAS_ERR_STRIDEB;
        return;
    }
    if (strideC < M*N) {
        xerbla_batch(func_name, BBLAS_ERR_STRIDEC, -1);
        *info =  BBLAS_ERR_STRIDEC;
        return;
    }

    /* Skip subproblems where nothing needs to be done */
    if (M == 0 || N == 0     ||
	    ((alpha == (BBLAS_Complex64_t)0.0 || K == 0) &&
	    beta == (BBLAS_Complex64_t)1.0 ))
    {
        *info =  BBLAS_SUCCESS;
        return;
    }
    for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
    {
        /* Call to cblas_zgemm */
        cblas_zgemm(
        BblasColMajor,
        transA,
        transB,
        M,
        N,
        K,
        CBLAS_SADDR(alpha),
        //&arrayA[batch_iter*M*K],
        &arrayA[batch_iter*strideA],
        lda,
        //&arrayB[batch_iter*K*N],
        &arrayB[batch_iter*strideB],
        ldb,
        CBLAS_SADDR(beta),
        //&arrayC[batch_iter*M*N],
        &arrayC[batch_iter*strideC],
        ldc);
    }
    /* Successful */
    *info = BBLAS_SUCCESS;
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group