text:
/*! \file allvars.h
 *  \brief declares global variables.
 *
 *  This file declares all global variables. Further variables should be added here, and declared as
 *  'extern'. The actual existence of these variables is provided by the file 'allvars.c'. To produce
 *  'allvars.c' from 'allvars.h', do the following:
 *
 *     - Erase all #define statements
 *     - add #include "allvars.h"
 *     - delete all keywords 'extern'
 *     - delete all struct definitions enclosed in {...}, e.g.
 *        "extern struct global_data_all_processes {....} All;"
 *        becomes "struct global_data_all_processes All;"
 */

#ifndef ALLVARS_H
#define ALLVARS_H

#include <mpi.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#ifdef MPISENDRECV_CHECKSUM
#define MPI_Sendrecv MPI_Check_Sendrecv
#endif

#ifdef MPISENDRECV_SIZELIMIT
#define MPI_Sendrecv MPI_Sizelimited_Sendrecv
#endif

#include "tags.h"
#ifdef CHEMISTRY
#include "chemistry.h"
#endif
#include "assert.h"

#define  GADGETVERSION   "3.0"	/*!< code version string */

#define  GENERATIONS     2	/*!< Number of star particles that may be created per gas particle */

#define  TIMEBINS        29

#define  TIMEBASE        (1<<TIMEBINS)	/*!< The simulated timespan is mapped onto the integer interval [0,TIMESPAN],
					 *   where TIMESPAN needs to be a power of 2. Note that (1<<28) corresponds
					 *   to 2^29
					 */
#ifndef  MULTIPLEDOMAINS
#define  MULTIPLEDOMAINS     1
#endif

#ifndef  TOPNODEFACTOR
#define  TOPNODEFACTOR       2.5
#endif

#define  NODELISTLENGTH      8

typedef unsigned long long peanokey;

#define  BITS_PER_DIMENSION 21	/* for Peano-Hilbert order. Note: Maximum is 10 to fit in 32-bit integer ! */
#define  PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION))

#ifndef  GAMMA
#define  GAMMA         (5.0/3)	/*!< adiabatic index of simulated gas */
#endif

#define  GAMMA_MINUS1  (GAMMA-1)

#define  HYDROGEN_MASSFRAC 0.76	/*!< mass fraction of hydrogen, relevant only for radiative cooling */

#define  METAL_YIELD       0.02	/*!< effective metal yield for star formation */

#define  MAX_REAL_NUMBER  1e37
#define  MIN_REAL_NUMBER  1e-37

#define  RNDTABLE 8192

/* ... often used physical constants (cgs units) */

#define  GRAVITY     6.672e-8
#define  SOLAR_MASS  1.989e33
#define  SOLAR_LUM   3.826e33
#define  RAD_CONST   7.565e-15
#define  AVOGADRO    6.0222e23
#define  BOLTZMANN   1.3806e-16
#define  GAS_CONST   8.31425e7
#define  C           2.9979e10
#define  PLANCK      6.6262e-27
#define  CM_PER_MPC  3.085678e24
#define  PROTONMASS  1.6726e-24
#define  ELECTRONMASS 9.10953e-28
#define  THOMPSON     6.65245e-25
#define  ELECTRONCHARGE  4.8032e-10
#define  HUBBLE          3.2407789e-18	/* in h/sec */
#define  LYMAN_ALPHA      1215.6e-8	/* 1215.6 Angstroem */
#define  LYMAN_ALPHA_HeII  303.8e-8	/* 303.8 Angstroem */
#define  OSCILLATOR_STRENGTH       0.41615
#define  OSCILLATOR_STRENGTH_HeII  0.41615

#ifdef NAVIERSTOKES
#define  LOG_LAMBDA      37.8	/* logarithmic Coulomb factor */
#endif

#ifdef CHEMISTRY
#define  T_CMB0      2.728	/* present-day CMB temperature */
#endif

#define  SEC_PER_MEGAYEAR   3.155e13
#define  SEC_PER_YEAR       3.155e7

/*Determines the maximum size of arrays related to the number of CR populations */
#ifndef NUMCRPOP   /*!< Number of CR populations pressent in parameter file */
#define NUMCRPOP 1
#endif

#ifndef FOF_PRIMARY_LINK_TYPES
#define FOF_PRIMARY_LINK_TYPES 2
#endif

#ifndef FOF_SECONDARY_LINK_TYPES
#define FOF_SECONDARY_LINK_TYPES 0
#endif

/* some flags for the field "flag_ic_info" in the file header */
#define FLAG_ZELDOVICH_ICS     1
#define FLAG_SECOND_ORDER_ICS  2
#define FLAG_EVOLVED_ZELDOVICH 3
#define FLAG_EVOLVED_2LPT      4
#define FLAG_NORMALICS_2LPT    5

#ifndef ASMTH
/*! ASMTH gives the scale of the short-range/long-range force split in units of FFT-mesh cells */
#define ASMTH 1.25
#endif
#ifndef RCUT
/*! RCUT gives the maximum distance (in units of the scale used for the force split) out to which short-range
 * forces are evaluated in the short-range tree walk.
 */
#define RCUT  4.5
#endif

#define COND_TIMESTEP_PARAMETER 0.25
#define VISC_TIMESTEP_PARAMETER 0.25

#define MAXLEN_OUTPUTLIST 350	/*!< maxmimum number of entries in output list */

#define DRIFT_TABLE_LENGTH  1000	/*!< length of the lookup table used to hold the drift and kick factors */

#define MAXITER 150

#ifndef LINKLENGTH
#define LINKLENGTH 0.2
#endif

#ifndef FOF_GROUP_MIN_LEN
#define FOF_GROUP_MIN_LEN 32
#endif

#define MINRESTFAC 0.05

#ifndef LONGIDS
typedef unsigned int MyIDType;
#else
typedef unsigned long long MyIDType;
#endif

#ifndef DOUBLEPRECISION     /* default is single-precision */
typedef float  MyFloat;
typedef float  MyDouble;
#else
#if (DOUBLEPRECISION == 2)   /* mixed precision */
typedef float   MyFloat;
typedef double  MyDouble;
#else                        /* everything double-precision */
typedef double  MyFloat;
typedef double  MyDouble;
#endif
#endif

#ifdef OUTPUT_IN_DOUBLEPRECISION
typedef double MyOutputFloat;
#else
typedef float MyOutputFloat;
#endif

#ifdef INPUT_IN_DOUBLEPRECISION
typedef double MyInputFloat;
#else
typedef float MyInputFloat;
#endif

struct unbind_data
{
  int index;
};

#ifdef FIX_PATHSCALE_MPI_STATUS_IGNORE_BUG
extern MPI_Status mpistat;
#undef MPI_STATUS_IGNORE
#define MPI_STATUS_IGNORE &mpistat
#endif

#ifdef FLTROUNDOFFREDUCTION
#define FLT(x) ((MyFloat)(x))
#ifdef SOFTDOUBLEDOUBLE      /* this requires a C++ compilation */
#include "dd.h"
typedef dd MyLongDouble;
#else
typedef long double MyLongDouble;
#endif
#else  /* not enabled */
#define FLT(x) (x)
typedef MyFloat MyLongDouble;
#endif  /* end FLTROUNDOFFREDUCTION */

#define CPU_ALL            0
#define CPU_TREEWALK1      1
#define CPU_TREEWALK2      2
#define CPU_TREEWAIT1      3
#define CPU_TREEWAIT2      4
#define CPU_TREESEND       5
#define CPU_TREERECV       6
#define CPU_TREEMISC       7
#define CPU_TREEBUILD      8
#define CPU_TREEUPDATE     9
#define CPU_TREEHMAXUPDATE 10
#define CPU_DOMAIN         11
#define CPU_DENSCOMPUTE    12
#define CPU_DENSWAIT       13
#define CPU_DENSCOMM       14
#define CPU_DENSMISC       15
#define CPU_HYDCOMPUTE     16
#define CPU_HYDWAIT        17
#define CPU_HYDCOMM        18
#define CPU_HYDMISC        19
#define CPU_DRIFT          20
#define CPU_TIMELINE       21
#define CPU_POTENTIAL      22
#define CPU_MESH           23
#define CPU_PEANO          24
#define CPU_COOLINGSFR     25
#define CPU_SNAPSHOT       26
#define CPU_FOF            27
#define CPU_BLACKHOLES     28
#define CPU_MISC           29
#define CPU_SMTHCOMPUTE    30
#define CPU_SMTHWAIT       31
#define CPU_SMTHCOMM       32
#define CPU_SMTHMISC       33
#define CPU_HOTNGBS        34
#define CPU_WEIGHTS_HOT    35
#define CPU_ENRICH_HOT     36
#define CPU_WEIGHTS_COLD   37
#define CPU_ENRICH_COLD    38
#define CPU_CSMISC         39

#define CPU_PARTS          40  /* this gives the number of parts above (must be last) */

#define CPU_STRING_LEN 120

#ifndef  TWODIMS
#define  NUMDIMS 3		/*!< For 3D-normalized kernel */
#define  KERNEL_COEFF_1  2.546479089470	/*!< Coefficients for SPH spline kernel and its derivative */
#define  KERNEL_COEFF_2  15.278874536822
#define  KERNEL_COEFF_3  45.836623610466
#define  KERNEL_COEFF_4  30.557749073644
#define  KERNEL_COEFF_5  5.092958178941
#define  KERNEL_COEFF_6  (-15.278874536822)
#define  NORM_COEFF      4.188790204786	/*!< Coefficient for kernel normalization. Note:  4.0/3 * PI = 4.188790204786 */
#else
#define  NUMDIMS 2		/*!< For 2D-normalized kernel */
#define  KERNEL_COEFF_1  (5.0/7*2.546479089470)	/*!< Coefficients for SPH spline kernel and its derivative */
#define  KERNEL_COEFF_2  (5.0/7*15.278874536822)
#define  KERNEL_COEFF_3  (5.0/7*45.836623610466)
#define  KERNEL_COEFF_4  (5.0/7*30.557749073644)
#define  KERNEL_COEFF_5  (5.0/7*5.092958178941)
#define  KERNEL_COEFF_6  (5.0/7*(-15.278874536822))
#define  NORM_COEFF      M_PI	/*!< Coefficient for kernel normalization. */
#endif

#if defined (BLACK_HOLES) || defined(CS_MODEL) || defined(RADTRANSFER) || defined(SNIA_HEATING)
#define PPP P
#else
#define PPP SphP
#endif

#define DMAX(a,b) (dmax1=(a),dmax2=(b),(dmax1>dmax2)?dmax1:dmax2)
#define DMIN(a,b) (dmin1=(a),dmin2=(b),(dmin1<dmin2)?dmin1:dmin2)
#define IMAX(a,b) (imax1=(a),imax2=(b),(imax1>imax2)?imax1:imax2)
#define IMIN(a,b) (imin1=(a),imin2=(b),(imin1<imin2)?imin1:imin2)

#ifdef PERIODIC
extern MyDouble boxSize, boxHalf;
#ifdef LONG_X
extern MyDouble boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
extern MyDouble boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
extern MyDouble boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif

#ifdef PERIODIC
#define NGB_PERIODIC_LONG_X(x) (xtmp=fabs(x),(xtmp>boxHalf_X)?(boxSize_X-xtmp):xtmp)
#define NGB_PERIODIC_LONG_Y(x) (xtmp=fabs(x),(xtmp>boxHalf_Y)?(boxSize_Y-xtmp):xtmp)
#define NGB_PERIODIC_LONG_Z(x) (xtmp=fabs(x),(xtmp>boxHalf_Z)?(boxSize_Z-xtmp):xtmp)
#else
#define NGB_PERIODIC_LONG_X(x) fabs(x)
#define NGB_PERIODIC_LONG_Y(x) fabs(x)
#define NGB_PERIODIC_LONG_Z(x) fabs(x)
#endif

#define FACT1 0.366025403785	/* FACT1 = 0.5 * (sqrt(3)-1) */

/*********************************************************/
/*  Global variables                                     */
/*********************************************************/

extern int FirstActiveParticle;
extern int *NextActiveParticle;
extern unsigned char *ProcessedFlag;

extern int TimeBinCount[TIMEBINS];
extern int TimeBinCountSph[TIMEBINS];
extern int TimeBinActive[TIMEBINS];

extern int FirstInTimeBin[TIMEBINS];
extern int LastInTimeBin[TIMEBINS];
extern int *NextInTimeBin;
extern int *PrevInTimeBin;

#ifdef SFR
extern double TimeBinSfr[TIMEBINS];
#endif

#ifdef BLACK_HOLES
extern double TimeBin_BH_mass[TIMEBINS];
extern double TimeBin_BH_dynamicalmass[TIMEBINS];
extern double TimeBin_BH_Mdot[TIMEBINS];
extern double TimeBin_BH_Medd[TIMEBINS];
#endif

extern int ThisTask;		/*!< the number of the local processor  */
extern int NTask;		/*!< number of processors */
extern int PTask;		/*!< note: NTask = 2^PTask */

#ifdef INVARIANCETEST
extern int World_ThisTask;
extern int World_NTask;
extern int Color;
extern MPI_Comm MPI_CommLocal;
#ifndef DO_NOT_REDEFINE_MPI_COMM_WORLD
#undef  MPI_COMM_WORLD
#define MPI_COMM_WORLD MPI_CommLocal
#endif
#endif

extern double CPUThisRun;	/*!< Sums CPU time of current process */

