text:
/******************************************************************************
 *                                                                            *
 * DECS.H                                                                     *
 *                                                                            *
 * GLOBAL MACROS, FUNCTION DEFINITIONS, INCLUDES, AND DECLARATIONS            *
 *                                                                            *
 ******************************************************************************/

#include <gsl/gsl_eigen.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_vector.h>
#include <limits.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include "constants.h"
#include "params.h"

/*******************************************************************************
      PRE-PROCESSOR MAGIC :
*******************************************************************************/
#define XSTR(x) STR(x)
#define STR(x) #x

/*******************************************************************************
      CONSTANTS :
*******************************************************************************/

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923132169164
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.4142135623730950488016887242
#endif

/*******************************************************************************
      COMPILE-TIME PARAMETERS :
*******************************************************************************/

// Number of active zones on each MPI process
#define N1 (N1TOT / N1CPU)
#define N2 (N2TOT / N2CPU)
#define N3 (N3TOT / N3CPU)

// Max size for 1D slice is NMAX
#define N12 (N1 > N2 ? N1 : N2)
#define NMAX (N12 > N3 ? N12 : N3)

#define NDIM (4) // Number of total dimensions
#define NPG (5)  // Number of positions on grid for grid functions
#define NG (3)   // Number of ghost zones
#define N123G                      \
  ((N1 + 2 * NG) * (N2 + 2 * NG) * \
      (N3 + 2 * NG)) // total number of cells on a cpu

// Fixup parameters
// rho_flor = RHOMIN / pow(r, -FLOOR_POWER)
#define FLR_POWER1 (2.)
#define FLR_POWER2 (2.)
#define FLR_R0 (50.) // TODO change this?
#define RHOMINLIMIT (1.e-17)
#define UUMINLIMIT (1.e-3 * RHOMINLIMIT)
#define RHOMIN (1.e-5)
#define UUMIN (1.e-3 * RHOMIN)
#define BSQORHOMAX (50.)
#define BSQOUMAX (2500.)
#define RHOEPS (0.9)
#define UORHOMAX (50.)
#define COLD_FLOORS (1)
#define ATM_THRESH (0.1)

// Numerical convenience to represent a small (<< 1) non-zero quantity
#define SMALL (1.e-20)
#define ABSURDLY_LARGE (1.e100)

// Maximum value of gamma, the Lorentz factor
#define GAMMAMAX (50.)

// Maximum fractional increase in timestep per timestep
#define SAFE (1.3)

// Superphoton diagnostics
#define MAXNSCATT (3)

// Chosen to give some wiggle room
// in scattering stability criterion
#define SCATT_BIAS_SAFETY (0.75) // (0.9/RAD_SCATT_TYPES)

//#define NUMIN     (1.e10)
//#define NUMAX     (1.e25)
//#define NU_BINS_SPEC    (200)
//#define NTH (8)
//#define NPHI (8)

// Whether to move polar axis slightly off of coordinate singularity
#define COORDSINGFIX 1
#define SINGSMALL (1.E-20)

// Whether to record a representative sample of superphoton positions
#define TRACK_PH (0)

// I/O format strings
#define FMT_DBL_OUT "%28.18e"
#define FMT_INT_OUT "%10d"
#define STRLEN (2048)

// Reconstruction algorithms
#define LINEAR (0)
#define PPM (1)
#define WENO (2)
#define MP5 (3)

// Primitive and conserved variables
#define RHO (0)
#define UU (1)
#define U1 (2)
#define U2 (3)
#define U3 (4)
#define B1 (5)
#define B2 (6)
#define B3 (7)
#define NVAR_BASE (B3 + 1)

// Passive variables (if present)
#define PASSIVE_START (NVAR_BASE)
#define PASSIVE_STOP (NVAR_BASE + NVAR_PASSIVE)
#define PASSTYPE_INTRINSIC (0)
#define PASSTYPE_NUMBER (1)
#if EOS == EOS_TYPE_TABLE
#define YE (PASSIVE_START)
#define YE_EM (YE + 1)
#if METRIC == MKS
#define ATM (YE + 2)
#define IS_ATM (0.0)
#define NOT_ATM (1.0)
#endif // METRIC
#endif // EOS_TYPE_TABLE

#if ELECTRONS
#define KEL (B3 + NVAR_PASSIVE + 1)
#define KTOT (B3 + NVAR_PASSIVE + 2)
#define NVAR (NVAR_BASE + NVAR_PASSIVE + 2)
#else
#define NVAR (NVAR_BASE + NVAR_PASSIVE)
#endif

// NETRINOS
#define RATYPE_NONE (0)
#define RADTYPE_LIGHT (1)
#define RADTYPE_NEUTRINOS (2)
#define RAD_TYPE_START (0)
#define TYPE_TRACER (-1)
#if RADIATION == RADTYPE_NEUTRINOS
#define RAD_NUM_TYPES (3)
#define NU_ELECTRON (0)
#define ANTINU_ELECTRON (1)
#define NU_HEAVY (2)
#if MULTISCATT_TEST
#define RAD_SCATT_TYPES (3)
#else
#define RAD_SCATT_TYPES (4) // TODO: Should be 5, including electrons
#define RSCATT_TYPE_P (0)
#define RSCATT_TYPE_N (1)
#define RSCATT_TYPE_A (2)
#define RSCATT_TYPE_ALPHA (3)
#define RSCATT_TYPE_E (4) // TOOD: implement me
#endif
#define NRADCOMP (2)
#define RADG_YE (4)
#define RADG_YE_EM (5)
#elif RADIATION == RADTYPE_LIGHT
#define RAD_NUM_TYPES (1)
#define RAD_SCATT_TYPES (1)
#define NRADCOMP (0)
#define PHOTON (0)
#else
#define RAD_NUM_TYPES (0)
#define RAD_SCATT_TYPES (0)
#define NRADCOMP (0)
#endif
#define DLEPTON_THRESH (1e-10)

// EOS
#define EOS_TYPE_GAMMA (0)
#define EOS_TYPE_POLYTROPE (1)
#define EOS_TYPE_TABLE (2)
#if EOS == EOS_TYPE_GAMMA
#define EOS_NUM_EXTRA (0)
#define POLYTROPE_FALLBACK (0)
#elif EPS == EOS_TYPE_POLYTROPE
#define EOS_NUM_EXTRA (0)
#define POLYTROPE_FALLBACK (0)
#elif EOS == EOS_TYPE_TABLE
#if GAMMA_FALLBACK
#define POLYTROPE_FALLBACK (0)
// #define POLYTROPE_FALLBACK (1)
#else
#define POLYTROPE_FALLBACK (1)
#endif // GAMMA_FALLBACK
#define EOS_NUM_EXTRA (3)
#define EOS_LRHO (0)
#define EOS_LT (1)
#define EOS_YE (2)
// mass fractions
#define NUM_MASS_FRACTIONS (4)
#define MF_XA (0)
#define MF_XH (1)
#define MF_XN (2)
#define MF_XP (3)
#endif // EOS

// Centering of grid functions
#define FACE1 (0)
#define FACE2 (1)
#define FACE3 (2)
#define CORN (3)
#define CENT (4)
#define FACESTART (FACE1)
#define FACEEND (FACE3 + 1)
#define PGSTART (FACE1)
#define PGEND (NPG)

// Slope limiter
#define MC (0)
#define VANL (1)
#define MINM (2)

// Fluid and radiation boundaries
#define BC_OUTFLOW (0)
#define BC_PERIODIC (1)
#define BC_POLAR (2)
#define BC_PROB (3)
#define BC_ESCAPE (4)
#define BC_CAMERA (5)
#define BC_EQUILIB (6)

// Metric
#define MINKOWSKI (0)
#define SPHERICAL (1)
#define MKS (2)

// Diagnostic calls
#define DIAG_INIT (0)
#define DIAG_DUMP (1)
#define DIAG_LOG (2)
#define DIAG_FINAL (3)

// Types of restarts
#define RESTART_TEMP (0)
#define RESTART_PERM (1)

// Failure modes
#define FAIL_UTOPRIM (0)
#define FAIL_VCHAR_DISCR (1)
#define FAIL_COEFF_NEG (2)
#define FAIL_COEFF_SUP (3)
#define FAIL_GAMMA (4)
#define FAIL_METRIC (5)

// Geodesic integration and interpolation
#define PUSH_FAIL (0)
#define PUSH_SUCCESS (1)
#define SPH_INTERP_FAIL (0)
#define SPH_INTERP_SUCCESS (1)

// Root finding
#define ROOT_SUCCESS (1)
#define ROOT_FAIL (0)
#define FCOUNT_NBINS (6)
#define FCOUNT_MORE (FCOUNT_NBINS - 1)

// Timers
#define TIMER_UPDATE (0)
#define TIMER_FLUXCALC (1)
#define TIMER_FIXUP (2)
#define TIMER_BOUND (3)
#define TIMER_DIAG (4)
#define TIMER_OUT (5)
#define TIMER_ELECTRON (6)
#define TIMER_MAKE (7)
#define TIMER_PUSH (8)
#define TIMER_INTERACT (9)
#define TIMER_MICRO (10)
#define TIMER_ALL (11)
#define NUM_TIMERS (12)

// Units
#define NEED_UNITS (RADIATION || EOS == EOS_TYPE_TABLE)

/*******************************************************************************
    GLOBAL ARRAYS
*******************************************************************************/
typedef double grid_prim_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG][NVAR];
typedef double grid_double_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG];
typedef double grid_radtype_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG]
                                [RAD_NUM_TYPES];
typedef double grid_fourvector_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG]
                                   [NDIM];
typedef double grid_radg_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG]
                             [NDIM + NRADCOMP];
typedef double grid_tensor_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG][NDIM]
                               [NDIM];
typedef int    grid_int_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG];
typedef double grid_eosvar_type[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG]
                               [EOS_NUM_EXTRA];
typedef void (*passive_init_ftype)(int, int, int, double *, double *);
typedef double (*hc_ftype)(double, double);

extern grid_prim_type       P;     // Primitive variables
extern grid_prim_type       F1;    // X1 fluxes
extern grid_prim_type       F2;    // X2 fluxes
extern grid_prim_type       F3;    // X3 fluxes
extern grid_prim_type       Ph;    // Half-step primitives
extern grid_prim_type       Psave; // Half-step primitives
extern grid_int_type        pflag; // Failure points
extern grid_int_type        fail_save;
extern grid_int_type        fixup_required;
extern grid_fourvector_type jcon;
extern grid_eosvar_type     extra; // extra variables needed by EOS
#if RADIATION
extern grid_radg_type   radG;     // Radiation four-force
extern grid_radg_type   radG_int; // ...integrated
extern grid_radg_type   radG_buf; // ...buffer for communication
extern grid_tensor_type Rmunu;    // Radiation stress-energy tensor
extern grid_int_type    Nsph;
extern grid_double_type nph;

extern struct of_photon **photon_lists;
extern struct of_photon **photon_mpi_lists;

#if DIAGNOSTICS_USE_RADTYPES
#define NULNU_IDX0 (RAD_NUM_TYPES)
#else
#define NULNU_IDX0 (MAXNSCATT + 1)
#endif
extern double nuLnu[NULNU_IDX0][NTH][NPHI][NU_BINS_SPEC];