extern int NumForceUpdate;	/*!< number of active particles on local processor in current timestep  */
extern long long GlobNumForceUpdate;

extern int NumSphUpdate;	/*!< number of active SPH particles on local processor in current timestep  */

extern int MaxTopNodes;	        /*!< Maximum number of nodes in the top-level tree used for domain decomposition */

extern int RestartFlag;		/*!< taken from command line used to start code. 0 is normal start-up from
				   initial conditions, 1 is resuming a run from a set of restart files, while 2
				   marks a restart from a snapshot file. */
extern int RestartSnapNum;

extern int *Exportflag;	        /*!< Buffer used for flagging whether a particle needs to be exported to another process */
extern int *Exportnodecount;
extern int *Exportindex;

extern int *Send_offset, *Send_count, *Recv_count, *Recv_offset, *Sendcount_matrix;

extern size_t AllocatedBytes;
extern size_t FreeBytes;

extern double CPU_Step[CPU_PARTS];
extern char CPU_Symbol[CPU_PARTS];
extern char CPU_SymbolImbalance[CPU_PARTS];
extern char CPU_String[CPU_STRING_LEN + 1];

extern double WallclockTime;    /*!< This holds the last wallclock time measurement for timings measurements */

extern int Flag_FullStep;	/*!< Flag used to signal that the current step involves all particles */

extern int TreeReconstructFlag;
extern int GlobFlag;

extern char DumpFlag;

extern int NumPart;		/*!< number of particles on the LOCAL processor */
extern int N_gas;		/*!< number of gas particles on the LOCAL processor  */
extern long long Ntype[6];	/*!< total number of particles of each type */
extern int NtypeLocal[6];	/*!< local number of particles of each type */

extern gsl_rng *random_generator;	/*!< the random number generator used */

#ifdef SFR
extern int Stars_converted;	/*!< current number of star particles in gas particle block */
#endif

extern double TimeOfLastTreeConstruction;	/*!< holds what it says */

extern int *Ngblist;		/*!< Buffer to hold indices of neighbours retrieved by the neighbour search
				   routines */

extern double *R2ngblist;

extern double DomainCorner[3], DomainCenter[3], DomainLen, DomainFac;
extern int *DomainStartList, *DomainEndList;

extern double *DomainWork;
extern int *DomainCount;
extern int *DomainCountSph;
extern int *DomainTask;
extern int *DomainNodeIndex;
extern int *DomainList, DomainNumChanged;

extern peanokey *Key, *KeySorted;

extern struct topnode_data
{
  peanokey Size;
  peanokey StartKey;
  long long Count;
  MyFloat GravCost;
  int Daughter;
  int Pstart;
  int Blocks;
  int Leaf;
} *TopNodes;

extern int NTopnodes, NTopleaves;

extern double RndTable[RNDTABLE];

#ifdef SUBFIND
extern int GrNr;
extern int NumPartGroup;
#endif

/* variables for input/output , usually only used on process 0 */

extern char ParameterFile[100];	/*!< file name of parameterfile used for starting the simulation */

extern FILE *FdInfo,		/*!< file handle for info.txt log-file. */
 *FdEnergy,			/*!< file handle for energy.txt log-file. */
 *FdTimings,			/*!< file handle for timings.txt log-file. */
 *FdBalance,			/*!< file handle for balance.txt log-file. */
 *FdCPU;			/*!< file handle for cpu.txt log-file. */

#ifdef SFR
extern FILE *FdSfr;		/*!< file handle for sfr.txt log-file. */
#endif

#ifdef RADTRANSFER
extern FILE *FdEddington;     /*!< file handle for eddington.txt log-file. */
extern FILE *FdRadtransfer;		/*!< file handle for radtransfer.txt log-file. */
#endif

#ifdef DISTORTIONTENSORPS
#ifdef PMGRID
extern FILE *FdTidaltensor;     /*!< file handle for tidaltensor.txt log-file. */
#endif
extern FILE *FdCaustics;	/*!< file handle for Caustics.txt log-file. */
#endif

#ifdef BLACK_HOLES
extern FILE *FdBlackHoles;	/*!< file handle for blackholes.txt log-file. */
extern FILE *FdBlackHolesDetails;
#endif

#ifdef FORCETEST
extern FILE *FdForceTest;	/*!< file handle for forcetest.txt log-file. */
#endif

#ifdef DARKENERGY
extern FILE *FdDE;  /*!< file handle for darkenergy.txt log-file. */
#endif

#ifdef XXLINFO
extern FILE *FdXXL;		/*!< file handle for xxl.txt log-file. */

#ifdef MAGNETIC
extern double MeanB;

#ifdef TRACEDIVB
extern double MaxDivB;
#endif
#endif
#ifdef TIME_DEP_ART_VISC
extern double MeanAlpha;
#endif
#endif

/*! table for the cosmological drift factors */
extern double DriftTable[DRIFT_TABLE_LENGTH];

/*! table for the cosmological kick factor for gravitational forces */
extern double GravKickTable[DRIFT_TABLE_LENGTH];

/*! table for the cosmological kick factor for hydrodynmical forces */
extern double HydroKickTable[DRIFT_TABLE_LENGTH];

extern void *CommBuffer;	/*!< points to communication buffer, which is used at a few places */

/*! This structure contains data which is the SAME for all tasks (mostly code parameters read from the
 * parameter file).  Holding this data in a structure is convenient for writing/reading the restart file, and
 * it allows the introduction of new global variables in a simple way. The only thing to do is to introduce
 * them into this structure.
 */
extern struct global_data_all_processes
{
  long long TotNumPart;		/*!<  total particle numbers (global value) */
  long long TotN_gas;		/*!<  total gas particle number (global value) */

#ifdef BLACK_HOLES
  int TotBHs;
#endif

  int MaxPart;			/*!< This gives the maxmimum number of particles that can be stored on one
				   processor. */
  int MaxPartSph;		/*!< This gives the maxmimum number of SPH particles that can be stored on one
				   processor. */

  int ICFormat;			/*!< selects different versions of IC file-format */

  int SnapFormat;		/*!< selects different versions of snapshot file-formats */

  int DoDynamicUpdate;

  int NumFilesPerSnapshot;	/*!< number of files in multi-file snapshot dumps */
  int NumFilesWrittenInParallel;	/*!< maximum number of files that may be written simultaneously when
					   writing/reading restart-files, or when writing snapshot files */

  int BufferSize;		/*!< size of communication buffer in MB */
  int BunchSize;     	        /*!< number of particles fitting into the buffer in the parallel tree algorithm  */

  double PartAllocFactor;	/*!< in order to maintain work-load balance, the particle load will usually
				   NOT be balanced.  Each processor allocates memory for PartAllocFactor times
				   the average number of particles to allow for that */

  double TreeAllocFactor;	/*!< Each processor allocates a number of nodes which is TreeAllocFactor times
				   the maximum(!) number of particles.  Note: A typical local tree for N
				   particles needs usually about ~0.65*N nodes. */

  double TopNodeAllocFactor;	/*!< Each processor allocates a number of nodes which is TreeAllocFactor times
				   the maximum(!) number of particles.  Note: A typical local tree for N
				   particles needs usually about ~0.65*N nodes. */

#ifdef SCALARFIELD
  double ScalarBeta;
  double ScalarScreeningLength;
#endif

  /* some SPH parameters */

  int DesNumNgb;		/*!< Desired number of SPH neighbours */
#ifdef SUBFIND
  int DesLinkNgb;
  double ErrTolThetaSubfind;
#endif

  double MaxNumNgbDeviation;	/*!< Maximum allowed deviation neighbour number */
#ifdef START_WITH_EXTRA_NGBDEV
  double MaxNumNgbDeviationStart;    /*!< Maximum allowed deviation neighbour number to start with*/
#endif

  double ArtBulkViscConst;	/*!< Sets the parameter \f$\alpha\f$ of the artificial viscosity */
  double InitGasTemp;		/*!< may be used to set the temperature in the IC's */
  double InitGasU;		/*!< the same, but converted to thermal energy per unit mass */
  double MinGasTemp;		/*!< may be used to set a floor for the gas temperature */
  double MinEgySpec;		/*!< the minimum allowed temperature expressed as energy per unit mass */

  /* some force counters  */

  long long TotNumOfForces;	/*!< counts total number of force computations  */

  long long NumForcesSinceLastDomainDecomp;	/*!< count particle updates since last domain decomposition */

  /* some variable for dynamic work-load adjustment based on CPU measurements */

  double Cadj_Cost;
  double Cadj_Cpu;

  /* system of units  */

  double UnitTime_in_s,		/*!< factor to convert internal time unit to seconds/h */
    UnitMass_in_g,		/*!< factor to convert internal mass unit to grams/h */
    UnitVelocity_in_cm_per_s,	/*!< factor to convert intqernal velocity unit to cm/sec */
    UnitLength_in_cm,		/*!< factor to convert internal length unit to cm/h */
    UnitPressure_in_cgs,	/*!< factor to convert internal pressure unit to cgs units (little 'h' still
				   around!) */
    UnitDensity_in_cgs,		/*!< factor to convert internal length unit to g/cm^3*h^2 */
    UnitCoolingRate_in_cgs,	/*!< factor to convert internal cooling rate to cgs units */
    UnitEnergy_in_cgs,		/*!< factor to convert internal energy to cgs units */
    UnitTime_in_Megayears,	/*!< factor to convert internal time to megayears/h */
    GravityConstantInternal,	/*!< If set to zero in the parameterfile, the internal value of the
				   gravitational constant is set to the Newtonian value based on the system of
				   units specified. Otherwise the value provided is taken as internal gravity
				   constant G. */
    G;				/*!< Gravity-constant in internal units */
#ifdef ANNIHILATION_RADIATION
  double UnitDensity_in_Gev_per_cm3; /*!< factor to convert internal density unit to GeV/c^2 / cm^3 */
#endif
  /* Cosmology */

  double Hubble;		/*!< Hubble-constant in internal units */
  double Omega0,		/*!< matter density in units of the critical density (at z=0) */
    OmegaLambda,		/*!< vaccum energy density relative to crictical density (at z=0) */
    OmegaBaryon,		/*!< baryon density in units of the critical density (at z=0) */
    HubbleParam;		/*!< little `h', i.e. Hubble constant in units of 100 km/s/Mpc.  Only needed to get absolute
				 * physical values for cooling physics
				 */

  double BoxSize;		/*!< Boxsize in case periodic boundary conditions are used */

  /* Code options */

  int ComovingIntegrationOn;	/*!< flags that comoving integration is enabled */
  int PeriodicBoundariesOn;	/*!< flags that periodic boundaries are enabled */
  int ResubmitOn;		/*!< flags that automatic resubmission of job to queue system is enabled */
  int TypeOfOpeningCriterion;	/*!< determines tree cell-opening criterion: 0 for Barnes-Hut, 1 for relative
				   criterion */
  int TypeOfTimestepCriterion;	/*!< gives type of timestep criterion (only 0 supported right now - unlike
				   gadget-1.1) */
  int OutputListOn;		/*!< flags that output times are listed in a specified file */
  int CoolingOn;		/*!< flags that cooling is enabled */
  int StarformationOn;		/*!< flags that star formation is enabled */

  /* parameters determining output frequency */

  int SnapshotFileCount;	/*!< number of snapshot that is written next */
  double TimeBetSnapshot,	/*!< simulation time interval between snapshot files */
    TimeOfFirstSnapshot,	/*!< simulation time of first snapshot files */
    CpuTimeBetRestartFile,	/*!< cpu-time between regularly generated restart files */
    TimeLastRestartFile,	/*!< cpu-time when last restart-file was written */
    TimeBetStatistics,		/*!< simulation time interval between computations of energy statistics */
    TimeLastStatistics;		/*!< simulation time when the energy statistics was computed the last time */
  int NumCurrentTiStep;		/*!< counts the number of system steps taken up to this point */

  /* Current time of the simulation, global step, and end of simulation */

  double Time,			/*!< current time of the simulation */
    TimeBegin,			/*!< time of initial conditions of the simulation */
    TimeStep,			/*!< difference between current times of previous and current timestep */
    TimeMax;			/*!< marks the point of time until the simulation is to be evolved */

  /* variables for organizing discrete timeline */

  double Timebase_interval;	/*!< factor to convert from floating point time interval to integer timeline */
  int Ti_Current;		/*!< current time on integer timeline */
  int Ti_nextoutput;		/*!< next output time on integer timeline */

#ifdef PMGRID
  int PM_Ti_endstep, PM_Ti_begstep;
  double Asmth[2], Rcut[2];
  double Corner[2][3], UpperCorner[2][3], Xmintot[2][3], Xmaxtot[2][3];
  double TotalMeshSize[2];
#endif

#ifdef CHEMISTRY
  double Epsilon;
#endif

  int Ti_nextlineofsight;
#ifdef OUTPUTLINEOFSIGHT
  double TimeFirstLineOfSight;
#endif

  /* variables that keep track of cumulative CPU consumption */

  double TimeLimitCPU;
  double CPU_Sum[CPU_PARTS];    /*!< sums wallclock time/CPU consumption in whole run */

  /* tree code opening criterion */

  double ErrTolTheta;		/*!< BH tree opening angle */
  double ErrTolForceAcc;	/*!< parameter for relative opening criterion in tree walk */

  /* adjusts accuracy of time-integration */

  double ErrTolIntAccuracy;	/*!< accuracy tolerance parameter \f$ \eta \f$ for timestep criterion. The
				   timesteps is \f$ \Delta t = \sqrt{\frac{2 \eta eps}{a}} \f$ */

  double MinSizeTimestep,	/*!< minimum allowed timestep. Normally, the simulation terminates if the
				   timestep determined by the timestep criteria falls below this limit. */
    MaxSizeTimestep;		/*!< maximum allowed timestep */

  double MaxRMSDisplacementFac;	/*!< this determines a global timestep criterion for cosmological simulations
				   in comoving coordinates.  To this end, the code computes the rms velocity
				   of all particles, and limits the timestep such that the rms displacement
				   is a fraction of the mean particle separation (determined from the
				   particle mass and the cosmological parameters). This parameter specifies
				   this fraction. */

  double CourantFac;		/*!< SPH-Courant factor */

  /* frequency of tree reconstruction/domain decomposition */

  double TreeDomainUpdateFrequency;	/*!< controls frequency of domain decompositions  */

  /* gravitational and hydrodynamical softening lengths (given in terms of an `equivalent' Plummer softening
   * length)
   *
   * five groups of particles are supported 0=gas,1=halo,2=disk,3=bulge,4=stars
   */
  double MinGasHsmlFractional,	/*!< minimum allowed SPH smoothing length in units of SPH gravitational
				   softening length */
    MinGasHsml;			/*!< minimum allowed SPH smoothing length */

  double SofteningGas,		/*!< for type 0 */
    SofteningHalo,		/*!< for type 1 */
    SofteningDisk,		/*!< for type 2 */
    SofteningBulge,		/*!< for type 3 */
    SofteningStars,		/*!< for type 4 */
    SofteningBndry;		/*!< for type 5 */

  double SofteningGasMaxPhys,	/*!< for type 0 */
    SofteningHaloMaxPhys,	/*!< for type 1 */
    SofteningDiskMaxPhys,	/*!< for type 2 */
    SofteningBulgeMaxPhys,	/*!< for type 3 */
    SofteningStarsMaxPhys,	/*!< for type 4 */
    SofteningBndryMaxPhys;	/*!< for type 5 */

  double SofteningTable[6];	/*!< current (comoving) gravitational softening lengths for each particle type */
  double ForceSoftening[6];	/*!< the same, but multiplied by a factor 2.8 - at that scale the force is Newtonian */

  /*! If particle masses are all equal for one type, the corresponding entry in MassTable is set to this
   *  value, * allowing the size of the snapshot files to be reduced
   */
  double MassTable[6];

  /* some filenames */
  char InitCondFile[100],
    OutputDir[100],
    SnapshotFileBase[100],
    EnergyFile[100],
    CpuFile[100],
    InfoFile[100], TimingsFile[100], RestartFile[100], ResubmitCommand[100], OutputListFilename[100];

  /*! table with desired output times */
  double OutputListTimes[MAXLEN_OUTPUTLIST];
  char OutputListFlag[MAXLEN_OUTPUTLIST];
  int OutputListLength;		/*!< number of times stored in table of desired output times */

#if defined(ADAPTIVE_GRAVSOFT_FORGAS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS_HSML)
  double ReferenceGasMass;
#endif

#ifdef VORONOI_MESHRELAX
  double MeanMass;
  double MeanPressure;
#endif

#ifdef RADTRANSFER
  double IonizingLumPerSolarMass;
  double IonizingLumPerSFR;
  int Radiation_Ti_begstep;
  int Radiation_Ti_endstep;
#endif

#if defined(SIM_ADAPTIVE_SOFT) || defined(REINIT_AT_TURNAROUND)
  double CurrentTurnaroundRadius;
  double SIM_epsilon;
#endif

#ifdef ANNIHILATION_RADIATION
  /* present day velocity dispersion of DM particle in cm/s (e.g. Neutralino = 0.03 cm/s) */
  double DM_velocity_dispersion;
#endif

#ifdef SFR		/* star formation and feedback sector */
  double CritOverDensity;
  double CritPhysDensity;
  double OverDensThresh;
  double PhysDensThresh;
  double EgySpecSN;
  double FactorSN;
  double EgySpecCold;
  double FactorEVP;
  double FeedbackEnergy;
  double TempSupernova;
  double TempClouds;
  double MaxSfrTimescale;
  double WindEfficiency;
  double WindEnergyFraction;
  double WindFreeTravelLength;
  double WindFreeTravelDensFac;
  double FactorForSofterEQS;
#endif

#ifdef CS_MODEL
  double FactorSFR;
  double DecouplingParam;
  double MinTlifeSNI;
  double MaxTlifeSNI;
  double TlifeSNII;
  int    Raiteri_TlifeSNII;
  double RateSNI;
  double SN_Energy_cgs;
  double Tcrit_Phase;
  double DensFrac_Phase;
  double SN_Energy_frac_cold;
  double MaxHotHsmlParam;
  double InitialHotHsmlFactor;
  double DensityTailThreshold;
  double MaxNumHotNgbDeviation;	/*!< Maximum allowed deviation HOT neighbour number */
#endif

#ifdef DARKENERGY
  double DarkEnergyParam;	/*!< fixed w for equation of state */
#ifdef TIMEDEPDE
  char DarkEnergyFile[100];	/*!< tabelized w for equation of state */
#ifdef TIMEDEPGRAV
  double Gini;
#endif
#endif
#endif

#ifdef RESCALEVINI
  double VelIniScale;		/*!< Scale the initial velocities by this amount */
#endif

#if defined(SNIA_HEATING)
  double SnIaHeatingRate;
#endif

#ifdef TIME_DEP_ART_VISC
  double ViscSource0;		/*!< Given sourceterm in viscosity evolution */
  double DecayLength;		/*!< Number of h for the viscosity decay */
  double ViscSource;		/*!< Reduced sourceterm in viscosity evolution */
  double DecayTime;		/*!< Calculated decaytimescale */
  double AlphaMin;		/*!< Minimum of allowed viscosity parameter */
#endif

#ifdef CONDUCTION
  double ConductionCoeff;	/*!< Thermal Conductivity */
#ifdef CONDUCTION_SATURATION
  double ElectronFreePathFactor;	/*!< Factor to get electron mean free path */
#endif

  int Conduction_Ti_begstep, Conduction_Ti_endstep;
  double MaxSizeConductionStep;
#endif

#ifdef MAGNETIC
#ifdef BINISET
  double BiniX, BiniY, BiniZ;	/*!< Initial values for B */
#endif

#ifdef BSMOOTH
  int BSmoothInt;
  double BSmoothFrac;
  int MainTimestepCounts;
#ifdef SETMAINTIMESTEPCOUNT
  int MainTimestepCountIni;
#endif
#endif

#ifdef MAGNETIC_DISSIPATION
  double ArtMagDispConst;	/*!< Sets the parameter \f$\alpha\f$ of the artificial magnetic disipation */
#ifdef TIME_DEP_MAGN_DISP
  double ArtMagDispMin;
  double ArtMagDispSource;
  double ArtMagDispTime;
#endif
#endif

#ifdef DIVBCLEANING_DEDNER
  double DivBcleanParabolicSigma;
  double DivBcleanHyperbolicSigma;
#endif

#ifdef HEALPIX
  //change this to read in the Parameterfile
  int Nside;
#define NSIDE2NPIX(nside)  (12*nside*nside)
  float *healpixmap;
#endif

#ifdef MAGNETIC_DIFFUSION
  double MagneticEta;
#endif
#endif

#ifdef BLACK_HOLES
  double TimeNextBlackHoleCheck;
  double TimeBetBlackHoleSearch;
  double BlackHoleAccretionFactor;	/*!< Fraction of BH bondi accretion rate */
  double BlackHoleFeedbackFactor;	/*!< Fraction of the black luminosity feed into thermal feedback */
  double SeedBlackHoleMass;	/*!< Seed black hole mass */
  double MinFoFMassForNewSeed;	/*!< Halo mass required before new seed is put in */
  double BlackHoleNgbFactor;	/*!< Factor by which the normal SPH neighbour should be increased/decreased */
  double BlackHoleMaxAccretionRadius;
  double BlackHoleEddingtonFactor;	/*! Factor above Eddington */
#ifdef FOF
  double massDMpart;
#endif
#ifdef MODIFIEDBONDI
  double BlackHoleRefDensity;
  double BlackHoleRefSoundspeed;
#endif
#endif

#ifdef COSMIC_RAYS
  double CR_Alpha[NUMCRPOP];	/*!< Cosmic ray spectral index [2..3] */
  double CR_SNEff;		/*!< SN injection efficiency [0..1] */
  double CR_SNAlpha;		/*!< SN injection spectral index [2..3] */
  int bDebugFlag;		/*!< enables debug outputs after triggered */

#if defined(CR_DIFFUSION)
  double CR_DiffusionCoeff;	/*!< (temporary) fixed value for CR diffusivity */

  double CR_DiffusionDensScaling;	/*!< grade of density dependence of diffusivity */
  double CR_DiffusionDensZero;	/*!< Reference point density for diffusivity */

  double CR_DiffusionEntropyScaling;	/*!< grade of specific energy dependence of diffusivity */

  double CR_DiffusionEntropyZero;	/*!< Reference Entropic function for diffusivity */

  double CR_DiffusionMaxSizeTimestep;
  int CR_Diffusion_Ti_begstep, CR_Diffusion_Ti_endstep;
#endif				/* CR_DIFFUSION */

#if defined(CR_SHOCK)
#if (CR_SHOCK == 1)
  double CR_ShockAlpha;		/*!< spectral index to be used in shock injection */
#else
  double CR_ShockCutoff;	/*!< Cutoff factor x_inj for CR accel */
#endif
  double CR_ShockEfficiency;	/*!< energy fraction of shock energy fed into CR */
#endif				/* CR_SHOCK */

#ifdef FIX_QINJ
  double Shock_Fix_Qinj;	/*!< inject only CRps with threshold cutoff Shock_Fix_Qinj */
#endif

#ifdef CR_BUBBLES
  double CR_AGNEff;               /*!< AGN injection efficiency [0..1] */
#endif
#endif				/* COSMIC_RAYS */

#ifdef MACHNUM
  double Shock_Length;		/*!< length scale on which the shock is smoothed out */
  double Shock_DeltaDecayTimeMax;	/*!< maximum time interval (Dloga) for which the
					   Mach number is kept at its maximum */
#endif

#ifdef REIONIZATION
  int not_yet_reionized;	/*!< flag that makes sure that there is only one reionization */
#endif

#ifdef BUBBLES
  double BubbleDistance;
  double BubbleRadius;
  double BubbleTimeInterval;
  double BubbleEnergy;
  double TimeOfNextBubble;
  double FirstBubbleRedshift;
#ifdef FOF
  int BiggestGroupLen;
  float BiggestGroupCM[3];
  double BiggestGroupMass;
#endif
#endif

#ifdef BH_BUBBLES
  double BubbleDistance;
  double BubbleRadius;
  double BubbleEnergy;
  double BlackHoleRadioTriggeringFactor;
  double DefaultICMDensity;
  double RadioFeedbackFactor;
#ifdef UNIFIED_FEEDBACK
  double RadioThreshold;
#endif
#endif

#if defined(MULTI_BUBBLES) && defined(FOF)
#ifndef BLACK_HOLES
  double MinFoFMassForNewSeed;	/*!< Halo mass required before new seed is put in */
  double massDMpart;
#endif
  double BubbleDistance;
  double BubbleRadius;
  double BubbleTimeInterval;
  double BubbleEnergy;
  double TimeOfNextBubble;
  double ClusterMass200;
  double FirstBubbleRedshift;
#endif

#ifdef NAVIERSTOKES
  double NavierStokes_ShearViscosity;
  double FractionSpitzerViscosity;
  double ShearViscosityTemperature;
#endif
#ifdef NAVIERSTOKES_BULK
  double NavierStokes_BulkViscosity;
#endif
#ifdef VISCOSITY_SATURATION
  double IonMeanFreePath;
#endif

#ifdef EOS_DEGENERATE
  char EosTable[100];
#endif
}
All;

/*! This structure holds all the information that is
 * stored for each particle of the simulation.
 */