#if RADIATION == RADTYPE_NEUTRINOS
extern double rad_type_counts[RAD_NUM_TYPES];
extern double lepton_tot, lepton_last, dlepton_tot, dlepton_perc;
extern double lepton_gas, lepton_rad;
extern double lepton_lost, lepton_lost_step, lepton_lost_local;
#endif

extern double Jrad[MAXNSCATT + 2][N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG];
extern double Jrad_buf[MAXNSCATT + 2][N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG];
extern double dtau_avg[RAD_SCATT_TYPES + 1][N1 + 2 * NG][N2 + 2 * NG]
                      [N3 + 2 * NG];
extern double en_int_avg[RAD_SCATT_TYPES + 1][N1 + 2 * NG][N2 + 2 * NG]
                        [N3 + 2 * NG];
extern grid_int_type     Nem, Nabs, Nsc;
extern grid_radtype_type Nem_phys, Nabs_phys, radtype_buf;

extern grid_int_type    Nsuper;
extern grid_double_type Esuper;
extern grid_prim_type   psupersave;
#endif // RADIATION

// Default initialization is 0, which in this case is
// PASSTYPE_INTRINSIC, i.e., volume densities.
extern int                passive_type[NVAR_PASSIVE];
extern char               passive_name[NVAR_PASSIVE][STRLEN];
extern passive_init_ftype do_passive_fixup[NVAR_PASSIVE];

#if ELECTRONS
extern grid_double_type Qvisc, Qcoul;
#endif // ELECTRONS

/*******************************************************************************
    GLOBAL VARIABLES SECTION
*******************************************************************************/
// Command line arguments
extern char outputdir[STRLEN], dumpdir[STRLEN], restartdir[STRLEN];
extern char xmfdir[STRLEN];
#if RADIATION && TRACERS
extern char tracerdir[STRLEN];
#endif
extern char init_from_grmhd[STRLEN];
extern char metric[STRLEN], reconstruction[STRLEN];
extern char eos[STRLEN], nulnutype[STRLEN];
extern int  tracers;

// path for eos
#if EOS == EOS_TYPE_TABLE
extern char eospath[STRLEN];
#endif

// opacity table paths
#if RADIATION
#if (RADIATION == RADTYPE_NEUTRINOS)
#if BURROWS_OPACITIES
extern char opac_param_file[STRLEN];
extern char opac_file[STRLEN];
#endif // BURROWS
#if HDF5_OPACITIES
extern char opac_file[STRLEN];
#endif // hdf5
#endif // neutrinos
#endif // radiation

// Physics parameters
extern double a;
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
extern double gam;
#endif
#if EOS == EOS_TYPE_POLYTROPE
extern double poly_K, poly_gam;
#endif
#if POLYTROPE_FALLBACK
extern double rho_poly_thresh;
#endif
extern double M_unit;
extern double Reh;
extern double Risco;
#if NEED_UNITS
extern double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit;
#endif
#if EOS == EOS_TYPE_TABLE
extern double TEMP_unit;
#endif
#if RADIATION
extern double Ne_unit, Thetae_unit, kphys_to_num;
extern double tp_over_te, thetae_max, sigma_max, kdotk_tol;
#endif

// Numerical parameters
extern double Rin, Rout, Rout_vis, hslope;
extern double poly_norm, poly_xt, poly_alpha, mks_smooth;
#if RADIATION
extern double Rout_rad;
extern double nph_per_proc;
extern double tune_emiss, t_tune_emiss, dt_tune_emiss;
extern double tune_scatt, t_tune_scatt, dt_tune_scatt;
extern int    made_tune_proc, abs_tune_proc, scatt_tune_proc;
extern double numin, numax;
extern double kappa;
extern double startx_rad[NDIM], stopx_rad[NDIM];
extern double wgtC;
extern int    step_made, step_abs, step_scatt, step_lost, step_rec, step_tot;
extern int    tracer_tot, tracer_tot_all, ntcr_per_proc;
extern int    step_sent, step_rcvd, step_fail;
extern int    step_made_all, step_abs_all, step_scatt_all, step_lost_all;
extern int    step_rec_all, step_sent_all, step_rcvd_all, step_tot_all;
extern int    step_fail_all;
extern int    step_tot_max, step_tot_min;
extern double load_imbalance;
extern double Nph_to_track;
extern double sim_vol;
#if FLATEMISS
extern double cnu_flat;
#endif // FLATEMISS
#if MULTISCATT_TEST
extern double ms_theta_nu0, ms_delta0;
#endif // MULTISCATT_TEST
#endif // RADIATION
#if ELECTRONS
extern double tptemin, tptemax;
#endif
extern double cour;
#if RADIATION
extern double cour_cool;
#endif
extern double dV, dx[NDIM], startx[NDIM], stopx[NDIM], startx_proc[NDIM],
    stopx_proc[NDIM];
extern double x1Min, x1Max, x2Min, x2Max, x3Min, x3Max;
extern double dt, dtsave;
extern double t, tf;
extern int    nstep;
extern int    is_restart;

// Output parameters
extern double DTd;
extern double DTl;
extern double DTr;
extern int    DNr;
extern int    DTp;
extern int    DTf;
extern double DTw;
extern int    dump_cnt;
extern int    rdump_cnt;
extern double tdump, trestart, tlog;
extern double root_fcount[FCOUNT_NBINS];

// Global flags
extern int failed;
extern int lim;

// Diagnostics
extern double mdot, mdot_eh;
extern double edot, edot_eh;
extern double ldot, ldot_eh;
extern int    icurr, jcurr, kcurr;

// Parallelism
extern int nthreads;

// Electrons
#if ELECTRONS
extern double game, gamp;
extern double fel0;
#endif

// Set global variables that indicate current local metric, etc.
struct of_geom {
  double gcon[NDIM][NDIM];
  double gcov[NDIM][NDIM];
  double g;
  double alpha;
};

struct of_state {
  double ucon[NDIM];
  double ucov[NDIM];
  double bcon[NDIM];
  double bcov[NDIM];
};

#if RADIATION
// WARNING: if you change struct_of_photon, be sure to change
// the the photon hdf5 and MPI  types in io.c and mpi.c
#define NSUP 3
struct of_photon {
  // NSUP >=3 X^{\mu}, K^{\mu}, K_{\mu} so photon data always available anywhere
  // between n and n+1
  double X[NSUP][NDIM];
  double Kcov[NSUP][NDIM];
  double Kcon[NSUP][NDIM];
  double w;
  double KdotKprev;
  // radiation type. For neutrinos, flavor. Always active.
  // Not always important.
  // TODO: make sure to always set type when it's needed.
  int               type;
  int               nscatt;
  int               origin[NDIM];
  double            t0;
  int               is_tracked;
  struct of_photon *next;
};

#define PH_ELEM (11)
struct of_track_photon {
  double X1;
  double X2;
  double X3;
  int    nscatt;
};
struct of_microphysics {
#if RADIATION == RADTYPE_NEUTRINOS
  double rho;
  double T;
  double Ye;
  double Abar, Zbar;
  double Xi[NUM_MASS_FRACTIONS];
#if BURROWS_OPACITIES || HDF5_OPACITIES
  double jgrp[RAD_NUM_TYPES][NU_BINS + 1];
  double alphagrp[RAD_NUM_TYPES][NU_BINS + 1];
#endif
#else // RADTYPE_LIGHT
  double Thetae;
  double Ne;
#endif
  // needed for Bk-angle (which isn't really needed for neutrinos, but
  // this is the easy engineering solution)
  double B;
};
typedef double (*dsdom_ftype)(
    double, double, const struct of_microphysics *, int, int);
#endif // RADIATION

#if EOS == EOS_TYPE_TABLE
struct of_adiabat {
  double  s, ye;
  double  lrho_min, lrho_max;
  int     imin, imax;
  double  hm1_min, hm1_max;
  double *lT;
};
#endif

#if EOS == EOS_TYPE_TABLE || HDF_OPACITIES || BURROWS_OPACITIES
struct of_tablebounds {
  int    Nrho, NT, NYe;
  double lrho_min, lrho_max, dlrho;
  double lT_min, lT_max, dlT;
  double Ye_min, Ye_max, dYe;
};
#endif

// More grid functions. Axisymmetry assumed.
extern double         conn[N1 + 2 * NG][N2 + 2 * NG][NDIM][NDIM][NDIM];
extern struct of_geom ggeom[N1 + 2 * NG][N2 + 2 * NG][NPG];
#if RADIATION
// extern double dt_light, dt_light_min;
extern double                 dt_light[N1 + 2 * NG][N2 + 2 * NG], dt_light_min;
extern struct of_microphysics m_grd[N1 + 2 * NG][N2 + 2 * NG][N3 + 2 * NG];
extern grid_fourvector_type   Ucon_grd, Ucov_grd, Bcon_grd, Bcov_grd;
#endif

// MPI-specific stuff
extern int global_start[NDIM];
extern int global_stop[NDIM];

/*******************************************************************************
    MACROS
*******************************************************************************/
#define ILOOP for (int i = 0 + NG; i < N1 + NG; i++)
#define ILOOPALL for (int i = 0; i < N1 + 2 * NG; i++)
#define JLOOP for (int j = 0 + NG; j < N2 + NG; j++)
#define JLOOPALL for (int j = 0; j < N2 + 2 * NG; j++)
#define KLOOP for (int k = 0 + NG; k < N3 + NG; k++)
#define KLOOPALL for (int k = 0; k < N3 + 2 * NG; k++)
#define ZLOOP                              \
  for (int i = 0 + NG; i < N1 + NG; i++)   \
    for (int j = 0 + NG; j < N2 + NG; j++) \
      for (int k = 0 + NG; k < N3 + NG; k++)
#define ZLOOPALL ILOOPALL JLOOPALL KLOOPALL
#define ISLOOP(istart, istop) for (int i = istart + NG; i <= istop + NG; i++)
#define JSLOOP(jstart, jstop) for (int j = jstart + NG; j <= jstop + NG; j++)
#define KSLOOP(kstart, kstop) for (int k = kstart + NG; k <= kstop + NG; k++)
#define ZSLOOP(istart, istop, jstart, jstop, kstart, kstop) \
  for (int i = istart + NG; i <= istop + NG; i++)           \
    for (int j = jstart + NG; j <= jstop + NG; j++)         \
      for (int k = kstart + NG; k <= kstop + NG; k++)
// Loop over faces
#define FACELOOP for (int face = FACESTART; face < FACEEND; face++)
// Loop over all locations
#define LOCLOOP for (int loc = PGSTART; loc < PGEND; loc++)
// Loop over primitive variables
#define PLOOP for (int ip = 0; ip < NVAR; ip++)
#define BASELOOP for (int ip = 0; ip < NVAR_BASE; ip++)

// Loop over spacetime indices
#define DLOOP1 for (int mu = 0; mu < NDIM; mu++)
#define DLOOP2                      \
  for (int mu = 0; mu < NDIM; mu++) \
    for (int nu = 0; nu < NDIM; nu++)
#define DLOOP3                        \
  for (int mu = 0; mu < NDIM; mu++)   \
    for (int nu = 0; nu < NDIM; nu++) \
      for (int sigma = 0; sigma < NDIM; sigma++)
// spacelike indices
#define SDLOOP for (int mu = 1; mu < NDIM; mu++)

// Loop over extra variables
// TODO: Figure out how to make this conditionally defined. ~JMM
#define EOS_ELOOP for (int e = 0; e < EOS_NUM_EXTRA; e++)

// Loop over passive scalars
#define PASSLOOP for (int ipass = PASSIVE_START; ipass < PASSIVE_STOP; ipass++)
#define PASSELEM(ipass) (ipass - PASSIVE_START)
#define PASSTYPE(ipass) (passive_type[PASSELEM(ipass)])
#define PASSNAME(ipass) (passive_name[PASSELEM(ipass)])

#define TYPELOOP for (int itp = RAD_TYPE_START; itp < RAD_NUM_TYPES; itp++)
#define SCATTLOOP for (int iscatt = 0; iscatt < RAD_SCATT_TYPES; iscatt++)
#define JRADLOOP for (int n = 0; n < MAXNSCATT + 2; n++)
#define NULOOP for (int inu = 0; inu < NU_BINS + 1; inu++)

#define MY_MIN(fval1, fval2) (((fval1) < (fval2)) ? (fval1) : (fval2))
#define MY_MAX(fval1, fval2) (((fval1) > (fval2)) ? (fval1) : (fval2))
#define MY_SIGN(fval) (((fval) < 0.) ? -1. : 1.)

#define delta(i, j) ((i == j) ? 1. : 0.)
#define dot(a, b) (a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3])

/*******************************************************************************
    FUNCTION DECLARATIONS
*******************************************************************************/

// bounds.c
void bound_prim(grid_prim_type prim);
void fix_flux(grid_prim_type F1, grid_prim_type F2, grid_prim_type F3);
#if RADIATION
void bound_superphotons(grid_prim_type P, double t, double dt);
void polar_fix(double X[NDIM], struct of_photon *ph);
int  rad_error_check(struct of_photon **ph);
// int rad_boundary_transport(double X[NDIM], double Kcov[NDIM]);
void bound_rad_transport(
    double X[NDIM], struct of_photon *ph, int is_transported);
int bound_rad_isactive(double X[NDIM], struct of_photon *ph);
int rad_mpi_transport(struct of_photon **ph, struct of_photon **prev,
    struct of_photon **head, double X[NDIM], int active);
#endif

// coord.c
void   coord(int i, int j, int k, int loc, double X[NDIM]);
double r_of_X(const double X[NDIM]);
double th_of_X(const double X[NDIM]);
void   jac_harm_to_bl(
      const double X[NDIM], double Jcov[NDIM][NDIM], double Jcon[NDIM][NDIM]);
void jac_bl_to_cart(
    const double X[NDIM], double Jcov[NDIM][NDIM], double Jcon[NDIM][NDIM]);
void jac_harm_to_cart(
    const double X[NDIM], double Jcov[NDIM][NDIM], double Jcon[NDIM][NDIM]);
void set_dxdX(double X[NDIM], double dxdX[NDIM][NDIM]);
void bl_coord(const double X[NDIM], double *r, double *th);
int  bl_i_of_r(double r);
void cart_coord(const double X[NDIM], double Xcart[NDIM]);
void set_gcov(double X[NDIM], double gcov[NDIM][NDIM]);
void set_metric(double X[NDIM], struct of_geom *g);
void set_points();
void zero_arrays(void);
void set_grid(void);
void cart_to_sph(const double X[NDIM], double *r, double *th, double *ph);

// current.c
void current_calc();

// diag.c
void   reset_log_variables();
void   reset_dump_variables();
void   diag(int call_code);
void   fail(int fail_type);
void   area_map(int i, int j, int k, grid_prim_type prim);
void   diag_flux(grid_prim_type F1, grid_prim_type F2, grid_prim_type F3);
double flux_ct_divb(int i, int j, int k);
#if RADIATION
void record_superphoton(double X[NDIM], struct of_photon *ph);
void bin_all_superphotons();
void report_load_imbalance();
#if RADIATION == RADTYPE_NEUTRINOS
void print_rad_types();
void count_leptons(grid_prim_type P, double dt, int nstep);
#endif
#endif

// electrons.c
#if ELECTRONS
void init_electrons();
void heat_electrons(
    grid_prim_type Pi, grid_prim_type Ps, grid_prim_type Pf, double Dt);
double get_fel(int i, int j, int k, double p[NVAR]);
#if RADIATION
void coulomb(
    grid_prim_type Pi, grid_prim_type Ps, grid_prim_type Pf, double Dt);
void apply_rad_force_e(
    grid_prim_type Prh, grid_prim_type Pr, grid_radg_type radG, double Dt);
#endif // RADIATION
void fixup_electrons(grid_prim_type p);
#endif

// emissivity.c
#if RADIATION
double jnu(double nu, int type, const struct of_microphysics *m, double theta);
double Jnu(double nu, int type, const struct of_microphysics *m);
double get_J(struct of_microphysics *m);
void   init_emissivity();
#endif

// eos.c
void   init_EOS();
double EOS_bad_eos_error();
double EOS_get_gamma(const double *extra);
double EOS_pressure_rho0_u(double rho, double u, const double *extra);
double EOS_pressure_rho0_w(double rho, double w, double gamma,
    const struct of_geom *geom, double *extra);
double EOS_enthalpy_rho0_u(double rho, double u, const double *extra);
double EOS_sound_speed_rho0_u(double rho, double u, const double *extra);
void   EOS_set_floors(double scale, double rho, double u, double bsq,
      double *rhoflr, double *uflr, const double *extra);
double EOS_entropy_rho0_u(double rho, double u, const double *extra);
double EOS_adiabatic_constant(double rho, double u, const double *extra);
double EOS_temperature(double rho, double u, const double *extra);
double EOS_u_press(double press, double rho, double *extra);
#if RADIATION
double EOS_u_N_Theta(double rho, double N, double Theta, double *extra);
double EOS_pressure_N_Theta(double N, double Theta);
double EOS_Theta_unit();
#endif // RADIATION

// eos_gamma.c
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
double EOS_Gamma_pressure_rho0_u(double rho, double u);
double EOS_Gamma_pressure_rho0_w(double rho, double w);
double EOS_Gamma_entropy_rho0_u(double rho, double u);
double EOS_Gamma_enthalpy_rho0_u(double rho, double u);
double EOS_Gamma_adiabatic_constant_rho0_u(double rho, double u);
double EOS_Gamma_sound_speed_rho0_u(double rho, double u);
void   EOS_Gamma_set_floors(double scale, double rho, double u, double bsq,
      double *rhoflr, double *uflr);
double EOS_Gamma_rho_floor(double scale, double bsq);
double EOS_Gamma_u_floor(double scale, double bsq);
double EOS_Gamma_u_scale(double rho);
double EOS_Gamma_u_press(double press);
double EOS_Gamma_temp(double rho, double u);
#if RADIATION
double EOS_Gamma_Theta_unit();
#endif // RADIATION
#endif // EOS_TYPE_GAMMA

// eos_poly.c
#if EOS == EOS_TYPE_POLYTROPE || POLYTROPE_FALLBACK
double EOS_Poly_pressure_rho0_u(double rho, double u, double K, double Gam);
double EOS_Poly_pressure_rho0_w(double rho, double w, double K, double Gam);
double EOS_Poly_enthalpy_rho0_u(double rho, double u, double K, double Gam);
double EOS_Poly_entropy_rho0_u(double rho, double u, double K, double Gam);
double EOS_Poly_sound_speed_rho0_u(double rho, double u, double K, double Gam);
void   EOS_Poly_set_floors(double scale, double rho, double u, double bsq,
      double *rhoflr, double *uflr);
double EOS_Poly_rho_floor(double scale, double bsq);
double EOS_Poly_u_floor(double scale, double bsq);
double EOS_Poly_adiabatic_constant(double rho, double u, double K, double Gam);
#endif // EOS_TYPE_POLY

// eos_stellar_collapse.c
#if EOS == EOS_TYPE_TABLE
void   EOS_SC_init(char *name);
void   EOS_SC_fill(double *restrict p, double *restrict eos);
double EOS_SC_pressure_rho0_u(double lrho, double lT, double ye);
double EOS_SC_pressure_rho0_w(double rho, double w, double ye, double *lTold);
double EOS_SC_specific_enthalpy_rho0_u(double lrho, double lT, double ye);
double EOS_SC_sound_speed(double lrho, double lT, double ye);
double EOS_SC_entropy(double lrho, double lT, double ye);
double EOS_SC_gamma(double lrho, double lT, double ye);
double EOS_SC_temperature(double lT);
double EOS_SC_get_u_of_T(double rho, double T, double ye);
double EOS_SC_u_press(double press, double rho, double ye, double *lTold);
void EOS_SC_mass_fractions(double Xi[NUM_MASS_FRACTIONS], const double *extra);
void EOS_SC_avg_ions(double *Abar, double *Zbar, const double *extra);
void EOS_SC_set_floors(double scale, double rho, double u, double ye,
    double bsqr, double *rhoflr, double *uflr);
double EOS_SC_rho_floor(double scale, double bsq);
double EOS_SC_u_floor(double scale, double bsq, double ye);
double EOS_SC_get_min_lrho();
double EOS_SC_get_min_rho();
double EOS_SC_get_min_lT();
double EOS_SC_get_minu(double rho, double ye, double scale);
void   EOS_SC_get_polytrope(
      double lrho, double lT, double ye, double *poly_K, double *poly_gamma);
double EOS_SC_hm1_min_adiabat(const struct of_adiabat *a);
int    EOS_SC_find_adiabat_1d(double s, double ye, struct of_adiabat *a);
void   EOS_SC_print_adiabat(const struct of_adiabat *a);
void   EOS_SC_adiabat_free(struct of_adiabat *a);
void   EOS_SC_isoentropy_hm1(double hm1, const struct of_adiabat *a,
      double *lrho_guess, double *rho, double *u);
void   do_ye_fixup(
      int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]);
void do_ye_em_fixup(
    int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]);
void do_atm_fixup(
    int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]);
void EOS_SC_get_bounds(struct of_tablebounds *b);
void EOS_SC_print_table_mins();
#endif // EOS_TYPE_TABLE

// estimate_thetae.c
#if RADIATION
double get_Thetae_est(int i, int j, int k);
void   estimate_Thetae(
      grid_prim_type P, grid_eosvar_type extra, double t, double dt);
#endif

// fixup.c
void   fixup(grid_prim_type pv, grid_eosvar_type extra);
double get_scale(int i, int j, int k);
void   fixup1zone(
      int i, int j, int k, double pv[NVAR], double extra[EOS_NUM_EXTRA]);
void fixup_utoprim(grid_prim_type pv, grid_eosvar_type extra);

// interact.c
#if RADIATION
void interact(grid_prim_type P, grid_eosvar_type extra, double t, double dt);
#endif

// input.c
void init_params(char *pfname);
void set_param(char *key, void *data);

// io.c
// void set_core_params();
// void set_param(char *key, void *data);
// void read_params(char *pfname);
void init_io();
void init_fluid_restart();
#if RADIATION
void track_ph();
#if TRACERS
void dump_tracers();
#endif
#endif
void dump_grid();
void dump();
void restart_write(int restart_type);
void restart_read(char *fname);
int  restart_init();