extern struct particle_data
{
  MyDouble Pos[3];   /*!< particle position at its current time */
  MyDouble Vel[3];   /*!< particle velocity at its current time */
  MyDouble Mass;     /*!< particle mass */
  MyIDType ID;

  union
  {
    MyFloat       GravAccel[3];		/*!< particle acceleration due to gravity */
    MyLongDouble dGravAccel[3];
  } g;
#ifdef PMGRID
  MyFloat GravPM[3];		/*!< particle acceleration due to long-range PM gravity force */
#endif
#ifdef FORCETEST
  MyFloat GravAccelDirect[3];	/*!< particle acceleration calculated by direct summation */
#endif
#if defined(EVALPOTENTIAL) || defined(COMPUTE_POTENTIAL_ENERGY) || defined(OUTPUTPOTENTIAL)
  union
  {
    MyFloat       Potential;		/*!< gravitational potential */
    MyLongDouble dPotential;
  } p;
#endif

#ifdef DISTORTIONTENSORPS
  MyLongDouble distortion_tensorps[6][6];   /*!< Phase Space Distortion tensor */
  MyLongDouble tidal_tensorps[3][3];        /*!< tidal tensor (=second derivatives of grav. potential) */
  MyLongDouble V_matrix[3][3];              /*!< initial orientation of CDM sheet the particle is embedded in */
  MyDouble init_density;                    /*!< initial stream density */
  MyFloat caustic_counter;                  /*!< caustic counter */
  MyDouble last_stream_determinant;         /*!< last stream density determinant, needed to identify caustics */
#ifdef REINIT_AT_TURNAROUND
  int turnaround_flag;                      /*!< mark when a particle turned around */
#endif
#ifdef ANNIHILATION_RADIATION
 MyDouble annihilation;                     /*!< integrated annihilation rate */
 MyDouble rho_normed_cutoff_current;        /*!< current and last normed_cutoff density in rho_max/rho_init * sqrt(sigma) */
 MyDouble rho_normed_cutoff_last;
 MyDouble stream_density;                   /*!< physical stream density that is going to be integrated (in terms of rho_crit) */
 MyFloat analytic_caustics;                 /*!< number of caustics that were integrated analyticall,
                                                 i.e. where the physical caustic density was higher
                                                 than the numerical GDE density */
#endif
#ifdef OUTPUT_LAST_CAUSTIC
 MyDouble lc_Time;                          /*!< time of caustic passage */
 MyDouble lc_Pos[3];                        /*!< position of caustic */
 MyDouble lc_Vel[3];                        /*!< particle velocity when passing through caustic */
 MyDouble lc_rho_normed_cutoff;             /*!< normed_cutoff density at caustic */
 MyDouble lc_Dir_x[3];                      /*!< principal axis frame of smear out */
 MyDouble lc_Dir_y[3];
 MyDouble lc_Dir_z[3];
 MyDouble lc_smear_x;                       /*!< smear out length */
 MyDouble lc_smear_y;
 MyDouble lc_smear_z;
#endif
#ifdef PMGRID
  MyLongDouble tidal_tensorpsPM[3][3];	    /*!< for TreePM simulations, long range tidal field */
#endif
#endif

  MyFloat OldAcc;			/*!< magnitude of old gravitational force. Used in relative opening
                                          criterion */
#if defined(EVALPOTENTIAL) && defined(PMGRID)
  MyFloat PM_Potential;
#endif

#ifdef STELLARAGE
  MyFloat StellarAge;		/*!< formation time of star particle */
#endif
#ifdef METALS
  MyFloat Metallicity;		/*!< metallicity of gas or star particle */
#endif				/* closes METALS */

#if defined (BLACK_HOLES) || defined(CS_MODEL) || defined(RADTRANSFER) || defined(SNIA_HEATING)
  MyFloat Hsml;

  union
  {
    MyFloat       NumNgb;
    MyLongDouble dNumNgb;
  } n;
#if defined(RADTRANSFER) || defined(SNIA_HEATING)
  MyFloat DensAroundStar;
#endif
#endif

#ifdef BLACK_HOLES
  int SwallowID;
#ifdef BH_COUNTPROGS
  int BH_CountProgs;
#endif
  MyFloat BH_Mass;
  MyFloat BH_Mdot;
#ifdef BH_BUBBLES
  MyFloat BH_Mass_bubbles;
  MyFloat BH_Mass_ini;
#ifdef UNIFIED_FEEDBACK
  MyFloat BH_Mass_radio;
#endif
#endif
  union
  {
    MyFloat BH_Density;
    MyLongDouble dBH_Density;
  } b1;
  union
  {
    MyFloat BH_Entropy;
    MyLongDouble dBH_Entropy;
  } b2;
  union
  {
    MyFloat BH_SurroundingGasVel[3];
    MyLongDouble dBH_SurroundingGasVel[3];
  } b3;
  union
  {
    MyFloat BH_accreted_Mass;
    MyLongDouble dBH_accreted_Mass;
  } b4;
  union
  {
    MyFloat BH_accreted_BHMass;
    MyLongDouble dBH_accreted_BHMass;
  } b5;
  union
  {
    MyFloat BH_accreted_momentum[3];
    MyLongDouble dBH_accreted_momentum[3];
  } b6;
#ifdef BH_BUBBLES
  union
  {
    MyFloat BH_accreted_BHMass_bubbles;
    MyLongDouble dBH_accreted_BHMass_bubbles;
  } b7;
#ifdef UNIFIED_FEEDBACK
  union
  {
    MyFloat BH_accreted_BHMass_radio;
    MyLongDouble dBH_accreted_BHMass_radio;
  } b8;
#endif
#endif
#ifdef REPOSITION_ON_POTMIN
  MyFloat BH_MinPotPos[3];
  MyFloat BH_MinPot;
#endif
#ifdef BH_KINETICFEEDBACK
  MyFloat ActiveTime;
  MyFloat ActiveEnergy;
#endif
#endif

#ifdef SUBFIND
  unsigned int GrNr;
  int SubNr;
  int DM_NumNgb;
  int targettask, origintask, submark, origindex;
  MyFloat DM_Hsml;
  union
  {
    MyFloat DM_Density;
    MyFloat DM_Potential;
  } u;
  union
  {
    MyFloat DM_VelDisp;
    MyFloat DM_BindingEnergy;
  } v;
#ifdef DENSITY_SPLIT_BY_TYPE
  union
  {
    MyFloat int_energy;
    MyFloat density_sum;
  } w;
#endif

#ifdef SUBFIND_RESHUFFLE_CATALOGUE_WITH_VORONOI
  MyFloat DM_Hsml_V;
  MyFloat DM_Density_V;
#endif

#ifdef SAVE_HSML_IN_IC_ORDER
  MyIDType ID_ic_order;
#endif
#ifdef SUBFIND_ALTERNATIVE_COLLECTIVE
  peanokey Key;
#endif
#endif

#if defined(ORDER_SNAPSHOTS_BY_ID) && !defined(SUBFIND)
  int     GrNr;
  int     SubNr;
#endif

#ifdef SHELL_CODE
  MyDouble radius;
  MyDouble enclosed_mass;
  MyDouble dMdr;
#endif

#ifdef CS_MODEL
  MyFloat Zm[12];
  MyFloat ZmReservoir[12];
#ifdef CS_FEEDBACK
  MyFloat EnergySN;
  MyFloat EnergySNCold;
#endif
#endif

  float GravCost;		/*!< weight factor used for balancing the work-load */

  int Ti_begstep;		/*!< marks start of current timestep of particle on integer timeline */
  int Ti_current;		/*!< current time of the particle */

  short int Type;		/*!< flags particle type.  0=gas, 1=halo, 2=disk, 3=bulge, 4=stars, 5=bndry */
  short int TimeBin;

#ifdef WAKEUP
  int dt_step;
#endif
}
 *P,				/*!< holds particle data on local processor */
 *DomainPartBuf;		/*!< buffer for particle data used in domain decomposition */

/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
 * variables.
 */
extern struct sph_particle_data
{
  MyDouble Entropy;		/*!< current value of entropy (actually entropic function) of particle */
  MyFloat  Pressure;		/*!< current pressure */
  MyFloat  VelPred[3];		/*!< predicted SPH particle velocity at the current time */
#ifdef ALTERNATIVE_VISCOUS_TIMESTEP
  MyFloat MinViscousDt;
#else
  MyFloat MaxSignalVel;           /*!< maximum signal velocity */
#endif
#ifdef VOLUME_CORRECTION
  MyFloat DensityOld;
  MyFloat DensityStd;
#endif

#ifdef VORONOI
  MyFloat MaxDelaunayRadius;
  MyFloat Volume;
  MyFloat Center[3];
#endif

  union
  {
    MyFloat       Density;		/*!< current baryonic mass density of particle */
    MyLongDouble dDensity;
  } d;
  union
  {
    MyFloat       DtEntropy;		/*!< rate of change of entropy */
    MyLongDouble dDtEntropy;
  } e;
  union
  {
    MyFloat       HydroAccel[3];	/*!< acceleration due to hydrodynamical force */
    MyLongDouble dHydroAccel[3];
  } a;
  union
  {
    MyFloat       DhsmlDensityFactor;	/*!< correction factor needed in entropy formulation of SPH */
    MyLongDouble dDhsmlDensityFactor;
  } h;
  union
  {
    MyFloat       DivVel;		/*!< local velocity divergence */
    MyLongDouble dDivVel;
  } v;
#ifndef NAVIERSTOKES
  union
  {
    MyFloat CurlVel;     	        /*!< local velocity curl */
    MyFloat       Rot[3];		/*!< local velocity curl */
    MyLongDouble dRot[3];
  } r;
#else
  union
  {
    MyFloat DV[3][3];
    struct
    {
      MyFloat DivVel;
      MyFloat CurlVel;
      MyFloat StressDiag[3];
      MyFloat StressOffDiag[3];
#ifdef NAVIERSTOKES_BULK
      MyFloat StressBulk;
#endif
    } s;
  } u;
#endif

#if !(defined(BLACK_HOLES) || defined(CS_MODEL) || defined(RADTRANSFER) || defined(SNIA_HEATING))
  MyFloat Hsml;			/*!< current smoothing length */
  union
  {
    MyFloat       NumNgb;
    MyLongDouble dNumNgb;
  } n;
#endif

#if defined(BH_THERMALFEEDBACK) || defined(BH_KINETICFEEDBACK)
  union
  {
    MyFloat       Injected_BH_Energy;
    MyLongDouble dInjected_BH_Energy;
  } i;
#endif

#ifdef COOLING
  MyFloat Ne;  /*!< electron fraction, expressed as local electron number
		    density normalized to the hydrogen number density. Gives
		    indirectly ionization state and mean molecular weight. */
#endif
#ifdef SFR
  MyFloat Sfr;
#endif
#ifdef WINDS
  MyFloat DelayTime;		/*!< remaining maximum decoupling time of wind particle */
#endif

#ifdef MAGNETIC
  MyFloat BPred[3];
#ifdef EULERPOTENTIALS
  MyFloat EulerA,EulerB;
  MyFloat dEulerA[3],dEulerB[3];
#else
  MyFloat B[3];
  MyFloat DtB[3];
#endif
#if defined(TRACEDIVB) || defined(TIME_DEP_MAGN_DISP)
  MyFloat divB;
#endif
#if defined(BSMOOTH) || defined(BFROMROTA)
  MyFloat BSmooth[3];
#endif
#ifdef TIME_DEP_MAGN_DISP
  MyFloat Balpha, DtBalpha;
#endif
#ifdef DIVBCLEANING_DEDNER
  MyFloat Phi, PhiPred, DtPhi;
#ifdef SMOOTH_PHI
  MyFloat SmoothPhi;
  MyFloat SmoothDivB;
#endif
#endif
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
  MyFloat RotB[3];
#ifdef SMOOTH_ROTB
  MyFloat SmoothedRotB[3];
#endif
#endif
#if defined(BSMOOTH) || defined(SMOOTH_ROTB) || defined(SMOOTH_PHI)
  MyFloat DensityNorm;
#endif
#endif
#ifdef TIME_DEP_ART_VISC
  MyFloat alpha, Dtalpha;
#endif
#ifdef NS_TIMESTEP
  MyFloat ViscEntropyChange;
#endif
#ifdef CONDUCTION_SATURATION
  MyFloat GradEntr[3];
#endif

#ifdef MHM
  MyFloat FeedbackEnergy;
#endif

#ifdef COSMIC_RAYS
  MyFloat CR_C0[NUMCRPOP];			/*!< Cosmic ray amplitude adiabatic invariable */
  MyFloat CR_q0[NUMCRPOP];			/*!< Cosmic ray cutoff adiabatic invariable */
  MyFloat CR_E0[NUMCRPOP];			/*!< Specific Energy at Rho0 */
  MyFloat CR_n0[NUMCRPOP];			/*!< baryon fraction in cosmic rays */

  MyFloat CR_DeltaE[NUMCRPOP];		/*!< Specific Energy growth during timestep */
  MyFloat CR_DeltaN[NUMCRPOP];		/*!< baryon fraction growth during timestep */
#ifdef MACHNUM
  MyFloat CR_Gamma0[NUMCRPOP];
#endif

#ifdef CR_OUTPUT_INJECTION
  MyFloat CR_Specific_SupernovaHeatingRate;
#endif
#endif				/* COSMIC_RAYS */

#ifdef MACHNUM
  MyFloat Shock_MachNumber;	/*!< Mach number */
  MyFloat Shock_DecayTime;	/*!< Shock decay time */
#ifdef COSMIC_RAYS
  MyFloat Shock_DensityJump;	/*!< Density jump at the shock */
  MyFloat Shock_EnergyJump;	/*!< Energy jump at the shock */
  MyFloat PreShock_PhysicalDensity;	/*!< Specific energy in the preshock regime */
  MyFloat PreShock_PhysicalEnergy;	/*!< Density in the preshock regime */
  MyFloat PreShock_XCR;		/*!< XCR = PCR / Pth in the preshock regime */
#endif
#ifdef MACHSTATISTIC
  MyFloat Shock_DtEnergy;		/*!< Change of thermal specific energy at Shocks */
#endif
#ifdef OUTPUT_PRESHOCK_CSND
  MyFloat PreShock_PhysicalSoundSpeed;	/*!< Sound speed in the preshock regime */
  MyFloat PreShock_PhysicalDensity;	/*!< Specific energy in the preshock regime */
#endif
#endif				/* Mach number estimate */

#ifdef CHEMISTRY
  MyFloat elec;
  MyFloat HI;
  MyFloat HII;

  MyFloat HeI;
  MyFloat HeII;
  MyFloat HeIII;

  MyFloat H2I;
  MyFloat H2II;

  MyFloat HM;

  MyFloat Gamma;
  MyFloat t_elec, t_cool;
#endif

#ifdef RADTRANSFER
  MyFloat ET[6];                /* eddington tensor - symmetric -> only 6 elements needed */
  MyFloat Je;                   /* emissivity */
  MyFloat nHI;                  /* HI number density */
  MyFloat nHII;                 /* HII number density */
  MyFloat n_elec;               /* electron number density */
  MyFloat n_gamma;              /* photon number density */
#ifdef RADTRANSFER_FLUXLIMITER
  MyFloat Grad_ngamma[3];
#endif
#ifndef CG
  MyFloat n_gamma_old;
#endif
#endif

#if defined CS_MODEL
  MyFloat DensityOld;
#ifdef CS_FEEDBACK
  union
  {
    MyFloat       DensityAvg;		/*!< current baryonic mass density of particle */
    MyLongDouble dDensityAvg;
  } da;
  union
  {
    MyFloat       EntropyAvg;		/*!< current baryonic mass density of particle */
    MyLongDouble dEntropyAvg;
  } ea;
  MyFloat HotHsml;
  int     HotNgbNum;
  MyFloat DensPromotion;
  MyFloat TempPromotion;
#endif
#endif

#ifdef EOS_DEGENERATE
  MyFloat u;                            /* internal energy density */
  MyFloat temp;                         /* temperature */
  MyFloat xnuc[EOS_NSPECIES];           /* nuclear mass fractions */
#endif

#ifdef WAKEUP
  short int wakeup;             /*!< flag to wake up particle */
#endif
}
  *SphP,				/*!< holds SPH particle data on local processor */
  *DomainSphBuf;			/*!< buffer for SPH particle data in domain decomposition */