// make_superphotons.c
#if RADIATION
void make_superphotons(
    grid_prim_type Prad, grid_eosvar_type extra, double t, double dt);
void set_weight(grid_prim_type Prad, grid_eosvar_type extra);
void get_dnz(grid_prim_type Prad, grid_eosvar_type extra);
#endif

// metric.c
double gcon_func(double lgcov[][NDIM], double lgcon[][NDIM]);
void   conn_func(double *X, struct of_geom *geom, double conn[][NDIM][NDIM]);
void   lower(double ucon[NDIM], double gcov[NDIM][NDIM], double ucov[NDIM]);
void   raise(double ucov[NDIM], double gcon[NDIM][NDIM], double ucon[NDIM]);
struct of_geom *get_geometry(int ii, int jj, int kk, int loc);
void            blgset(int i, int j, struct of_geom *geom);
double          bl_gdet_func(double r, double th);
void            bl_gcov_func(double r, double th, double gcov[][NDIM]);
void            bl_gcon_func(double r, double th, double gcon[][NDIM]);
double MINOR(double m[16], int r0, int r1, int r2, int c0, int c1, int c2);
void   adjoint(double m[16], double adjOut[16]);
double determinant(double m[16]);
double invert(double *m, double *invOut);

// mpi.c
void init_mpi();
void sync_mpi_boundaries_X1L(grid_prim_type Pr);
void sync_mpi_boundaries_X1R(grid_prim_type Pr);
void sync_mpi_boundaries_X2L(grid_prim_type Pr);
void sync_mpi_boundaries_X2R(grid_prim_type Pr);
void sync_mpi_boundaries_X3L(grid_prim_type Pr);
void sync_mpi_boundaries_X3R(grid_prim_type Pr);
#if RADIATION
void sync_radG();
void sync_Jrad();
void sync_radtype_vec(grid_radtype_type v);
void sync_mpi_photons(
    struct of_photon **ph_mpi, grid_prim_type P, double t, double dt);
void mpi_reduce_nuLnu();
#endif
int    mpi_nprocs();
double mpi_max(double f);
double mpi_min(double f);
int    mpi_max_int(int f);
double mpi_reduce(double f);
int    mpi_reduce_int(int f);
void   mpi_dbl_allreduce_array(double *A, int size);
int    mpi_io_proc();
void   mpi_int_broadcast(int *val);
void   mpi_int_broadcast_array(int *val, int size);
void   mpi_int_broadcast_proc(int *val, int root);
void   mpi_dbl_broadcast(double *val);
void   mpi_allgather_int1(int *buffer, int element);
int    mpi_accumulate_int(int my_val);
double mpi_io_reduce(double val);
double mpi_io_max(double val);
int    mpi_myrank();
void   mpi_sync_output();
void   mpi_barrier();
int    mpi_is_periodic(int dir);

// opac_emis_neutrino.c
#if RADIATION
#if RADIATION == RADTYPE_NEUTRINOS && BURROWS_OPACITIES
void   init_opac_emis_burrows();
void   fill_opac_emis_burrows(struct of_microphysics *m);
double jnu_burrows(double nu, int type, const struct of_microphysics *m);
double Jnu_burrows(double nu, int type, const struct of_microphysics *m);
double int_jnudnudOmega_burrows(const struct of_microphysics *m);
double alpha_nu_burrows(double nu, int type, const struct of_microphysics *m);
void   test_opac_emis();
#if EOS == EOS_TYPE_TABLE
void opac_emis_to_hdf(const char *name, const struct of_tablebounds *b);
#endif // EOS
#endif // Burrows opacities
#if RADIATION == RADTYPE_NEUTRINOS && HDF5_OPACITIES
void   init_opac_emis_hdf(char *name);
void   fill_opac_emis_hdf(struct of_microphysics *m);
double jnu_hdf(double nu, int type, const struct of_microphysics *m);
double Jnu_hdf(double nu, int type, const struct of_microphysics *m);
double int_jnudnudOmega_hdf(const struct of_microphysics *m);
double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m);
#endif // HDF opacities
#endif // RADIATION

// passive.c
//#if NVAR_PASSIVE > 0
void fixup_passive(
    int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]);
void init_passives();
void name_passives();
//#endif

// phys.c
void   primtoflux(double *pr, struct of_state *q, int dir, int magnetic,
      struct of_geom *geom, double *flux);
void   bcon_calc(double *pr, double *ucon, double *ucov, double *bcon);
void   mhd_calc(double *pr, int dir, int magnetic, struct of_state *q,
      double *extra, double *mhd);
void   source(double *ph, struct of_geom *geom, int ii, int jj, double *dU,
      double Dt, double *extra);
double bsq_calc(double *pr, struct of_geom *geom);
void   get_state(double *pr, struct of_geom *geom, struct of_state *q);
void   ucon_calc(double *pr, struct of_geom *geom, double *ucon);
int    mhd_gamma_calc(double *pr, struct of_geom *geom, double *gamma);
void   mhd_vchar(double *pr, struct of_state *q, struct of_geom *geom, int js,
      double *vmax, double *vmin);

// problem.c
void set_problem_params();
void init_prob();
void bound_gas_prob_x1l(int i, int j, int k, grid_prim_type P);
void bound_gas_prob_x1r(int i, int j, int k, grid_prim_type P);
void bound_gas_prob_x2l(int i, int j, int k, grid_prim_type P);
void bound_gas_prob_x2r(int i, int j, int k, grid_prim_type P);
void bound_gas_prob_x3l(int i, int j, int k, grid_prim_type P);
void bound_gas_prob_x4r(int i, int j, int k, grid_prim_type P);

// push_superphotons.c
#if RADIATION
int push_X_K(double X[NDIM], double Kcov[NDIM], double Kcon[NDIM],
    grid_prim_type P, grid_prim_type Prad, double KdotKprev, int type,
    double dtpush);
int push_superphoton(
    struct of_photon *ph, grid_prim_type P, grid_prim_type Prad, double dtpush);
void push_superphotons(grid_prim_type P, grid_prim_type Prad, double dt);
#endif

// rad_utils.c
#if RADIATION
void init_rad(grid_prim_type Prad);
void init_superphoton_resolution();
void update_superphoton_resolution(grid_prim_type Prad, grid_eosvar_type extra);
double linear_interp_log(double x, double *table, double lmin, double dl);
void   list_remove(struct of_photon **ph, struct of_photon **ph_head,
      struct of_photon **ph_prev);
double get_Thetae(double P[NVAR]);
double scatterer_dimensionless_temp(
    int radtype, int interaction, const struct of_microphysics *m);
double scatterer_number_density(
    int radtype, int interaction, const struct of_microphysics *m);
void precompute_microphysics();
void get_fluid_zone(int i, int j, int k, grid_prim_type Prad,
    grid_eosvar_type extra, struct of_microphysics *m, double Ucon[NDIM],
    double Ucov[NDIM], double Bcon[NDIM], double Bcov[NDIM]);
int  is_null(double Kcov[NDIM], double Kcon[NDIM], double K0, double KdotKprev,
     double *KdotK);
void set_Rmunu();
// int is_null(struct of_photon *ph, double *KdotK);
void   Xtoijk(double X[NDIM], int *i, int *j, int *k);
void   copy_photon(struct of_photon *ph, struct of_photon *phc);
void   print_ph_diag(struct of_photon *ph);
int    get_X_K_interp(struct of_photon *ph, double t_interp, grid_prim_type P,
       double X[NDIM], double Kcov[NDIM], double Kcon[NDIM]);
int    to_be_pushed(double t, double dt, struct of_photon *ph);
double get_dtpush(struct of_photon *ph, double dt);
void   swap_ph(struct of_photon **donor, struct of_photon **recipient);
void   get_nuLnu_bin(double X[NDIM], int *thbin, int *phibin);
void   bin_superphoton_direction(const struct of_photon *ph);
double get_min_dt_cool(grid_prim_type P, grid_eosvar_type extra);
void   set_cooling_time(
      grid_double_type tau_cool, grid_prim_type P, grid_eosvar_type extra);
#if RADIATION == RADTYPE_NEUTRINOS
void record_lepton_flux(const struct of_photon *ph);
void check_nu_type(const char *location);
int  get_lepton_sign(const struct of_photon *ph);
#endif // NEUTRINOS
#endif // RADIATION

// radiation.c
#if RADIATION
double Bnu_inv(double nu, const struct of_microphysics *m); // TODO?
double jnu_inv(
    double nu, int type, const struct of_microphysics *m, double theta);
double alpha_inv_scatt(
    double nu, int type, int interaction, const struct of_microphysics *m);
double alpha_inv_abs(
    double nu, int type, const struct of_microphysics *m, double theta);
double get_fluid_nu(double X[NDIM], double Kcov[NDIM], double Ucon[NDIM]);
double get_bk_angle(double X[NDIM], double K[NDIM], double Ucov[NDIM],
    double Bcov[NDIM], double B);
#endif

// random.c
void   init_random(int seed);
double get_rand();
double get_chisq(double nu);
double get_gaussian(double mu, double sigma);
void   get_ran_dir_3d(double *nx, double *ny, double *nz);

// reconstruction.c
void reconstruct(double ptmp[NMAX + 2 * NG][NVAR], int N,
    double p_l[NMAX + 2 * NG][NVAR], double p_r[NMAX + 2 * NG][NVAR]);

// root_finding.c
// counts root_finding iterations
void initialize_root_fcounts();
void print_root_fcounts();
// solves for f(x,params) - ytarget = 0
int find_root(double (*f)(const double, const void *), const void *params,
    const double ytarget, double xguess, const double xmin, const double xmax,
    const double xtol, const double ytol, double *xroot);

// scattering.c
#if RADIATION
int    scatt_temp_too_small(const struct of_microphysics *m);
int    scatter_superphoton(grid_prim_type P, grid_eosvar_type extra,
       struct of_photon *ph, double X[NDIM], double Kcov[NDIM], double Kcon[NDIM],
       int interaction);
void   init_all_hotcross();
double total_cross_lkup(
    double w, int type, int interaction, const struct of_microphysics *m);
#endif

// step.c
void step();

// tetrads.c
#if RADIATION
void coord_to_tetrad(
    double Ecov[NDIM][NDIM], double Kcoord[NDIM], double Ktetrad[NDIM]);
void tetrad_to_coord(
    double Econ[NDIM][NDIM], double Ktetrad[NDIM], double Kcoord[NDIM]);
void make_tetrad(int i, int j, int k, double Ucon[NDIM], double trial[NDIM],
    double Gcov[NDIM][NDIM], double Econ[NDIM][NDIM], double Ecov[NDIM][NDIM]);
void normalize_null(double Gcov[NDIM][NDIM], double K[NDIM]);
void normalize_null_cov(double Gcon[NDIM][NDIM], double K[NDIM]);
#endif

// timing.c
void   time_set(int n, double val);
double time_read(int n);
void   time_init();
void   timer_start(int timerCode);
void   timer_stop(int timerCode);
double get_time_per_step(int timerCode);
void   report_performance();
void   timers_reset();

// tracers.c
#if RADIATION && TRACERS
double get_total_tracer_mass();
int    tracer_max_id(long unsigned int ntracers);
int    count_tracers_local();
void   lagrange_interp_grid_3d(
      double v_interp[], double X[NDIM], double *v, int size);
void lagrange_interp_prim_3d(
    double P_interp[NVAR], double X[NDIM], grid_prim_type P);
void   push_tracers(double X[NDIM], double Kcov[NDIM], double Kcon[NDIM],
      grid_prim_type P, grid_prim_type Ph, double dt);
int    tracer_get_id(struct of_photon *ph);
double tracer_get_mass(struct of_photon *ph);
void   tracer_get_X_u(struct of_photon *ph, double t_interp, double X[NDIM],
      double ucov[NDIM], double ucon[NDIM]);
void   set_tracer(struct of_photon *ph, int id, int nstep, double t,
      double X[NDIM], double mass, grid_prim_type P);
void   make_tracer(struct of_photon **head, int id, int nstep, double t,
      double X[NDIM], double mass, grid_prim_type P);
void   sample_tracers_in_cell(struct of_photon **head, int nstep, int i, int j,
      int k, int *last_id, int ntracers_per_cell, double t, grid_prim_type P);
void   sample_all_tracers(long unsigned int ntracers_tot, int nstep, double t,
      grid_int_type tcrs_in_cell, grid_prim_type P);
void   prune_tracers();
#endif

// util.c
double find_min(const double *array, int size);
double find_max(const double *array, int size);
double find_median(double *array, int size);
double interp_1d(double x, const double xmin, const double xmax, const int imin,
    const int imax, const double *restrict tab_x, const double *restrict tab_y);
int    find_index(double value, const double *array, int size);
void * safe_malloc(size_t size);
void   safe_system(const char *command);
void   safe_fscanf(FILE *stream, const char *format, ...);
int    is_practically_nan(double v);
#if NEED_UNITS
void set_units();
#endif

// utop.c
int Utoprim(double U[NVAR], struct of_geom *geom, double prim[NVAR]);

// xdmf_output.c
void write_xml_file(int dump_id, double t, const char *vnams[NVAR]);


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

int main()
{
  int MAT_DIM = 4;
  int i,j;
  int myArray[16] = {4, 3, 2, 1, 1,10, 3, 4, 5, 3, 2, -4, 4, 8, 7, 9};

  gsl_matrix *m = gsl_matrix_alloc(MAT_DIM,MAT_DIM);

  for (i=0; i<MAT_DIM; i++)
    for (j=0; j<MAT_DIM; j++)
      gsl_matrix_set (m, i, j, (double) myArray[i*MAT_DIM + j]);

  for (i=0; i<MAT_DIM; i++)
    for (j=0; j<MAT_DIM; j++)
      printf("m(%d,%d) = %g\n", i, j,
              gsl_matrix_get (m, i, j));

  double det;
  int signum;
  gsl_permutation *p = gsl_permutation_alloc(m->size1);

  gsl_linalg_LU_decomp(m, p, &signum);
  det = gsl_linalg_LU_det(m, signum);
  gsl_permutation_free(p);
  gsl_matrix_free(m);

  printf("Determinant: %g\n", det);

  return 0;
}

ext: c

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar


text:
/** Spiking Neuron Neural Simulator
* Written by Francis Jeanson 2012
*/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

/* Code Utils */
#define TRUE 1
#define FALSE 0

/* Simulation Defs */
#define TIME_SCALE_FACTOR 1
#define SIMULTIME 1.0   // Simulation Duration
#define STEPSIZE 0.0001 // 0.00001
#define tau 0.4         // 0.04 /* integration constant */
#define TOTAL_STEPS (int)(SIMULTIME/STEPSIZE)

/* Network Model */
#define NET_SIZE 24
#define SHEET_WIDTH 5
#define SHEET_HEIGHT 5

#define INHIB_PERCENT 0.0               // percentage of inhibitory nodes
#define DEFAULT_THRESHOLD 12            // Number of necessary synchronous spikes to reach threshold on target
#define WEIGHT (0.12/DEFAULT_THRESHOLD) // Default connection weight // 0.12 minimum weight for excitation by one afferent

#define CONNECT_RANDOM TRUE
#define CONNECTION_PROB 1.0
#define MIN_DELAY (int)(0.01/STEPSIZE)   // in slots: seconds/resolution -- with tau 0.4 STEPSIZE 0.0001: relative refractory 92ts,  absolute refractory 88ts
#define MAX_DELAY (int)(0.03/STEPSIZE)   // in slots: seconds/resolution

#define CONNECT_TOPOGRAPHIC FALSE
#define CONNECT_EXCIT_PROB 0.05         // probability that an excitatory connection is made
#define CONNECT_INHIB_PROB 0.125        // probability that an inhibitory connection is made
#define DELAY_EXCIT_DECAY_CONST 140     // Should be proportional to SHEET_WIDTH and SHEET_HEIGHT
#define DELAY_INHIB_DECAY_CONST DELAY_EXCIT_DECAY_CONST//(3*DELAY_EXCIT_DECAY_CONST);
#define MIN_TOPO_DELAY (int)(0.0005/STEPSIZE)
#define MAX_TOPO_DELAY (int)(0.01/STEPSIZE)

#define CONNECT_SMALL_WORLD FALSE   // connect network using preferential attachement
#define SMALL_WORLD_CONSTANT 0.02   // a smaller constant means smaller std

#define STDP_WINDOW (int)(0.03/STEPSIZE)
#define STDP_WEIGHT_FACTOR 0.03
#define MAX_WEIGHT 0.08

typedef int boolean;

typedef struct {
    float v;
    float r;
    float I;
    float xpos;
    float ypos;
    float axon[MAX_DELAY];
    int is_excit;
    int dpush; // the circular index to push output voltage on this axon
    float v_hist[STDP_WINDOW];
    int I_duration; //used to determine if input is long enough
} Neuron;

typedef struct {
    float weight;
    boolean are_connected;
    int delay;
} Connection;

const gsl_rng * r; // libgsl random var

void save_nodeconnections(int node_id, Connection Connections[][NET_SIZE], Neuron *Neurons, char *file_name);
double std_int(int *value_list, int length);
double std_float(float *value_list, int length);
int on_sheet_area(int i, int x, int y, int rad);
int get_node_x(int i);
int get_node_y(int i);
void generate_network(Connection Connections[][NET_SIZE], Neuron *Neurons, int mindelay, int maxdelay, float connect_prob, float weight, char *file_name);
int load_network_file(char *file_name, Connection Connections[][NET_SIZE], Neuron *Neurons);

ext: h

max_stars_repo_name: FrancisJMR/SpikeNeuronFN

max_issues_repo_name: FrancisJMR/SpikeNeuronFN

max_forks_repo_name: FrancisJMR/SpikeNeuronFN


text:
/* integration/tests.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

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

#include "tests.h"

/* These are the test functions from table 4.1 of the QUADPACK book */

/* f1(x) = x^alpha * log(1/x) */
/* integ(f1,x,0,1) = 1/(alpha + 1)^2 */

double f1 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(x,alpha) * log(1/x) ;
}

/* f2(x) = 4^-alpha / ((x-pi/4)^2 + 16^-alpha) */
/* integ(f2,x,0,1) = arctan((4-pi)4^(alpha-1)) + arctan(pi 4^(alpha-1)) */

double f2 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(4.0,-alpha) / (pow((x-M_PI/4.0),2.0) + pow(16.0,-alpha)) ;
}

/* f3(x) = cos(2^alpha * sin(x)) */
/* integ(f3,x,0,pi) = pi J_0(2^alpha) */

double f3 (double x, void * params) {
  double alpha = *(double *) params ;
  return cos(pow(2.0,alpha) * sin(x)) ;
}

/* Functions 4, 5 and 6 are duplicates of functions  1, 2 and 3 */
/* ....                                                         */

/* f7(x) = |x - 1/3|^alpha */
/* integ(f7,x,0,1) = ((2/3)^(alpha+1) + (1/3)^(alpha+1))/(alpha + 1) */

double f7 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(fabs(x - (1.0/3.0)),alpha) ;
}

/* f8(x) = |x - pi/4|^alpha */
/* integ(f8,x,0,1) =
   ((1 - pi/4)^(alpha+1) + (pi/4)^(alpha+1))/(alpha + 1) */

double f8 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(fabs(x - (M_PI/4.0)),alpha) ;
}

/* f9(x) = sqrt(1 - x^2) / (x + 1 + 2^-alpha) */
/* integ(f9,x,-1,1) = pi/sqrt((1+2^-alpha)^2-1) */

double f9 (double x, void * params) {
  double alpha = *(double *) params ;
  return 1 / ((x + 1 + pow(2.0,-alpha)) * sqrt(1-x*x)) ;
}

/* f10(x) = sin(x)^(alpha - 1) */
/* integ(f10,x,0,pi/2) = 2^(alpha-2) ((Gamma(alpha/2))^2)/Gamma(alpha) */

double f10 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(sin(x), alpha-1) ;
}

/* f11(x) = log(1/x)^(alpha - 1) */
/* integ(f11,x,0,1) = Gamma(alpha) */

double f11 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(log(1/x), alpha-1) ;
}

/* f12(x) = exp(20*(x-1)) * sin(2^alpha * x) */
/* integ(f12,x,0,1) =
   (20 sin(2^alpha) - 2^alpha cos(2^alpha) + 2^alpha exp(-20))
   /(400 + 4^alpha) */

double f12 (double x, void * params) {
  double alpha = *(double *) params ;
  return exp(20*(x-1)) * sin(pow(2.0,alpha) * x) ;
}

/* f13(x) = cos(2^alpha * x)/sqrt(x(1 - x)) */
/* integ(f13,x,0,1) = pi cos(2^(alpha-1)) J_0(2^(alpha-1))  */

double f13 (double x, void * params) {
  double alpha = *(double *) params ;
  return cos(pow(2.0,alpha)*x)/sqrt(x*(1-x)) ;
}

double f14 (double x, void * params) {
  double alpha = *(double *) params ;
  return exp(-pow(2.0,-alpha)*x)*cos(x)/sqrt(x) ;
}

double f15 (double x, void * params) {
  double alpha = *(double *) params ;
  return x*x * exp(-pow(2.0,-alpha)*x) ;
}

double f16 (double x, void * params) {
  double alpha = *(double *) params ;
  if (x==0 && alpha == 1) return 1 ;  /* make the function continuous in x */
  if (x==0 && alpha > 1) return 0 ;   /* avoid problems with pow(0,1) */
  return pow(x,alpha-1)/pow((1+10*x),2.0) ;
}

double f17 (double x, void * params) {
  double alpha = *(double *) params ;
  return pow(2.0,-alpha)/(((x-1)*(x-1)+pow(4.0,-alpha))*(x-2)) ;
}

/* f454(x) = x^3 log|(x^2-1)(x^2-2)| */
/* integ(f454,x,0,inf) = 61 log(2) + (77/4) log(7) - 27 */

double f454 (double x, void * params) {
  double x2 = x * x;
  double x3 = x * x2;
  params = 0 ;
  return x3 * log(fabs((x2 - 1.0) * (x2 - 2.0))) ;
}

/* f455(x) = log(x)/(1+100*x^2) */
/* integ(f455,x,0,inf) = -log(10)/20 */