extern peanokey *DomainKeyBuf;

/* global state of system
*/
extern struct state_of_system
{
  double Mass,
    EnergyKin,
    EnergyPot,
    EnergyInt,
    EnergyTot,
    Momentum[4],
    AngMomentum[4],
    CenterOfMass[4],
    MassComp[6],
    EnergyKinComp[6],
    EnergyPotComp[6],
    EnergyIntComp[6], EnergyTotComp[6], MomentumComp[6][4], AngMomentumComp[6][4], CenterOfMassComp[6][4];
}
SysState, SysStateAtStart, SysStateAtEnd;

/* Various structures for communication during the gravity computation.
 */

extern struct data_index
{
  int Task;
  int Index;
  int IndexGet;
}
 *DataIndexTable;		/*!< the particles to be exported are grouped
                                  by task-number. This table allows the
                                  results to be disentangled again and to be
                                  assigned to the correct particle */

extern struct data_nodelist
{
  int NodeList[NODELISTLENGTH];
}
*DataNodeList;

extern struct gravdata_in
{
  MyFloat Pos[3];
#ifdef UNEQUALSOFTENINGS
  int Type;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
  MyFloat Soft;
#endif
#endif
  MyFloat OldAcc;
  int NodeList[NODELISTLENGTH];
}
 *GravDataIn,			/*!< holds particle data to be exported to other processors */
 *GravDataGet;			/*!< holds particle data imported from other processors */

extern struct gravdata_out
{
  MyLongDouble Acc[3];
#ifdef EVALPOTENTIAL
  MyLongDouble Potential;
#endif
#ifdef DISTORTIONTENSORPS
  MyLongDouble tidal_tensorps[3][3];
#endif
  int Ninteractions;
}
 *GravDataResult,		/*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
 *GravDataOut;			/*!< holds partial results received from other processors. This will overwrite the GravDataIn array */

extern struct potdata_out
{
  MyLongDouble Potential;
}
 *PotDataResult,		/*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
 *PotDataOut;			/*!< holds partial results received from other processors. This will overwrite the GravDataIn array */

/*! Header for the standard file format.
 */
extern struct io_header
{
  int npart[6];			/*!< number of particles of each type in this file */
  double mass[6];		/*!< mass of particles of each type. If 0, then the masses are explicitly
				   stored in the mass-block of the snapshot file, otherwise they are omitted */
#ifdef COSMIC_RAYS
  double SpectralIndex_CR_Pop[NUMCRPOP]; /*!< spectral indices of cosmic ray populations */
#endif
  double time;			/*!< time of snapshot file */
  double redshift;		/*!< redshift of snapshot file */
  int flag_sfr;			/*!< flags whether the simulation was including star formation */
  int flag_feedback;		/*!< flags whether feedback was included (obsolete) */
  unsigned int npartTotal[6];	/*!< total number of particles of each type in this snapshot. This can be
				   different from npart if one is dealing with a multi-file snapshot. */
  int flag_cooling;		/*!< flags whether cooling was included  */
  int num_files;		/*!< number of files in multi-file snapshot */
  double BoxSize;		/*!< box-size of simulation in case periodic boundaries were used */
  double Omega0;		/*!< matter density in units of critical density */
  double OmegaLambda;		/*!< cosmological constant parameter */
  double HubbleParam;		/*!< Hubble parameter in units of 100 km/sec/Mpc */
  int flag_stellarage;		/*!< flags whether the file contains formation times of star particles */
  int flag_metals;		/*!< flags whether the file contains metallicity values for gas and star
				   particles */
  unsigned int npartTotalHighWord[6];	/*!< High word of the total number of particles of each type */
  int flag_entropy_instead_u;	/*!< flags that IC-file contains entropy instead of u */
  int flag_doubleprecision;	/*!< flags that snapshot contains double-precision instead of single precision */

  int flag_ic_info;             /*!< flag to inform whether IC files are generated with ordinary Zeldovich approximation,
                                     or whether they ocontains 2nd order lagrangian perturbation theory initial conditions.
                                     For snapshots files, the value informs whether the simulation was evolved from
                                     Zeldoch or 2lpt ICs. Encoding is as follows:
                                        FLAG_ZELDOVICH_ICS     (1)   - IC file based on Zeldovich
                                        FLAG_SECOND_ORDER_ICS  (2)   - Special IC-file containing 2lpt masses
                                        FLAG_EVOLVED_ZELDOVICH (3)   - snapshot evolved from Zeldovich ICs
                                        FLAG_EVOLVED_2LPT      (4)   - snapshot evolved from 2lpt ICs
                                        FLAG_NORMALICS_2LPT    (5)   - standard gadget file format with 2lpt ICs
                                     All other values, including 0 are interpreted as "don't know" for backwards compatability.
                                 */
  float lpt_scalingfactor;      /*!< scaling factor for 2lpt initial conditions */

#ifdef COSMIC_RAYS
  char fill[48-8*NUMCRPOP];	/*!< fills to 256 Bytes */
#else
  char fill[48];		/*!< fills to 256 Bytes */
#endif

}
header;				/*!< holds header for snapshot files */

enum iofields
{ IO_POS,
  IO_VEL,
  IO_ID,
  IO_MASS,
  IO_SECONDORDERMASS,
  IO_U,
  IO_RHO,
  IO_NE,
  IO_NH,
  IO_HSML,
  IO_SFR,
  IO_AGE,
  IO_Z,
  IO_BHMASS,
  IO_BHMDOT,
  IO_BHPROGS,
  IO_BHMBUB,
  IO_BHMINI,
  IO_BHMRAD,
  IO_POT,
  IO_ACCEL,
  IO_CR_C0,
  IO_CR_Q0,
  IO_CR_P0,
  IO_CR_E0,
  IO_CR_n0,
  IO_CR_ThermalizationTime,
  IO_CR_DissipationTime,
  IO_ELECT,
  IO_HI,
  IO_HII,
  IO_HeI,
  IO_HeII,
  IO_HeIII,
  IO_H2I,
  IO_H2II,
  IO_HM,
  IO_DTENTR,
  IO_STRESSDIAG,
  IO_STRESSOFFDIAG,
  IO_STRESSBULK,
  IO_SHEARCOEFF,
  IO_TSTP,
  IO_BFLD,
  IO_BSMTH,
  IO_DBDT,
  IO_DIVB,
  IO_ABVC,
  IO_AMDC,
  IO_PHI,
  IO_ROTB,
  IO_SROTB,
  IO_COOLRATE,
  IO_CONDRATE,
  IO_DENN,
  IO_EGYPROM,
  IO_EGYCOLD,
  IO_MACH,
  IO_DTENERGY,
  IO_PRESHOCK_CSND,
  IO_PRESHOCK_DENSITY,
  IO_PRESHOCK_ENERGY,
  IO_PRESHOCK_XCR,
  IO_DENSITY_JUMP,
  IO_ENERGY_JUMP,
  IO_CRINJECT,
  IO_TIDALTENSORPS,
  IO_DISTORTIONTENSORPS,
  IO_EULERA,
  IO_EULERB,
  IO_FLOW_DETERMINANT,
  IO_PHASE_SPACE_DETERMINANT,
  IO_ANNIHILATION_RADIATION,
  IO_STREAM_DENSITY,
  IO_EOSTEMP,
  IO_EOSXNUC,
  IO_PRESSURE,
  IO_nHII,
  IO_RADGAMMA,
  IO_LAST_CAUSTIC,
  IO_SHEET_ORIENTATION,
  IO_INIT_DENSITY,
  IO_CAUSTIC_COUNTER,
  IO_SHELL_INFO,
  IO_DMHSML,                    /* for 'SUBFIND_RESHUFFLE_CATALOGUE' option */
  IO_DMDENSITY,
  IO_DMVELDISP,
  IO_DMHSML_V,                 /* for 'SUBFIND_RESHUFFLE_CATALOGUE_WITH_VORONOI' option */
  IO_DMDENSITY_V,

  IO_LASTENTRY			/* This should be kept - it signals the end of the list */
};

/*
 * Variables for Tree
 * ------------------
 */

extern struct NODE
{
  MyFloat len;			/*!< sidelength of treenode */
  MyFloat center[3];		/*!< geometrical center of node */

#ifdef RADTRANSFER
  MyFloat stellar_mass;         /*!< mass in stars in the node*/
  MyFloat stellar_s[3];         /*!< enter of mass for the stars in the node*/
#endif

#ifdef ADAPTIVE_GRAVSOFT_FORGAS
  MyFloat maxsoft;		/*!< hold the maximum gravitational softening of particles in the
				   node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
  union
  {
    int suns[8];		/*!< temporary pointers to daughter nodes */
    struct
    {
      MyFloat s[3];		/*!< center of mass of node */
      MyFloat mass;		/*!< mass of node */
      unsigned int bitflags;	/*!< flags certain node properties */
      int sibling;		/*!< this gives the next node in the walk in case the current node can be used */
      int nextnode;		/*!< this gives the next node in case the current node needs to be opened */
      int father;		/*!< this gives the parent node of each node (or -1 if we have the root node) */
    }
    d;
  }
  u;
#ifdef SCALARFIELD
  MyFloat s_dm[3];
  MyFloat mass_dm;
#endif
  int Ti_current;
}
 *Nodes_base,			/*!< points to the actual memory allocted for the nodes */
 *Nodes;			/*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
				   gives the first allocated node */

extern struct extNODE
{
  MyLongDouble dp[3];
#ifdef SCALARFIELD
  MyLongDouble dp_dm[3];
  MyFloat vs_dm[3];
#endif
#ifdef FLTROUNDOFFREDUCTION
  MyFloat s_base[3];
  MyFloat len_base;
#ifdef SCALARFIELD
  MyFloat s_dm_base[3];
#endif
#endif
  MyFloat vs[3];
  MyFloat vmax;
  MyFloat divVmax;
  MyFloat hmax;			/*!< maximum SPH smoothing length in node. Only used for gas particles */
  int Ti_lastkicked;
  int Flag;
}
 *Extnodes, *Extnodes_base;

extern int MaxNodes;		/*!< maximum allowed number of internal nodes */
extern int Numnodestree;	/*!< number of (internal) nodes in each tree */

extern int *Nextnode;		/*!< gives next node in tree walk  (nodes array) */
extern int *Father;		/*!< gives parent node in tree (Prenodes array) */

#ifdef STATICNFW
extern double Rs, R200;
extern double Dc;
extern double RhoCrit, V200;
extern double fac;
#endif

#ifdef CHEMISTRY
/* ----- chemistry part ------- */

#define H_number_fraction 0.76
#define He_number_fraction 0.06

/* ----- Tables ------- */
extern double T[N_T], J0_nu[N_nu], J_nu[N_nu], nu[N_nu];
extern double k1a[N_T], k2a[N_T], k3a[N_T], k4a[N_T], k5a[N_T], k6a[N_T], k7a[N_T], k8a[N_T], k9a[N_T],
  k10a[N_T], k11a[N_T];
extern double k12a[N_T], k13a[N_T], k14a[N_T], k15a[N_T], k16a[N_T], k17a[N_T], k18a[N_T], k19a[N_T],
  k20a[N_T], k21a[N_T];
extern double ciHIa[N_T], ciHeIa[N_T], ciHeIIa[N_T], ciHeISa[N_T], reHIIa[N_T], brema[N_T];
extern double ceHIa[N_T], ceHeIa[N_T], ceHeIIa[N_T], reHeII1a[N_T], reHeII2a[N_T], reHeIIIa[N_T];

/* cross-sections */
#ifdef RADIATION
extern double sigma24[N_nu], sigma25[N_nu], sigma26[N_nu], sigma27[N_nu], sigma28[N_nu], sigma29[N_nu],
  sigma30[N_nu], sigma31[N_nu];
#endif
#endif

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

/* This file is part of phonopy. */

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

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

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

/* * Neither the name of the phonopy project nor the names of its */
/*   contributors may be used to endorse or promote products derived */
/*   from this software without specific prior written permission. */

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