double f455 (double x, void * params) {
  params = 0 ;
  return log(x) / (1.0 + 100.0 * x * x) ;
}

/* f456(x) = log(x) */
/* integ(f456*sin(10 pi x),x,0,1) = -(gamma + log(10pi) - Ci(10pi))/(10pi) */

double f456 (double x, void * params) {
  params = 0 ;
  if (x == 0.0)
    {
      return 0;
    }
  return log(x) ;
}

/* f457(x) = 1/sqrt(x) */
/* integ(f457*cos(pi x / 2),x,0,+inf) = 1 */

double f457 (double x, void * params) {
  params = 0 ;
  if (x == 0.0)
    {
      return 0;
    }
  return 1/sqrt(x) ;
}

/* f458(x) = 1/(1 + log(x)^2)^2 */
/* integ(log(x) f458(x),x,0,1) = (Ci(1) sin(1) + (pi/2 - Si(1)) cos(1))/pi
                               = -0.1892752 */

double f458 (double x, void * params) {
  params = 0 ;

  if (x == 0.0)
    {
      return 0;
    }
  else
    {
      double u = log(x);
      double v = 1 + u * u;

      return 1.0 / (v * v) ;
    }
}

/* f459(x) = 1/(5 x^3 + 6) */
/* integ(f459/(x-0),x,-1,5) = log(125/631)/18 */

double f459 (double x, void * params) {
  params = 0 ;
  return 1.0 / (5.0 * x * x * x + 6.0) ;
}

/* myfn1(x) = exp(-x - x^2) */
/* integ(myfn1,x,-inf,inf) = sqrt(pi) exp(-1/4) */

double myfn1 (double x, void * params) {
  params = 0;
  return exp(-x - x*x) ;
}

/* myfn2(x) = exp(alpha*x) */
/* integ(myfn2,x,-inf,b) = exp(alpha*b)/alpha */

double myfn2 (double x, void * params) {
  double alpha = *(double *) params ;
  return exp(alpha*x) ;
}


text:
static char help[] =
"ODE system solver example using TS.  Solves a 2D system\n"
"    dy/dt = G(t,y)\n"
"with y(t0) = y0 to compute y(tf).  Sets problem type to NONLINEAR and\n"
"TS type to Runge-Kutta.  No Jacobian is supplied; compare odejac.c.\n"
"Exact solution is known.\n\n";

#include <petsc.h>

extern PetscErrorCode ExactSolution(PetscReal, Vec);
extern PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void*);

//STARTMAIN
int main(int argc,char **argv) {
  PetscErrorCode ierr;
  PetscInt   steps;
  PetscReal  t0 = 0.0, tf = 20.0, dt = 0.1, err;
  Vec        y, yexact;
  TS         ts;

  ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

  ierr = VecCreate(PETSC_COMM_WORLD,&y); CHKERRQ(ierr);
  ierr = VecSetSizes(y,PETSC_DECIDE,2); CHKERRQ(ierr);
  ierr = VecSetFromOptions(y); CHKERRQ(ierr);
  ierr = VecDuplicate(y,&yexact); CHKERRQ(ierr);

  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,NULL); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);

  // set time axis
  ierr = TSSetTime(ts,t0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,tf); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);

  // set initial values and solve
  ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
  ierr = ExactSolution(t0,y); CHKERRQ(ierr);
  ierr = TSSolve(ts,y); CHKERRQ(ierr);

  // compute error and report
  ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
  ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
  ierr = ExactSolution(tf,yexact); CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,yexact); CHKERRQ(ierr);    // y <- y - yexact
  ierr = VecNorm(y,NORM_INFINITY,&err); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
              "error at tf = %.3f with %d steps:  |y-y_exact|_inf = %g\n",
              tf,steps,err); CHKERRQ(ierr);

  VecDestroy(&y);  VecDestroy(&yexact);  TSDestroy(&ts);
  return PetscFinalize();
}
//ENDMAIN

//STARTCALLBACKS
PetscErrorCode ExactSolution(PetscReal t, Vec y) {
    PetscReal *ay;
    VecGetArray(y,&ay);
    ay[0] = t - PetscSinReal(t);
    ay[1] = 1.0 - PetscCosReal(t);
    VecRestoreArray(y,&ay);
    return 0;
}

PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec y, Vec g,
                               void *ptr) {
    const PetscReal *ay;
    PetscReal       *ag;
    VecGetArrayRead(y,&ay);
    VecGetArray(g,&ag);
    ag[0] = ay[1];            // = g_1(t,y)
    ag[1] = - ay[0] + t;      // = g_2(t,y)
    VecRestoreArrayRead(y,&ay);
    VecRestoreArray(g,&ag);
    return 0;
}
//ENDCALLBACKS


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

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_ieee_utils.h>

int test_nist (void);

/* Test program for mean.c.  JimDavies 7.96 */

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "test_float_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "test_float_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "test_float_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "test_int_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

int
main (void)
{
  size_t s1, s2;

  gsl_ieee_env_setup ();

  for (s1 = 1; s1 < 4 ; s1++)
    {
      s2 = (s1 < 3) ? 1 : (s1 - 1) ;

      test_func (s1,s2);
      test_float_func (s1,s2);
      test_long_double_func (s1,s2);

      test_ulong_func (s1,s2);
      test_long_func (s1,s2);
      test_uint_func (s1,s2);
      test_int_func (s1,s2);
      test_ushort_func (s1,s2);
      test_short_func (s1,s2);
      test_uchar_func (s1,s2);
      test_char_func (s1,s2);
    }

  test_nist();

  exit (gsl_test_summary ());
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include "BoundingVolume.h"
#include "Library.h"
#include "TileContentLoadResult.h"
#include "TileContentLoader.h"
#include "TileID.h"
#include "TileRefine.h"

#include <CesiumGltf/GltfReader.h>

#include <glm/mat4x4.hpp>
#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>

namespace Cesium3DTilesSelection {

class Tileset;

/**
 * @brief Creates {@link TileContentLoadResult} from glTF data.
 */
class CESIUM3DTILESSELECTION_API GltfContent final : public TileContentLoader {
public:
  /**
   * @copydoc TileContentLoader::load
   *
   * The result will only contain the `model`. Other fields will be
   * empty or have default values.
   */
  CesiumAsync::Future<std::unique_ptr<TileContentLoadResult>>
  load(const TileContentLoadInput& input) override;

  /**
   * @brief Create a {@link TileContentLoadResult} from the given data.
   *
   * (Only public to be called from `Batched3DModelContent`)
   *
   * @param pLogger Only used for logging
   * @param url The URL, only used for logging
   * @param data The actual glTF data
   * @return The {@link TileContentLoadResult}
   */
  static std::unique_ptr<TileContentLoadResult> load(
      const std::shared_ptr<spdlog::logger>& pLogger,
      const std::string& url,
      const gsl::span<const std::byte>& data);

  /**
   * @brief Creates texture coordinates for mapping {@link RasterOverlay} tiles
   * to {@link Tileset} tiles.
   *
   * Generates new texture coordinates for the `gltf` using the given
   * `projection`. The first new texture coordinate (`u` or `s`) will be 0.0 at
   * the `minimumX` of the given `rectangle` and 1.0 at the `maximumX`. The
   * second texture coordinate (`v` or `t`) will be 0.0 at the `minimumY` of
   * the given `rectangle` and 1.0 at the `maximumY`.
   *
   * Coordinate values for vertices in between these extremes are determined by
   * projecting the vertex position with the `projection` and then computing the
   * fractional distance of that projected position between the minimum and
   * maximum.
   *
   * Projected positions that fall outside the `rectangle` will be clamped to
   * the edges, so the coordinate values will never be less then 0.0 or greater
   * than 1.0.
   *
   * These texture coordinates are stored in the provided glTF, and a new
   * primitive attribute named `_CESIUMOVERLAY_n`, where `n` is the
   * `textureCoordinateID` passed to this function, is added to each primitive.
   *
   * @param gltf The glTF model.
   * @param transform The transformation of this glTF to ECEF coordinates.
   * @param textureCoordinateID The texture coordinate ID.
   * @param projection The projection. There is a linear relationship between
   * the coordinates of this projection and the generated texture coordinates.
   * @param rectangle The rectangle that all projected vertex positions are
   * expected to lie within.
   * @return The bounding region.
   */
  static CesiumGeospatial::BoundingRegion createRasterOverlayTextureCoordinates(
      CesiumGltf::Model& gltf,
      const glm::dmat4& transform,
      int32_t textureCoordinateID,
      const CesiumGeospatial::Projection& projection,
      const CesiumGeometry::Rectangle& rectangle);

private:
  static CesiumGltf::GltfReader _gltfReader;
};

} // namespace Cesium3DTilesSelection

max_stars_repo_name: JiangMuWen/cesium-native

max_issues_repo_name: JiangMuWen/cesium-native

max_forks_repo_name: JiangMuWen/cesium-native


text:
/*
 * Copyright (c) 2010-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include <lapacke.h>
#include "dplasma.h"
#include "dplasmatypes.h"

#include "map2.h"

static int
dplasma_zlacpy_sum_operator( parsec_execution_stream_t *es,
                             const parsec_tiled_matrix_dc_t *descA,
                             const parsec_tiled_matrix_dc_t *descB,
                             const void *_A, void *_B,
                             PLASMA_enum uplo, int m, int n,
                             void *args )
{
    int tempmm, tempnn, ldam, ldbm;
    const parsec_complex64_t *A = (const parsec_complex64_t*)_A;
    parsec_complex64_t       *B = (parsec_complex64_t*)_B;
    (void)es;
    (void)args;

    tempmm = ((m)==((descA->mt)-1)) ? ((descA->m)-(m*(descA->mb))) : (descA->mb);
    tempnn = ((n)==((descA->nt)-1)) ? ((descA->n)-(n*(descA->nb))) : (descA->nb);
    ldam = BLKLDD( descA, m );
    ldbm = BLKLDD( descB, m );

    parsec_complex64_t* unit_v = (parsec_complex64_t*)malloc(tempnn*sizeof(parsec_complex64_t));
    parsec_complex64_t* int_v = (parsec_complex64_t*)malloc(tempnn*sizeof(parsec_complex64_t));
    for(int i = 0; i < tempnn; i++ ) {
        unit_v[i] = (parsec_complex64_t)1.0;
        int_v[i] = (parsec_complex64_t)(1.0*(i+1));
    }
    parsec_complex64_t* unit_vt = (parsec_complex64_t*)malloc((tempmm+1)*sizeof(parsec_complex64_t));
    parsec_complex64_t* int_vt = (parsec_complex64_t*)malloc((tempmm+1)*sizeof(parsec_complex64_t));
    for(int i = 0; i < tempmm; i++ ) {
        unit_vt[i] = (parsec_complex64_t)1.0;
        int_vt[i] = (parsec_complex64_t)(1.0*(i+1));
    }

    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR, lapack_const( PlasmaUpperLower ), tempmm, tempnn, A, ldam, B, ldbm);

    CORE_zgemv(PlasmaNoTrans, tempmm, tempnn,
               (parsec_complex64_t)1.0, B, ldbm,
               unit_v, 1,
               (parsec_complex64_t)0.0, B+ldbm*(tempnn), 1);

    CORE_zgemv(PlasmaNoTrans, tempmm, tempnn,
               (parsec_complex64_t)1.0, B, ldbm,
               int_v, 1,
               (parsec_complex64_t)0.0, B+ldbm*(tempnn+1), 1);

    CORE_zgemv(PlasmaTrans, tempmm, tempnn+1,
               (parsec_complex64_t)1.0, B, ldbm,
               unit_vt, 1,
               (parsec_complex64_t)0.0, B+(tempmm), ldbm);

    CORE_zgemv(PlasmaTrans, tempmm, tempnn+1,
               (parsec_complex64_t)1.0, B, ldbm,
               int_vt, 1,
               (parsec_complex64_t)0.0, B+(tempmm+1), ldbm);

    free(unit_v);
    free(unit_vt);
    free(int_v);
    free(int_vt);

    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy_New - Generates an object that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2_New() for further information.
 *
 * WARNING: The computations are not done by this call.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval NULL if incorrect parameters are given.
 *          \retval The parsec taskpool describing the operation that can be
 *          enqueued in the runtime. It, then, needs to be
 *          destroy with the corresponding destruct function.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy_New
 * @sa dplasma_dlacpy_New
 * @sa dplasma_slacpy_New
 *
 ******************************************************************************/
parsec_taskpool_t*
dplasma_zlacpy_sum_New( PLASMA_enum uplo,
                        const parsec_tiled_matrix_dc_t *A,
                        parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t* tp;

    tp = dplasma_map2_New(uplo, PlasmaNoTrans, A, B,
                          dplasma_zlacpy_sum_operator, NULL );

    return tp;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 *  dplasma_zlacpy_Destruct - Free the data structure associated to an object
 *  created with dplasma_zlacpy_New().
 *
 *******************************************************************************
 *
 * @param[in,out] o
 *          On entry, the object to destroy.
 *          On exit, the object cannot be used anymore.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy
 *
 ******************************************************************************/
void
dplasma_zlacpy_sum_Destruct( parsec_taskpool_t *tp )
{
    parsec_taskpool_free(tp);
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlacpy - Generates an object that performs a copy of the matrix A
 * into the matrix B.
 *
 * See dplasma_map2() for further information.
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The PaRSEC context of the application that will run the operation.
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is copied:
 *          = PlasmaUpperLower: All matrix is referenced.
 *          = PlasmaUpper:      Only upper part is refrenced.
 *          = PlasmaLower:      Only lower part is referenced.
 *
 * @param[in] A
 *          Descriptor of the distributed original matrix A. Any tiled matrix
 *          descriptor can be used. However, if the data is stored in column
 *          major, the tile distribution must match the one of the matrix B.
 *
 * @param[in,out] B
 *          Descriptor of the distributed destination matrix B. Any tiled matrix
 *          descriptor can be used, with no specific storage.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlacpy_New
 * @sa dplasma_zlacpy_Destruct
 * @sa dplasma_clacpy
 * @sa dplasma_dlacpy
 * @sa dplasma_slacpy
 *
 ******************************************************************************/
int
dplasma_zlacpy_sum( parsec_context_t *parsec,
                PLASMA_enum uplo,
                const parsec_tiled_matrix_dc_t *A,
                parsec_tiled_matrix_dc_t *B)
{
    parsec_taskpool_t *parsec_zlacpy_sum = NULL;

    if ((uplo != PlasmaUpperLower) &&
        (uplo != PlasmaUpper)      &&
        (uplo != PlasmaLower))
    {
        dplasma_error("dplasma_zlacpy_sum", "illegal value of uplo");
        return -2;
    }

    parsec_zlacpy_sum = dplasma_zlacpy_sum_New(uplo, A, B);

    if ( parsec_zlacpy_sum != NULL )
    {
        parsec_context_add_taskpool(parsec, parsec_zlacpy_sum);
        dplasma_wait_until_completion(parsec);
        dplasma_zlacpy_sum_Destruct( parsec_zlacpy_sum );
    }
    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
#ifndef areaproblemdata_h
#define areaproblemdata_h

#include <ceed.h>
#include <petsc.h>
#include "../include/structs.h"
#include "../qfunctions/area/areacube.h"
#include "../qfunctions/area/areasphere.h"

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

// Problem options
typedef enum {
  CUBE = 0, SPHERE = 1
} ProblemType;

static BPData problem_options[6] = {
  [CUBE] = {
    .num_comp_x = 3,
    .num_comp_u = 1,
    .topo_dim = 2,
    .q_data_size = 1,
    .q_extra = 1,
    .setup_geo = SetupMassGeoCube,
    .apply = Mass,
    .setup_geo_loc = SetupMassGeoCube_loc,
    .apply_loc = Mass_loc,
    .in_mode = CEED_EVAL_INTERP,
    .out_mode = CEED_EVAL_INTERP,
    .q_mode = CEED_GAUSS,
    .enforce_bc = PETSC_FALSE,
  },
  [SPHERE] = {
    .num_comp_x = 3,
    .num_comp_u = 1,
    .topo_dim = 2,
    .q_data_size = 1,
    .q_extra = 1,
    .setup_geo = SetupMassGeoSphere,
    .apply = Mass,
    .setup_geo_loc = SetupMassGeoSphere_loc,
    .apply_loc = Mass_loc,
    .in_mode = CEED_EVAL_INTERP,
    .out_mode = CEED_EVAL_INTERP,
    .q_mode = CEED_GAUSS,
    .enforce_bc = PETSC_FALSE,
  }
};

#endif // areaproblemdata_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>

#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>

#include "common-driver-utils.h"
#include "stokes_block_scaling.h"
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
#include <Underworld/Underworld.h>
#include "Solvers/SLE/SLE.h" /* to give the AugLagStokes_SLE type */
#include "Solvers/KSPSolvers/KSPSolvers.h" /* for __KSP_COMMON */

#include "BSSCR.h"
#include "writeMatVec.h"

/* creates and builds the "checker-board" null-space vectors for pressure
   and sets the t and v "Vec pointers" on the bsscr struct to point to them */
#undef __FUNCT__
#define __FUNCT__ "KSPBuildPressure_CB_Nullspace_BSSCR"
PetscErrorCode KSPBuildPressure_CB_Nullspace_BSSCR(KSP ksp)
{
    KSP_BSSCR        *bsscr = (KSP_BSSCR *)ksp->data;
    FeEquationNumber *eq_num = bsscr->solver->st_sle->pSolnVec->eqNum;
    FeMesh           *feMesh = bsscr->solver->st_sle->pSolnVec->feVariable->feMesh; /* is the pressure mesh */
    unsigned          ijk[3];
    Vec               t, v;
    int numLocalNodes, globalNodeNumber, i, j, eq;
    MatStokesBlockScaling BA = bsscr->BA;
    PetscErrorCode ierr;
    Mat            Amat,Pmat, G;
    MatStructure   pflag;

    PetscFunctionBegin;
    /* get G matrix from Amat matrix operator on ksp */
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
    /* now create Vecs t and v to match size of G: i.e. pressure */ /* NOTE: not using "h" vector from ksp->vec_rhs because this part of the block vector doesn't always exist */
    MatGetVecs( G, &t, PETSC_NULL );/* t and v are destroyed in KSPDestroy_BSSCR */
    MatGetVecs( G, &v, PETSC_NULL );/* t and v such that can do G*t */

    numLocalNodes = Mesh_GetLocalSize( feMesh,  MT_VERTEX); /* number of nodes on current proc not counting any shadow nodes */
    for(j=0;j<numLocalNodes;j++){
	i = globalNodeNumber = Mesh_DomainToGlobal( feMesh, MT_VERTEX, j);
	RegularMeshUtils_Element_1DTo3D(feMesh, i, ijk);
	eq = eq_num->mapNodeDof2Eq[j][0];/* get global equation number -- 2nd arg is always 0 because pressure has only one dof */
	if(eq != -1){
	    if( (ijk[0]+ijk[1]+ijk[2])%2 ==0 ){
		VecSetValue(t,eq,1.0,INSERT_VALUES);
	    }
	    else{
		VecSetValue(v,eq,1.0,INSERT_VALUES);	    }}}

    VecAssemblyBegin( t );
    VecAssemblyEnd( t );
    VecAssemblyBegin( v );
    VecAssemblyEnd( v );

    /* Scaling the null vectors here because it easier at the moment *//* maybe should do this in the original scaling function */
    if( BA->scaling_exists == PETSC_TRUE ){
	Vec R2;
	/* Get the scalings out the block mat data */
	VecNestGetSubVec( BA->Rz, 1, &R2 );
	VecPointwiseDivide( t, t, R2); /* x <- x * 1/R2 */
	VecPointwiseDivide( v, v, R2);
    }

    // bsscr_writeVec( t, "t", "Writing t vector");
    // bsscr_writeVec( v, "v", "Writing v vector");

    bsscr->t=t;
    bsscr->v=v;

    PetscFunctionReturn(0);
}
/* creates and builds the "constant" null-space vector for pressure
   and sets the t "Vec pointer" on the bsscr struct to point to it */
#undef __FUNCT__
#define __FUNCT__ "KSPBuildPressure_Const_Nullspace_BSSCR"
PetscErrorCode KSPBuildPressure_Const_Nullspace_BSSCR(KSP ksp)
{
    KSP_BSSCR        *bsscr = (KSP_BSSCR *)ksp->data;
    FeEquationNumber *eq_num = bsscr->solver->st_sle->pSolnVec->eqNum;
    FeMesh           *feMesh = bsscr->solver->st_sle->pSolnVec->feVariable->feMesh; /* is the pressure mesh */
    unsigned          ijk[3];
    Vec               t;
    int numLocalNodes, globalNodeNumber, i, eq;
    MatStokesBlockScaling BA = bsscr->BA;
    PetscErrorCode ierr;
    Mat            Amat,Pmat, G;
    MatStructure   pflag;

    PetscFunctionBegin;
    /* get G matrix from Amat matrix operator on ksp */
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
    /* now create Vec t to match size of G: i.e. pressure */ /* NOTE: not using "h" vector from ksp->vec_rhs because this part of the block vector doesn't always exist */
    MatGetVecs( G, &t, PETSC_NULL );/* t is destroyed in KSPDestroy_BSSCR */

    VecSet(t, 1.0);

    VecAssemblyBegin( t );
    VecAssemblyEnd( t );

    /* Scaling the null vectors here because it easier at the moment */
    if( BA->scaling_exists == PETSC_TRUE ){
	Vec R2;
	/* Get the scalings out the block mat data */
	VecNestGetSubVec( BA->Rz, 1, &R2 );
	VecPointwiseDivide( t, t, R2); /* x <- x * 1/R2 */
    }
    // bsscr_writeVec( t, "t", "Writing t vector");
    bsscr->t=t;

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "KSPRemovePressureNullspace_BSSCR"
PetscErrorCode KSPRemovePressureNullspace_BSSCR(KSP ksp, Vec h_hat)
{
    KSP_BSSCR            *bsscr = (KSP_BSSCR *)ksp->data;
    MatStokesBlockScaling BA    = bsscr->BA;
    PetscErrorCode ierr;
    PetscScalar norm, a, a1, a2, hnorm, pnorm, gnorm;
    Vec t2, t,v;
    Mat            Amat,Pmat, D;
    MatStructure   pflag;
    double nstol   = bsscr->nstol; /* set in KSPCreate_BSSCR */

    PetscFunctionBegin;
    t=bsscr->t;
    v=bsscr->v;
    /* get G matrix from Amat matrix operator on ksp */
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    //MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
    MatNestGetSubMat( Amat, 1,0, &D );/* D should always exist */
    /* now create Vec t2 to match left hand size of G: i.e. velocity */
    MatGetVecs( D, &t2, PETSC_NULL );
    MatNorm(D,NORM_INFINITY,&gnorm); /* seems like not a bad estimate of the largest eigenvalue for this matrix */
    if(t){/* assumes that v and t are initially set to PETSC_NULL (in KSPCreate_BSSCR ) */
	MatMultTranspose( D, t, t2);
	VecNorm(t2, NORM_2, &norm);
	VecNorm(t, NORM_2, &a);
	norm=norm/a;/* so we are using unit null vector */
	if(norm < nstol*gnorm){/* then t in NS of G */
	    VecDot(t,h_hat, &a1);
	    VecDot(t,t, &a2);
	    a=-a1/a2;
	    VecAXPY(h_hat, a, t);
	    VecDot(t,h_hat, &a1);
	    PetscPrintf( PETSC_COMM_WORLD, "\n\t* Found t NS norm= %g : Dot = %g\n\n", norm, a1);
	}
    }
    if(v){
	MatMultTranspose( D, v, t2);
	VecNorm(t2, NORM_2, &norm);
	VecNorm(v, NORM_2, &a);
	norm=norm/a;/* so we are using unit null vector */
	if(norm < nstol*gnorm){/* then v in NS of G */
	    VecDot(v,h_hat, &a1);
	    VecDot(v,v, &a2);
	    a=-a1/a2;
	    VecAXPY(h_hat, a, v);
	    VecDot(v,h_hat, &a1);
	    PetscPrintf( PETSC_COMM_WORLD, "\n\t* Found v NS norm= %g : Dot = %g\n\n", norm, a1);
	}
    }
    // bsscr_writeVec( h_hat, "h", "Writing h_hat Vector in Solver");
    Stg_VecDestroy(&t2);

    PetscFunctionReturn(0);
}


text:
/*
   Compile as

   gcc compGraph2graph6.c -O2 -lgsl -lgslcblas -o compGraph2graph6
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>

/* Order of input graphs - note this code will only work for orders smaller than 62 */
#define N 19

/* Maximum order we expect for the the comparability graph */
#define MAX_ORDER 75000

/* Forbidden eigenvalue of the star complement */
#define SC_EIG 2

/* Specified precision used while doing floating point operations we treat EPS as zero */
#define EPS 0.000001

#define EXPECTED_CLIQUE 57

/* graph6 related things */

#define SIZELEN(n) ((n)<=SMALLN?1:((n)<=SMALLISHN?4:8))
#define G6LEN(n) (SIZELEN(n) \
         + ((size_t)(n)/12)*((size_t)(n)-1) + (((size_t)(n)%12)*((size_t)(n)-1)+11)/12)

#define SMALLN 62
#define BIAS6 63
#define TOPBIT6 32
#define SMALLISHN 258047
#define MAXBYTE 126
#define C6MASK 63

static gsl_matrix *mat, *mat_inv;
static gsl_permutation *perm;

static gsl_vector *vecs[1<<N]; //vecs and vecs_prod are arrays of vectors with size 2^N
static gsl_vector *vecs_prod[1<<N];

static gsl_vector *vec_j;

/* storing the current graph6 string of our graph */
static char line[G6LEN(N)+2];

static char gcode[G6LEN(MAX_ORDER)+3]; //gcode is string with enough room to store graph of order MAX_ORDER

/* number of currently processed graphs. */
static unsigned nproc = 0;
static unsigned skipped = 0;

static FILE *outFile;

static void stringtomat(char *s) {

	char *p;
	int i,j,k,x = 0;

    /* Clear the adjacency matrix */
    gsl_matrix_set_zero(mat);

    p = s + 1;

    k = 1;

    gsl_matrix_set(mat, 0, 0, SC_EIG);

    for (j = 1; j < N; ++j) {
        gsl_matrix_set(mat, j, j, SC_EIG);
        for (i = 0; i < j; ++i) {
            if (--k == 0) {
        		k = 6;
		        x = *(p++) - BIAS6;
            }

            if (x & TOPBIT6) {
                gsl_matrix_set(mat, i, j, -1);
                gsl_matrix_set(mat, j, i, -1);
            }
            x <<= 1;
        }
    }

    int signum;
    gsl_linalg_LU_decomp(mat, perm , &signum);
    gsl_linalg_LU_invert(mat, perm, mat_inv);

}

/* Some graph6 string thingie */
static void encodegraphsize(const int n, char **pp) {
    char *p;

    p = *pp;
    if (n <= SMALLN)
        *p++ = BIAS6 + n;
    else {
        *p++ = MAXBYTE;
        *p++ = BIAS6 + (n >> 12);
        *p++ = BIAS6 + ((n >> 6) & C6MASK);
        *p++ = BIAS6 + (n & C6MASK);
    }
    *pp = p;
}

/* Given that mat_inv is the M = SC_EIG*I - A we
   compute all binary vectors x such that x M x^t == SC_EIG
   and x M j^t == -1.

   Finally for any pair x,y of such vectors we add an edge to our
   graph if x M y^t is either 0 or -1.
*/

/* This is a global thingie since declaring it localy as
 * an array of size 1<<N breaks the stack limit */
gsl_vector *verts[1<<N];

static void constructGraph(void) {

    double res; //allocate memory for double-sized result
    unsigned i, j; //declare variables used for looping

    /* After the first iteration this holds the number of
       vertices of the obtained graph. The respective vertices
       are stored in vecs_prod[0]...vecs_prod[cache_size-1].
    */
    unsigned cache_size = 0; //initialize cache_size to 0

    for (i = 0; i < 1<<N; i++) { //loop through all 2^N vectors in vecs, looking to fill in the vertices
        gsl_blas_dsymv(CblasUpper, 1, mat_inv, vecs[i], 0, vecs_prod[cache_size]); //set vecs_prod[cache_size] to mat_inv*vecs[i], assuming mat_inv is symmetrical and using upper triangle ((rI-A_H)^-1*u), this is the bilinear form <u,_>

        gsl_blas_ddot(vecs_prod[cache_size], vec_j, &res); //calculate <u,j> and store it in space allocated for result

        if (fabs(res+1) < EPS) { //if <u,j> = -1 (within EPS precision)
            gsl_blas_ddot(vecs_prod[cache_size], vecs[i], &res); //calculate <u,u> and store it in space allocated for result

            if (fabs(res-SC_EIG) < EPS) { //if <u,u> = r (within EPS precision)
                verts[cache_size] = vecs[i]; //it's a vertex, store vecs[i] in verts
                cache_size+=1; //search for the next vector
            }
        }
    }
    //we now have matched arrays vert[] and vec_prod[], with valid vertices in the former and <u,_> in the latter
    if (cache_size < EXPECTED_CLIQUE) { //if there aren't enough elements to form target clique, no point in continuing, skip to next star complement candidate
        skipped++;
        return;
    }

    assert(cache_size < MAX_ORDER); //check if order is less than maximum allowed order
	//prepare to output compatibility graph
    char *p = gcode; //point *p at start of gcode
    encodegraphsize(cache_size, &p);

    int k = 6, x = 0; //initialize bit written counter and bit to be written

    for (i = 1; i < cache_size; i++) { //loop through all vertices u_i
        for (j = 0; j < i; j++) { //loop through <u_i,u_1> up to <u_i,u_i>, since <u,v> = <v,u> we don't have to check the rest as they will be checked in a later iteration
            x <<= 1; //ready next bit
            gsl_blas_ddot(vecs_prod[i], verts[j], &res); //calculate <u_i,u_j> and store it in space allocated for result

            /* We have an edge */
            if (fabs(res) <= EPS || fabs(res+1) <= EPS) { //if <u_i,u_j> = 0 or -1 (it's an edge)
                x |= 1; //set last bit of x to 1
            }
            if (--k == 0) { //decrease bit written counter
                *p++ = BIAS6 + x; //if we wrote 6 bits, add padding for ASCII character and write byte to where *p is pointing and advance pointer *p forward by 1 character
                k = 6; //reset counters
                x = 0;
            }

        }
    }
	//all full bytes written to gcode, prepare to write to file
    if (k != 6) { //if last byte is incomplete (still waiting to be written)
        *p++ = BIAS6 + (x << k); //push last bits to front, add padding, write to *p, advance *p
    }
    *p++ = '\n'; //stitch on newline, advance *p
    *p = '\0'; //stitch on null (to end string)

    fputs(gcode, outFile); //write to output file
}

static void init_vectors(void) {

    unsigned i,j;

    vec_j = gsl_vector_alloc(N); //allocate memory for vec_j with size N

    assert(vec_j); //check if vec_j was created successfully

    gsl_vector_set_all(vec_j, 1); //initialize vec_j as all 1's

    for (i = 0; i < 1<<N; i++) { //loop through all 2^N vectors in vecs[] and vecs_prod[]
        vecs[i] = gsl_vector_calloc(N); //initialize each vecs[i] to all 0
        vecs_prod[i] = gsl_vector_alloc(N); //allocate memory for vectors of size N

        assert(vecs[i] && vecs_prod[i]); //check if vectors were initialized successfully

        /* We fill the i'th vector of vecs */
        for (j = 0; j < N; j++) { //loop through length of vecs[i]
            if ( i & (1<<j) ) { //set vecs[i] to binary representation of i
                gsl_vector_set(vecs[i], j, 1);
            }
        }
    }
}

int main(int argc, char **argv) { //start here

    static FILE *infile; //declare pointer to input file

    assert (argc > 1); //check if arguments were provided (you can provide multiple, only the first will be used though)

    infile = fopen(argv[1], "r"); //open first argument in read mode and point *infile at it
    mat = gsl_matrix_alloc(N, N); //allocate memory for mat as NxN matrix (not initialized yet, contents are garbage)
    mat_inv = gsl_matrix_alloc(N, N); //allocate memory for mat_inv as NxN matrix
    perm = gsl_permutation_calloc(N); //initialize perm as identity permutation of length N

    char buf[512]; //declare string of 511+1 characters

    snprintf(buf, sizeof(buf), "%s.out", argv[1]); //fill buf with name of input file, append .out to the name
    outFile = fopen(buf, "w"); //create file with name contained in buf in write mode, point *outFile at it

    assert (outFile && infile && mat && mat_inv && perm); //check if everything exists

    init_vectors();

    while (1) {
        if (!fgets(line, sizeof(line), infile))
            break;
        stringtomat(line);
        constructGraph();
        nproc++;
	}

    printf("Successfuly processed: %u graphs. Skipped: %u\n" , nproc, skipped);

    return 0;
}

max_stars_repo_name: ThatEvilChickenNextDoor/GMSwitchedGraph

max_issues_repo_name: ThatEvilChickenNextDoor/GMSwitchedGraph

max_forks_repo_name: ThatEvilChickenNextDoor/GMSwitchedGraph