#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
#include <numpy/arrayobject.h>
#include <lapacke.h>
#include "phonoc_array.h"
#include "phonon4_h/fc4.h"
#include "phonon4_h/real_to_reciprocal.h"
#include "phonon4_h/frequency_shift.h"

static PyObject * py_get_fc4_normal_for_frequency_shift(PyObject *self, PyObject *args);
static PyObject * py_get_fc4_frequency_shifts(PyObject *self, PyObject *args);
static PyObject * py_real_to_reciprocal4(PyObject *self, PyObject *args);
static PyObject * py_reciprocal_to_normal4(PyObject *self, PyObject *args);
static PyObject * py_set_phonons_grid_points(PyObject *self, PyObject *args);
static PyObject * py_distribute_fc4(PyObject *self, PyObject *args);
static PyObject * py_rotate_delta_fc3s_elem(PyObject *self, PyObject *args);
static PyObject * py_set_translational_invariance_fc4(PyObject *self,
						      PyObject *args);
static PyObject * py_set_permutation_symmetry_fc4(PyObject *self,
						  PyObject *args);
static PyObject * py_get_drift_fc4(PyObject *self, PyObject *args);

static PyMethodDef functions[] = {
  {"fc4_normal_for_frequency_shift", py_get_fc4_normal_for_frequency_shift, METH_VARARGS, "Calculate fc4 normal for frequency shift"},
  {"fc4_frequency_shifts", py_get_fc4_frequency_shifts, METH_VARARGS, "Calculate fc4 frequency shift"},
  {"real_to_reciprocal4", py_real_to_reciprocal4, METH_VARARGS, "Transform fc4 of real space to reciprocal space"},
  {"reciprocal_to_normal4", py_reciprocal_to_normal4, METH_VARARGS, "Transform fc4 of reciprocal space to normal coordinate in special case for frequency shift"},
  {"phonons_grid_points", py_set_phonons_grid_points, METH_VARARGS, "Set phonons on grid points"},
  {"distribute_fc4", py_distribute_fc4, METH_VARARGS, "Distribute least fc4 to full fc4"},
  {"rotate_delta_fc3s_elem", py_rotate_delta_fc3s_elem, METH_VARARGS, "Rotate delta fc3s for a set of atomic indices"},
  {"translational_invariance_fc4", py_set_translational_invariance_fc4, METH_VARARGS, "Set translational invariance for fc4"},
  {"permutation_symmetry_fc4", py_set_permutation_symmetry_fc4, METH_VARARGS, "Set permutation symmetry for fc4"},
  {"drift_fc4", py_get_drift_fc4, METH_VARARGS, "Get drifts of fc4"},
  {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC init_phono4py(void)
{
  Py_InitModule3("_phono4py", functions, "C-extension for phono4py\n\n...\n");
  return;
}

static PyObject * py_get_fc4_normal_for_frequency_shift(PyObject *self,
							PyObject *args)
{
  PyArrayObject* fc4_normal_py;
  PyArrayObject* frequencies_py;
  PyArrayObject* eigenvectors_py;
  PyArrayObject* grid_points1_py;
  PyArrayObject* grid_address_py;
  PyArrayObject* mesh_py;
  PyArrayObject* fc4_py;
  PyArrayObject* shortest_vectors_py;
  PyArrayObject* multiplicity_py;
  PyArrayObject* masses_py;
  PyArrayObject* p2s_map_py;
  PyArrayObject* s2p_map_py;
  PyArrayObject* band_indicies_py;
  double cutoff_frequency;
  int grid_point0;

  if (!PyArg_ParseTuple(args, "OOOiOOOOOOOOOOd",
			&fc4_normal_py,
			&frequencies_py,
			&eigenvectors_py,
			&grid_point0,
			&grid_points1_py,
			&grid_address_py,
			&mesh_py,
			&fc4_py,
			&shortest_vectors_py,
			&multiplicity_py,
			&masses_py,
			&p2s_map_py,
			&s2p_map_py,
			&band_indicies_py,
			&cutoff_frequency)) {
    return NULL;
  }

  double* fc4_normal = (double*)fc4_normal_py->data;
  double* freqs = (double*)frequencies_py->data;
  /* npy_cdouble and lapack_complex_double may not be compatible. */
  /* So eigenvectors should not be used in Python side */
  lapack_complex_double* eigvecs =
    (lapack_complex_double*)eigenvectors_py->data;
  Iarray* grid_points1 = convert_to_iarray(grid_points1_py);
  const int* grid_address = (int*)grid_address_py->data;
  const int* mesh = (int*)mesh_py->data;
  double* fc4 = (double*)fc4_py->data;
  Darray* svecs = convert_to_darray(shortest_vectors_py);
  Iarray* multi = convert_to_iarray(multiplicity_py);
  const double* masses = (double*)masses_py->data;
  const int* p2s = (int*)p2s_map_py->data;
  const int* s2p = (int*)s2p_map_py->data;
  Iarray* band_indicies = convert_to_iarray(band_indicies_py);

  get_fc4_normal_for_frequency_shift(fc4_normal,
				     freqs,
				     eigvecs,
				     grid_point0,
				     grid_points1,
				     grid_address,
				     mesh,
				     fc4,
				     svecs,
				     multi,
				     masses,
				     p2s,
				     s2p,
				     band_indicies,
				     cutoff_frequency);

  free(grid_points1);
  free(svecs);
  free(multi);
  free(band_indicies);

  Py_RETURN_NONE;
}

static PyObject * py_get_fc4_frequency_shifts(PyObject *self, PyObject *args)
{
  PyArrayObject* frequency_shifts_py;
  PyArrayObject* fc4_normal_py;
  PyArrayObject* frequencies_py;
  PyArrayObject* grid_points1_py;
  PyArrayObject* temperatures_py;
  PyArrayObject* band_indicies_py;
  double unit_conversion_factor;

  if (!PyArg_ParseTuple(args, "OOOOOOd",
			&frequency_shifts_py,
			&fc4_normal_py,
			&frequencies_py,
			&grid_points1_py,
			&temperatures_py,
			&band_indicies_py,
			&unit_conversion_factor)) {
    return NULL;
  }

  double* freq_shifts = (double*)frequency_shifts_py->data;
  double* fc4_normal = (double*)fc4_normal_py->data;
  double* freqs = (double*)frequencies_py->data;
  Iarray* grid_points1 = convert_to_iarray(grid_points1_py);
  Darray* temperatures = convert_to_darray(temperatures_py);
  int* band_indicies = (int*)band_indicies_py->data;
  const int num_band0 = (int)band_indicies_py->dimensions[0];
  const int num_band = (int)frequencies_py->dimensions[1];

  get_fc4_frequency_shifts(freq_shifts,
			   fc4_normal,
			   freqs,
			   grid_points1,
			   temperatures,
			   band_indicies,
			   num_band0,
			   num_band,
			   unit_conversion_factor);

  free(grid_points1);
  free(temperatures);

  Py_RETURN_NONE;
}

static PyObject * py_real_to_reciprocal4(PyObject *self, PyObject *args)
{
  PyArrayObject* fc4_py;
  PyArrayObject* fc4_reciprocal_py;
  PyArrayObject* q_py;
  PyArrayObject* shortest_vectors;
  PyArrayObject* multiplicity;
  PyArrayObject* p2s_map;
  PyArrayObject* s2p_map;

  if (!PyArg_ParseTuple(args, "OOOOOOO",
			&fc4_reciprocal_py,
			&fc4_py,
			&q_py,
			&shortest_vectors,
			&multiplicity,
			&p2s_map,
			&s2p_map)) {
    return NULL;
  }

  double* fc4 = (double*)fc4_py->data;
  lapack_complex_double* fc4_reciprocal =
    (lapack_complex_double*)fc4_reciprocal_py->data;
  Darray* svecs = convert_to_darray(shortest_vectors);
  Iarray* multi = convert_to_iarray(multiplicity);
  const int* p2s = (int*)p2s_map->data;
  const int* s2p = (int*)s2p_map->data;
  const double* q = (double*)q_py->data;

  real_to_reciprocal4(fc4_reciprocal,
		      q,
		      fc4,
		      svecs,
		      multi,
		      p2s,
		      s2p);

  free(svecs);
  free(multi);

  Py_RETURN_NONE;
}

static PyObject * py_reciprocal_to_normal4(PyObject *self, PyObject *args)
{
  PyArrayObject* fc4_normal_py;
  PyArrayObject* fc4_reciprocal_py;
  PyArrayObject* frequencies_py;
  PyArrayObject* eigenvectors_py;
  PyArrayObject* grid_points_py;
  PyArrayObject* masses_py;
  PyArrayObject* band_indicies_py;
  double cutoff_frequency;

  if (!PyArg_ParseTuple(args, "OOOOOOOd",
			&fc4_normal_py,
			&fc4_reciprocal_py,
			&frequencies_py,
			&eigenvectors_py,
			&grid_points_py,
			&masses_py,
			&band_indicies_py,
			&cutoff_frequency)) {
    return NULL;
  }

  lapack_complex_double* fc4_normal =
    (lapack_complex_double*)fc4_normal_py->data;
  const lapack_complex_double* fc4_reciprocal =
    (lapack_complex_double*)fc4_reciprocal_py->data;
  const lapack_complex_double* eigenvectors =
    (lapack_complex_double*)eigenvectors_py->data;
  const double* frequencies = (double*)frequencies_py->data;
  const int* grid_points = (int*)grid_points_py->data;
  const double* masses = (double*)masses_py->data;
  const int* band_indices = (int*)band_indicies_py->data;
  const int num_band0 = (int)band_indicies_py->dimensions[0];
  const int num_band = (int)frequencies_py->dimensions[1];

  reciprocal_to_normal4(fc4_normal,
			fc4_reciprocal,
			frequencies + grid_points[0] * num_band,
			frequencies + grid_points[1] * num_band,
			eigenvectors + grid_points[0] * num_band * num_band,
			eigenvectors + grid_points[1] * num_band * num_band,
			masses,
			band_indices,
			num_band0,
			num_band,
			cutoff_frequency);

  Py_RETURN_NONE;
}

static PyObject * py_set_phonons_grid_points(PyObject *self, PyObject *args)
{
  PyArrayObject* frequencies;
  PyArrayObject* eigenvectors;
  PyArrayObject* phonon_done_py;
  PyArrayObject* grid_points_py;
  PyArrayObject* grid_address_py;
  PyArrayObject* mesh_py;
  PyArrayObject* shortest_vectors_fc2;
  PyArrayObject* multiplicity_fc2;
  PyArrayObject* fc2_py;
  PyArrayObject* atomic_masses_fc2;
  PyArrayObject* p2s_map_fc2;
  PyArrayObject* s2p_map_fc2;
  PyArrayObject* reciprocal_lattice;
  PyArrayObject* born_effective_charge;
  PyArrayObject* q_direction;
  PyArrayObject* dielectric_constant;
  double nac_factor, unit_conversion_factor;
  char uplo;

  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOdOOOOdc",
			&frequencies,
			&eigenvectors,
			&phonon_done_py,
			&grid_points_py,
			&grid_address_py,
			&mesh_py,
			&fc2_py,
			&shortest_vectors_fc2,
			&multiplicity_fc2,
			&atomic_masses_fc2,
			&p2s_map_fc2,
			&s2p_map_fc2,
			&unit_conversion_factor,
			&born_effective_charge,
			&dielectric_constant,
			&reciprocal_lattice,
			&q_direction,
			&nac_factor,
			&uplo)) {
    return NULL;
  }

  double* born;
  double* dielectric;
  double *q_dir;
  Darray* freqs = convert_to_darray(frequencies);
  /* npy_cdouble and lapack_complex_double may not be compatible. */
  /* So eigenvectors should not be used in Python side */
  Carray* eigvecs = convert_to_carray(eigenvectors);
  char* phonon_done = (char*)phonon_done_py->data;
  Iarray* grid_points = convert_to_iarray(grid_points_py);
  const int* grid_address = (int*)grid_address_py->data;
  const int* mesh = (int*)mesh_py->data;
  Darray* fc2 = convert_to_darray(fc2_py);
  Darray* svecs_fc2 = convert_to_darray(shortest_vectors_fc2);
  Iarray* multi_fc2 = convert_to_iarray(multiplicity_fc2);
  const double* masses_fc2 = (double*)atomic_masses_fc2->data;
  const int* p2s_fc2 = (int*)p2s_map_fc2->data;
  const int* s2p_fc2 = (int*)s2p_map_fc2->data;
  const double* rec_lat = (double*)reciprocal_lattice->data;
  if ((PyObject*)born_effective_charge == Py_None) {
    born = NULL;
  } else {
    born = (double*)born_effective_charge->data;
  }
  if ((PyObject*)dielectric_constant == Py_None) {
    dielectric = NULL;
  } else {
    dielectric = (double*)dielectric_constant->data;
  }
  if ((PyObject*)q_direction == Py_None) {
    q_dir = NULL;
  } else {
    q_dir = (double*)q_direction->data;
  }

  set_phonons_for_frequency_shift(freqs,
				  eigvecs,
				  phonon_done,
				  grid_points,
				  grid_address,
				  mesh,
				  fc2,
				  svecs_fc2,
				  multi_fc2,
				  masses_fc2,
				  p2s_fc2,
				  s2p_fc2,
				  unit_conversion_factor,
				  born,
				  dielectric,
				  rec_lat,
				  q_dir,
				  nac_factor,
				  uplo);

  free(freqs);
  free(eigvecs);
  free(grid_points);
  free(fc2);
  free(svecs_fc2);
  free(multi_fc2);

  Py_RETURN_NONE;
}

static PyObject * py_distribute_fc4(PyObject *self, PyObject *args)
{
  PyArrayObject* fc4_copy_py;
  PyArrayObject* fc4_py;
  int fourth_atom;
  PyArrayObject* rotation_cart_inv;
  PyArrayObject* atom_mapping_py;

  if (!PyArg_ParseTuple(args, "OOiOO",
			&fc4_copy_py,
			&fc4_py,
			&fourth_atom,
			&atom_mapping_py,
			&rotation_cart_inv)) {
    return NULL;
  }

  double* fc4_copy = (double*)fc4_copy_py->data;
  const double* fc4 = (double*)fc4_py->data;
  const double* rot_cart_inv = (double*)rotation_cart_inv->data;
  const int* atom_mapping = (int*)atom_mapping_py->data;
  const int num_atom = (int)atom_mapping_py->dimensions[0];

  return PyInt_FromLong((long) distribute_fc4(fc4_copy,
					      fc4,
					      fourth_atom,
					      atom_mapping,
					      num_atom,
					      rot_cart_inv));
}

static PyObject * py_rotate_delta_fc3s_elem(PyObject *self, PyObject *args)
{
  PyArrayObject* rotated_delta_fc3s_py;
  PyArrayObject* delta_fc3s_py;
  PyArrayObject* atom_mappings_of_rotations_py;
  PyArrayObject* site_symmetries_cartesian_py;
  int atom1, atom2, atom3;

  if (!PyArg_ParseTuple(args, "OOOOiii",
			&rotated_delta_fc3s_py,
			&delta_fc3s_py,
			&atom_mappings_of_rotations_py,
			&site_symmetries_cartesian_py,
			&atom1,
			&atom2,
			&atom3)) {
    return NULL;
  }

  double* rotated_delta_fc3s = (double*)rotated_delta_fc3s_py->data;
  const double* delta_fc3s = (double*)delta_fc3s_py->data;
  const int* rot_map_syms = (int*)atom_mappings_of_rotations_py->data;
  const double* site_syms_cart = (double*)site_symmetries_cartesian_py->data;
  const int num_rot = (int)site_symmetries_cartesian_py->dimensions[0];
  const int num_delta_fc3s = (int)delta_fc3s_py->dimensions[0];
  const int num_atom = (int)delta_fc3s_py->dimensions[1];

  return PyInt_FromLong((long) rotate_delta_fc3s_elem(rotated_delta_fc3s,
						      delta_fc3s,
						      rot_map_syms,
						      site_syms_cart,
						      num_rot,
						      num_delta_fc3s,
						      atom1,
						      atom2,
						      atom3,
						      num_atom));
}

static PyObject * py_set_translational_invariance_fc4(PyObject *self,
						      PyObject *args)
{
  PyArrayObject* fc4_py;
  int index;

  if (!PyArg_ParseTuple(args, "Oi",
			&fc4_py,
			&index)) {
    return NULL;
  }

  double* fc4 = (double*)fc4_py->data;
  const int num_atom = (int)fc4_py->dimensions[0];

  set_translational_invariance_fc4_per_index(fc4, num_atom, index);

  Py_RETURN_NONE;
}

static PyObject * py_set_permutation_symmetry_fc4(PyObject *self, PyObject *args)
{
  PyArrayObject* fc4_py;

  if (!PyArg_ParseTuple(args, "O",
			&fc4_py)) {
    return NULL;
  }

  double* fc4 = (double*)fc4_py->data;
  const int num_atom = (int)fc4_py->dimensions[0];

  set_permutation_symmetry_fc4(fc4, num_atom);

  Py_RETURN_NONE;
}

static PyObject * py_get_drift_fc4(PyObject *self, PyObject *args)
{
  PyArrayObject* fc4_py;

  if (!PyArg_ParseTuple(args, "O",
			&fc4_py)) {
    return NULL;
  }

  double* fc4 = (double*)fc4_py->data;
  const int num_atom = (int)fc4_py->dimensions[0];

  int i;
  double drift[4];
  PyObject* drift_py;

  get_drift_fc4(drift, fc4, num_atom);
  drift_py = PyList_New(4);

  for (i = 0; i < 4; i++) {
    PyList_SetItem(drift_py, i, PyFloat_FromDouble(drift[i]));
  }

  return drift_py;
}

ext: c

max_stars_repo_name: atztogo/forcefit

max_issues_repo_name: atztogo/forcefit

max_forks_repo_name: atztogo/forcefit


text:
/**
 *
 * @file sgecon.c
 *
 *  PLASMA computational routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Jakub Kurzak
 * @author Ichitaro Yamazaki
 * @author Mathieu Faverge
 * @date 2012-10-04
 * @generated s Tue Jan  7 11:45:10 2014
 *
 **/
#include <lapacke.h>
#include "common.h"
#include <math.h>
#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 * @ingroup float
 *
 *  PLASMA_sgecon - estimates the reciprocal of the condition number
 *  of a general complex matrix A, in either the 1-norm or the infinity-norm,
 *  using the LU factorization computed by PLASMA_sgetrf().
 *
 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the condition
 *  number is computed as
 *
 *    \f[ rcond = \frac{1}{\|\|A\-\| \times \-\|A^{-1}\|\|} \f]
 *
 *******************************************************************************
 *
 * @param[in] norm
 *          Specifies whether the 1-norm condition number
 *          or the infinity-norm condition number is required:
 *          = PlasmaOneNorm: One norm
 *          = PlasmaInfNorm: Infinity norm
 *
 * @param[in] N
 *          The order of the matrix A. N >= 0.
 *
 * @param[in] A
 *          The N-by-N matrix A.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[in] Anorm
 *          If norm = PlasmaOneNorm, the 1-norm of the original matrix A.
 *          If norm = PlasmaInfNorm, the infinity-norm of the original matrix A.
 *
 * \param[out] rcond
 *          The reciprocal of the condition number of the matrix A,
 *          computed as stated above.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgecon_Tile
 * @sa PLASMA_sgecon_Tile_Async
 * @sa PLASMA_cgecon
 * @sa PLASMA_dgecon
 * @sa PLASMA_sgecon
 *
 ******************************************************************************/
int PLASMA_sgecon(PLASMA_enum norm, int N,
                  float *A, int LDA, float Anorm, float *rcond)
{
    int NB;
    int status;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    PLASMA_desc descA;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgecon", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    /* Check input arguments */
    if (norm != PlasmaOneNorm && norm != PlasmaInfNorm) {
        plasma_error("PLASMA_sgecon", "illegal value of norm");
        return -1;
    }
    if (N < 0) {
        plasma_error("PLASMA_sgecon", "illegal value of N");
        return -2;
    }
    if (LDA < max(1,N)) {
        plasma_error("PLASMA_sgecon", "illegal value of LDA");
        return -4;
    }
    if (Anorm < 0.) {
        plasma_error("PLASMA_sgecon", "illegal value of Anorm");
        return -5;
    }

    /* Quick return */
    *rcond = (float)0.;
    if (N == 0) {
        *rcond = (float)1.;
        return PLASMA_SUCCESS;
    }
    else if (Anorm == 0.) {
        return PLASMA_SUCCESS;
    }

    /* Tune NB depending on M, N & NRHS; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_SGESV, N, N, 0);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgecon", "plasma_tune() failed");
        return status;
    }
    /* Set NT */
    NB = PLASMA_NB;
    plasma_sequence_create(plasma, &sequence);

    if (PLASMA_TRANSLATION == PLASMA_OUTOFPLACE) {
        plasma_sooplap2tile(
            descA, A, NB, NB, LDA, N, 0, 0, N, N,
            sequence, &request, plasma_desc_mat_free(&(descA)));
    } else {
        plasma_siplap2tile(
            descA, A, NB, NB, LDA, N, 0, 0, N, N, sequence, &request);
    }

    /* Call the tile interface */
    PLASMA_sgecon_Tile_Async(norm, &descA, Anorm, rcond, sequence, &request);

    if (PLASMA_TRANSLATION == PLASMA_OUTOFPLACE) {
        plasma_dynamic_sync();
        plasma_desc_mat_free(&descA);
    } else {
        plasma_siptile2lap(descA, A, NB, NB, LDA, N, sequence, &request);
        plasma_dynamic_sync();
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup float_Tile
 *
 *  PLASMA_sgecon_Tile - estimates the reciprocal of the condition number
 *  of a general complex matrix A, in either the 1-norm or the infinity-norm.
 *  Tile equivalent of PLASMA_sgecon().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] norm
 *          Specifies whether the 1-norm condition number
 *          or the infinity-norm condition number is required:
 *          = PlasmaOneNorm: One norm
 *          = PlasmaInfNorm: Infinity norm
 *
 * @param[in] A
 *          The N-by-N matrix A.
 *
 * @param[in] Anorm
 *          If norm = PlasmaOneNorm, the 1-norm of the original matrix A.
 *          If norm = PlasmaInfNorm, the infinity-norm of the original matrix A.
 *
 * \param[out] rcond
 *          The reciprocal of the condition number of the matrix A,
 *          computed as stated above.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgecon
 * @sa PLASMA_sgecon_Tile_Async
 * @sa PLASMA_cgecon_Tile
 * @sa PLASMA_dgecon_Tile
 * @sa PLASMA_sgecon_Tile
 *
 ******************************************************************************/
int PLASMA_sgecon_Tile(PLASMA_enum norm, PLASMA_desc *A, float Anorm, float *rcond)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgecon_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_sgecon_Tile_Async(norm, A, Anorm, rcond, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup float_Tile_Async
 *
 *  PLASMA_sgecon_Tile_Async - estimates the reciprocal of the condition number
 *  of a general complex matrix A, in either the 1-norm or the infinity-norm.
 *  Non-blocking equivalent of PLASMA_sgecon_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_sgecon
 * @sa PLASMA_sgecon_Tile
 * @sa PLASMA_cgecon_Tile_Async
 * @sa PLASMA_dgecon_Tile_Async
 * @sa PLASMA_sgecon_Tile_Async
 *
 ******************************************************************************/
int PLASMA_sgecon_Tile_Async(PLASMA_enum norm, PLASMA_desc *A, float Anorm, float *rcond,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    plasma_context_t *plasma;
    PLASMA_desc descA;
    PLASMA_desc descW;
    float *workN;
    float Ainvnorm;
    int kase, kase1;
    int isave[3], itrs = 0;
    int fallback = PLASMA_FALSE;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_sgecon_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_sgecon_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_sgecon_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if ( A->m != A->n ) {
        plasma_error("PLASMA_sgecon_Tile_Async", "invalid A descriptor (not square)");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_sgecon_Tile_Async", "invalid first descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descA = *A;
    }
    /* Check input arguments */
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_sgecon_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }

    /* Quick return */
    *rcond = (float)0.;
    if (descA.m == 0) {
        *rcond = (float)1.;
        return PLASMA_SUCCESS;
    }
    else if (Anorm == 0.) {
        return PLASMA_SUCCESS;
    }

    /* Estimate the norm of inv(A). */
    Ainvnorm = (float)0.;
    if (norm == PlasmaOneNorm)
        kase1 = 1;
    else
        kase1 = 2;
    kase = 0;

#if defined(REAL)
    int *isgn = (int*)plasma_shared_alloc(plasma, descA.m, PlasmaInteger);
#endif

    workN = (float*)plasma_shared_alloc(plasma, descA.m, PlasmaRealFloat);
    plasma_sdesc_alloc( descW, descA.mb, descA.nb,
                        descA.m, 1, 0, 0, descA.m, 1, plasma_desc_mat_free(&(descW)));

    do {
        itrs ++;
#if defined(REAL)
        LAPACKE_slacn2_work( descA.m, workN, descW.mat, isgn, &Ainvnorm, &kase, isave);
#else
        LAPACKE_slacn2_work( descA.m, workN, descW.mat,       &Ainvnorm, &kase, isave);
#endif

#define FALLBACK
#ifdef  FALLBACK
        /*
         * Fall back to LAPACK
         */
        if( isnan(Ainvnorm) || isinf(Ainvnorm) || Ainvnorm > LAPACKE_slamch('O') ) {
            int info;
            float *Atmp = (float*)malloc(descA.m * descA.n * sizeof(float));

            plasma_sooptile2lap( descA, Atmp, descA.mb, descA.nb, descA.m, descA.n, sequence, request);
            plasma_dynamic_sync();
            info = LAPACKE_sgecon(LAPACK_COL_MAJOR, lapack_const(norm), descA.n, Atmp, descA.m, Anorm, rcond);
            free(Atmp);
            fallback = PLASMA_TRUE;
            sequence->status = info;
            kase = 0;
        }
#endif

        if (kase != 0) {
            if (kase == kase1) {
                /* Multiply by inv(L). */
                plasma_parallel_call_9(plasma_pstrsm,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaLower,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_enum, PlasmaUnit,
                    float, 1.0,
                    PLASMA_desc, descA,
                    PLASMA_desc, descW,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);

                /* Multiply by inv(U). */
                plasma_parallel_call_9(plasma_pstrsm,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaUpper,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_enum, PlasmaNonUnit,
                    float, 1.0,
                    PLASMA_desc, descA,
                    PLASMA_desc, descW,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            else {
                /* Multiply by inv(U**T). */
                plasma_parallel_call_9(plasma_pstrsm,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaUpper,
                    PLASMA_enum, PlasmaTrans,
                    PLASMA_enum, PlasmaNonUnit,
                    float, 1.0,
                    PLASMA_desc, descA,
                    PLASMA_desc, descW,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);

                /* Multiply by inv(L**T). */
                plasma_parallel_call_9(plasma_pstrsm,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaLower,
                    PLASMA_enum, PlasmaTrans,
                    PLASMA_enum, PlasmaUnit,
                    float, 1.0,
                    PLASMA_desc, descA,
                    PLASMA_desc, descW,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
        }
        plasma_dynamic_sync();
    }
    while (kase != 0);

    /* Compute the estimate of the reciprocal condition number. */
    if ((Ainvnorm != 0.0) && (fallback == PLASMA_FALSE)) {
        *rcond = ((float)1.0 / Ainvnorm) / Anorm;
    }

#if defined(REAL)
    plasma_shared_free(plasma, isgn);
#endif
    plasma_shared_free(plasma, workN);
    plasma_desc_mat_free(&descW);

    return PLASMA_SUCCESS;
}


text:
#include <stdio.h>
#include <gsl/gsl_matrix.h>

#define MATRIX_SIZE 3

static const double matrix_1[MATRIX_SIZE][MATRIX_SIZE] = {
  {1., 4., 2.},
  {-1., -2., 1.},
  {3., 20., 19.},
};
static const double vector_1[] = {8., 3., 71.};

static const double matrix_2[MATRIX_SIZE][MATRIX_SIZE] = {
  {1.,    0.1,    0.225},
  {-0.15, 1.,     0.2},
  {-0.4,  -0.025, 1.}
};
static double const vector_2[] = {0.1, 0.9, 0.8};

void ulDecomposition(const gsl_matrix *A, const gsl_vector *b) {
  gsl_matrix *L = gsl_matrix_calloc(MATRIX_SIZE, MATRIX_SIZE);
  gsl_matrix *U = gsl_matrix_calloc(MATRIX_SIZE, MATRIX_SIZE);

  for (size_t i = 0; i < MATRIX_SIZE; i++) {
    gsl_matrix_set(L, i, i, 1.);
  }

}

int main(void) {
  printf("Gauss.\n\n");

  // Alloc matrix`s
  gsl_matrix_const_view matrix_1v = gsl_matrix_const_view_array((const double *) matrix_1, MATRIX_SIZE, MATRIX_SIZE);
  gsl_vector_const_view vector_1v = gsl_vector_const_view_array(vector_1, MATRIX_SIZE);

  const gsl_matrix *matrix1 = &matrix_1v.matrix;
  const gsl_vector *vector1 = &vector_1v.vector;

  ulDecomposition(matrix1, vector1);
  // gsl_matrix *orig_matrix = gsl_matrix_alloc(MATRIX_SIZE, MATRIX_SIZE);

  // gsl_vector *vector1 = gsl_vector_alloc(MATRIX_SIZE);
  // gsl_vector *orig_vector = gsl_vector_alloc(MATRIX_SIZE);

  // gsl_vector *result1 = gsl_vector_alloc(MATRIX_SIZE);
  // gsl_vector *result2 = gsl_vector_alloc(MATRIX_SIZE);

  // Init matrix

  // Init vector
}

void gaussSeidelSolve() {

}

max_stars_repo_name: unixs/numerical-la-matrix-forms

max_issues_repo_name: unixs/numerical-la-matrix-forms

max_forks_repo_name: unixs/numerical-la-matrix-forms


text:
/* Internal header file for cminpack, by Frederic Devernay. */
#ifndef __CMINPACKP_H__
#define __CMINPACKP_H__

#ifndef __CMINPACK_H__
#error "cminpackP.h in an internal cminpack header, and must be included after all other headers (including cminpack.h)"
#endif

#if (defined (USE_CBLAS) || defined (USE_LAPACK)) && !defined (__cminpack_double__) && !defined (__cminpack_float__)
#error "cminpack can use cblas and lapack only in double or single precision mode"
#endif

#ifdef USE_CBLAS
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif
#define __cminpack_enorm__(n,x) __cminpack_cblas__(nrm2)(n,x,1)
#else
#define __cminpack_enorm__(n,x) __cminpack_func__(enorm)(n,x)
#endif

#ifdef USE_LAPACK
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
#if defined(__LP64__) /* In LP64 match sizes with the 32 bit ABI */
typedef int 		__CLPK_integer;
typedef int 		__CLPK_logical;
typedef __CLPK_logical 	(*__CLPK_L_fp)();
typedef int 		__CLPK_ftnlen;
#else
typedef long int 	__CLPK_integer;
typedef long int 	__CLPK_logical;
typedef __CLPK_logical 	(*__CLPK_L_fp)();
typedef long int 	__CLPK_ftnlen;
#endif
int __cminpack_lapack__(lartg_)(
  __cminpack_real__ *f, __cminpack_real__ *g, __cminpack_real__ *cs,
  __cminpack_real__ *sn, __cminpack_real__ *r__);
int __cminpack_lapack__(geqp3_)(
  __CLPK_integer *m, __CLPK_integer *n, __cminpack_real__ *a, __CLPK_integer * lda,
  __CLPK_integer *jpvt, __cminpack_real__ *tau, __cminpack_real__ *work, __CLPK_integer *lwork,
  __CLPK_integer *info);
int __cminpack_lapack__(geqrf_)(
  __CLPK_integer *m, __CLPK_integer *n, __cminpack_real__ *a, __CLPK_integer * lda,
  __cminpack_real__ *tau, __cminpack_real__ *work, __CLPK_integer *lwork, __CLPK_integer *info);
#endif
#endif

#include "minpackP.h"

#endif /* !__CMINPACKP_H__ */

ext: h

max_stars_repo_name: ianormy/Pyto

max_issues_repo_name: ProjectZeroDays/Pyto

max_forks_repo_name: ProjectZeroDays/Pyto


text:
#include <pygsl/solver.h>
#include <gsl/gsl_min.h>

static const char min_f_type_name[] = "F-Minimizer";
static const char min_x_minimum_doc[]= "";
static const char min_x_lower_doc[]= "";
static const char min_x_upper_doc[]= "";
static const char min_f_minimum_doc[]= "";
static const char min_f_lower_doc[]= "";
static const char min_f_upper_doc[]= "";
static const char min_test_delta_doc[]= "";
static const char min_set_f_doc[]= "";

static const char min_init_brent_doc[]= "";
static const char min_init_goldensection_doc[]= "";

const char  * filename = __FILE__;
PyObject *module = NULL;

static PyObject*
PyGSL_min_f_minimum(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_f_minimum);
}

static PyObject*
PyGSL_min_f_lower(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_f_lower);
}

static PyObject*
PyGSL_min_f_upper(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_f_upper);
}

static PyObject*
PyGSL_min_x_minimum(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_x_minimum);
}

static PyObject*
PyGSL_min_x_lower(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_x_lower);
}

static PyObject*
PyGSL_min_x_upper(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_min_fminimizer_x_upper);
}

static PyObject*
PyGSL_min_solver_test_interval(PyGSL_solver * self, PyObject *args)
{
     double epsabs, epsrel;
     gsl_min_fminimizer *s = (gsl_min_fminimizer *) self->solver;
     if(!PyArg_ParseTuple(args, "dd", &epsabs, &epsrel))
	  return NULL;
     return PyInt_FromLong(gsl_min_test_interval(s->x_lower, s->x_upper, epsabs, epsrel));
}

static PyObject*
PyGSL_min_set_f(PyGSL_solver *self, PyObject *args, PyObject *kw)
{
     return PyGSL_solver_set_f(self, args, kw,  (void *)gsl_min_fminimizer_set, 0);
}

static PyMethodDef PyGSL_min_fmethods[] = {
     {"x_minimum", (PyCFunction)PyGSL_min_x_minimum,  METH_NOARGS, (char *)min_x_minimum_doc},
     {"x_lower",   (PyCFunction)PyGSL_min_x_lower,    METH_NOARGS, (char *)min_x_lower_doc},
     {"x_upper",   (PyCFunction)PyGSL_min_x_upper,    METH_NOARGS, (char *)min_x_upper_doc},
     {"f_minimum", (PyCFunction)PyGSL_min_f_minimum,  METH_NOARGS, (char *)min_f_minimum_doc},
     {"f_lower",   (PyCFunction)PyGSL_min_f_lower,    METH_NOARGS, (char *)min_f_lower_doc},
     {"f_upper",   (PyCFunction)PyGSL_min_f_upper,    METH_NOARGS, (char *)min_f_upper_doc},
     {"set",       (PyCFunction)PyGSL_min_set_f,    METH_VARARGS|METH_KEYWORDS, (char *)min_set_f_doc},
     {"test_interval",(PyCFunction)PyGSL_min_solver_test_interval, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

const struct _SolverStatic
min_solver_f = {{(void_m_t) gsl_min_fminimizer_free,
		 /* gsl_multimin_fminimizer_restart */  (void_m_t) NULL,
		 (name_m_t) gsl_min_fminimizer_name,
		 (int_m_t) gsl_min_fminimizer_iterate},
		1, PyGSL_min_fmethods, min_f_type_name};

static PyObject*
PyGSL_min_f_init(PyObject *self, PyObject *args,
		 const gsl_min_fminimizer_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_min_fminimizer_alloc,
			      &min_solver_f};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 0);
     FUNC_MESS_END();
     return tmp;
}

#define AMIN_F(name)                                                  \
static PyObject* PyGSL_min_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_min_f_init(self, args,  gsl_min_fminimizer_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}
AMIN_F(brent)
AMIN_F(goldensection)

static PyObject*
PyGSL_min_test_interval(PyObject * self, PyObject *args)
{
     double x_lower, x_upper, epsabs, epsrel;
     if(!PyArg_ParseTuple(args, "dddd", &x_lower, &x_upper, &epsabs, &epsrel))
	  return NULL;
     return PyInt_FromLong(gsl_min_test_interval(x_lower, x_upper, epsabs, epsrel));
}

static const char PyGSL_minimize_module_doc [] = "XXX Missing ";
static PyMethodDef mMethods[] = {
     /* solvers */
     {"brent",        PyGSL_min_init_brent, METH_NOARGS, (char *)min_init_brent_doc},
     {"goldensection",PyGSL_min_init_goldensection, METH_NOARGS, (char *)min_init_goldensection_doc},
     /* min methods */
     {"test_interval",   PyGSL_min_test_interval, METH_VARARGS, (char *)min_test_delta_doc},
     {NULL, NULL, 0, NULL}
};

void
initminimize(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("minimize", mMethods);
     import_pygsl_solver();
     assert(PyGSL_API);

     module = m;
     assert(m);
     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     if (!(item = PyString_FromString((char*)PyGSL_minimize_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }

     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();

 fail:
     FUNC_MESS("FAIL");
     return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include <cJSON.h>
#include <memcached/mcd_util-visibility.h>

#include <gsl/gsl>
#include <string>

namespace cb {
namespace breakpad {
/**
 * What information should breakpad minidumps contain?
 */
enum class Content {
    /**
     * Default content (threads+stack+env+arguments)
     */
    Default
};

/**
 * Settings for Breakpad crash catcher.
 */
struct MCD_UTIL_PUBLIC_API Settings {
    /**
     * Default constructor initialize the object to be in a disabled state
     */
    Settings() = default;

    /**
     * Initialize the Breakpad object from the specified JSON structure
     * which looks like:
     *
     *     {
     *         "enabled" : true,
     *         "minidump_dir" : "/var/crash",
     *         "content" : "default"
     *     }
     *
     * @param json The json to parse
     * @throws std::invalid_argument if the json dosn't look as expected
     */
    explicit Settings(gsl::not_null<const cJSON*> json);

    bool enabled{false};
    std::string minidump_dir;
    Content content{Content::Default};
};

} // namespace breakpad
} // namespace cb

MCD_UTIL_PUBLIC_API
std::string to_string(cb::breakpad::Content content);


text:
//
//  editor.hpp
//  PretendToWork
//
//  Created by tqtifnypmb on 14/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include "../crdt/Engine.h"
#include "../types.h"

#include <memory>
#include <gsl/gsl>
#include <map>

namespace brick
{

class Rope;
class Editor {
public:

    using DeltaList = std::vector<std::tuple<Revision, size_t, size_t>>;     // [(delta, begRow, endRow)]
    using SyncCb = std::function<void(const Editor::DeltaList& deltaList)>;

    Editor(size_t viewId, const detail::CodePointList& cplist, SyncCb syncCb);
    Editor(size_t viewId, SyncCb syncCb);

    template <class Converter>
    void insert(gsl::span<const char> bytes, size_t pos);
    void insert(const detail::CodePointList& cplist, size_t pos);
    void erase(Range range);
    void undo();

    void sync(Editor& editor);

    std::map<size_t, detail::CodePointList> region(size_t begLine, size_t endLine);

    void clearRevisions() {
        engine_.revisions().clear();
    }

private:
    Editor::DeltaList convertEngineDelta(const Engine::DeltaList& deltas);

    void updateLines(size_t pos, const detail::CodePointList& cplist);
    void updateLines(Range r);

    std::unique_ptr<Rope> rope_;
    Engine engine_;
    std::vector<size_t> linesIndex_;
    SyncCb sync_cb_;
};

template <class Converter>
void Editor::insert(gsl::span<const char> bytes, size_t pos) {
    insert(Converter::encode(bytes), pos);
}

}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick