text:
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
   @file pgp.c
   @brief Interface to pgplot

   Some usage of pgplot.Functionality goes like this: With filling the
   struct pgp_gdsc the user defines a layout of a plot, which consists
   of a number of viewgraphs with the same x-axis range. A default
   struct of that type is delivered by the routine pgp_default. A
   pgplot device is activated (and must be activated) by calling the
   routine pgp_opendev(). Then, a box is drawn by calling the routine
   pgp_openbox. With a call of that function an appropriate pgplot
   viewport is defined to enable the user to plot his viewgraph in the
   box by calling functions pgp_marker(), pgp_lines(), pgp_bars(), and
   pgp_errby(). With pgp_legend(), a string can be plotted to a legend
   line. The routine pgp_end() will terminate the plotting process
   properly.

   An example can be accessed with testpgp.c.

   $Source: /Volumes/DATA_J_II/data/CVS/tirific/src/pgp.c,v $
   $Date: 2011/05/25 22:25:26 $
   $Revision: 1.14 $
   $Author: jozsa $
   $Log: pgp.c,v $
   Revision 1.14  2011/05/25 22:25:26  jozsa
   Left work

   Revision 1.13  2011/05/10 00:30:16  jozsa
   Left work

   Revision 1.12  2009/05/26 07:56:41  jozsa
   Left work

   Revision 1.11  2007/08/22 15:58:43  gjozsa
   Left work

   Revision 1.10  2006/11/22 14:16:21  gjozsa
   Bugfix concerning RASH and horizontal/vertical lines

   Revision 1.9  2006/11/03 10:49:02  gjozsa
   Introduced logarithmic scaling, hms dms

   Revision 1.8  2006/10/05 12:39:58  gjozsa
   nasty bugfix

   Revision 1.7  2006/04/03 11:47:57  gjozsa
   Left work

   Revision 1.6  2006/02/08 16:13:17  gjozsa
   Extended graphics routines

   Revision 1.5  2005/10/12 14:51:25  gjozsa
   First point gets thick in the polar plot

   Revision 1.4  2005/10/12 09:53:09  gjozsa
   Included polar plot facility

   Revision 1.3  2005/06/24 12:00:29  gjozsa
   Left work

   Revision 1.2  2005/05/24 10:42:24  gjozsa
   Left work

   Revision 1.1  2005/05/17 12:58:15  gjozsa
   Added to cvs control

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* EXTERNAL INCLUDES */
/* ------------------------------------------------------------ */
#include <stdio.h>
#include <cpgplot.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include <stdlib.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* INTERNAL INCLUDES */
/* ------------------------------------------------------------ */
#include <pgp.h>
#include <maths.h>

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
   @def _MEMORY_HERE_ON
   @brief Controls the use of the memory_here module

   If you don't want to use the memory_here facility comment this
   define, otherways it will be included.

*/
/* ------------------------------------------------------------ */
/* #define _MEMORY_HERE_ON */
/* #include <memory_here.h> */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* PRIVATE SYMBOLIC CONSTANTS */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
   @define DEGTORAD
   @brief Conversion factor from deg to rad
*/
/* ------------------------------------------------------------ */
#define DEGTORAD 0.0174532925199433

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
   @define LOGFAULT
   @brief If a logarithm is asked for for a negative number or 0 this should be redurned
*/
/* ------------------------------------------------------------ */
#define LOGFAULT -40

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* PRIVATE MACROS */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* (PRIVATE) GLOBAL VARIABLES */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* PRIVATE TYPEDEFS */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* PRIVATE STRUCTS */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* PRIVATE FUNCTION DECLARATIONS */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
  @fn int putalabel(char where, char *overfrac, char underfrac, float charheight, float xmin, float xmax)
  @brief Puts a label at a viewgraph as a fraction

  Requires an open viewgraph, and the character height set to the axis
  numbering height. If both overfrac and underfrac are filled, then a fraction will be plotted. If one of both is an empty string (not NULL), then no fraction bar will be plotted.

  @param where      (char)   t, b, l, r position of the label
  @param overfrac   (char *) Text above the fraction bar
  @param underfrac  (char *) Text below the fraction bar
  @param charheight (float)  Character height relative to the actual character height, which should be set to the numbering height.
  @param dist       (float)  Distance from the axis baseline in units of the initial height
  @return (success) int putalabel: 1\\
          (error)   0

*/
/* ------------------------------------------------------------ */
int putalabel(char where, char *overfrac, char *underfrac, float charheight, float dist);

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
  @fn float pgp_logf10f(float x)
  @brief Returns the logarithm on basis 10 and LOGFAULT if x <= 0

  @param x (float) number to calculate the logarithm of
  @return logarithm of the number or, if impossible, LOGFAULT

*/
/* ------------------------------------------------------------ */
float pgp_logf10f(float x);

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/**
  @fn static int pgp_get_interparray(float *x, float *y, int npoints, int interptype, int nout, float **xout, float **yout)
  @brief Returns an array of interpolated values.

  Assumes an array of npoints x-y pairs as an input. Checks if x is strictly monotonically increasing. If *xout is not NULL, deallocates *xout. If *yout is not NULL, deallocates *yout. Allocates an array *xout of nout equally-spaced points with xout[0] = x[0] and xout[nout-1] = x[nout-1]. Allocates an array *yout of nout data points where the ith data point is an interpolation between x and y at the position xout[i]. interptype determines the interpolation type, which can be  PGP_I_LINEAR: linear; PGP_I_CSPLINE: cubic natural spline; PGP_I_AKIMA: Akima .

  @param x          (float * ) input data points, abscissae, must be strictly increasing
  @param y 	    (float * ) input data points, ordinate
  @param npoints    (  int   ) input number of data points
  @param interptype (  int   ) input interpolation type: PGP_I_LINEAR: linear; PGP_I_CSPLINE: cubic natural spline; PGP_I_AKIMA: Akima
  @param nout	    (  int   ) input number of output data points
  @param xout 	    (float **) output data array x-values (equally spaced)
  @param yout 	    (float **) output data array y-values

  @return logarithm of the number or, if impossible, LOGFAULT

*/
/* ------------------------------------------------------------ */
 static int pgp_get_interparray(float *x, float *y, int npoints, int interptype, int nout, float **xout, float **yout);

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* FUNCTION CODE */
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Opens the pgplot device as specified */
int pgp_opendev(char *device_name)
{
  int val;

  if (device_name)
    if (device_name[0]) {
      val = cpgbeg(0,device_name,1,1);

      /* Has to be set 0 for use in gipsy */
      cpgask(1);
      return val;
    }

  return 0;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Returns a default pgp_gdsc struct */
pgp_gdsc *pgp_gdsc_default(int nplots, int legendcols, int legendrows, float charheight)
{
  pgp_gdsc *gdsc;
  /* int i;  Control variable */
  float chrhtx; /* Horizontal character width in x-map units, 1 is whole page */
  float chrhty; /* Vertical character height in y-map units, 1 is whole page */
  float chrhtxref; /* Reference charheight */
  float scalexy; /* scale from x size of numbers to y size of numbers
		    (is the same physical size, but in terms of window
		    size they are different) */
  float pltheight;

  /* Allocate */
  if (!(gdsc = (pgp_gdsc *) malloc(sizeof(pgp_gdsc))))
    return NULL;

  /* Find the reference */
  cpgsch(1.0);
  cpgqcs(0, &chrhtxref, &chrhty);

  /* This is a scale */
  scalexy = chrhty/chrhtxref;

  /* Now, the character height in window units is like this */
  chrhtx = chrhtxref*charheight;

  /* Fill with default values */

  /* Number of plots */
  gdsc -> nplots = nplots;

  /* Symbols and descriptors in charheight, but at most one
     third??? of the plot height until the numbering is finished */

  /* The height of the plot is 1.0 at the start */
  pltheight = 1.0;

  /* The number is reduced by the axis description + 1 character */
  pltheight = pltheight-12.3*chrhty;

  /* Then, for each legend line, two characters */
  pltheight = pltheight-2.0*legendrows*chrhty;

  /* Then, each plot has only 1/nplots height */
  pltheight = pltheight/nplots;

  /* At the end, each plot has a margin of 4 heights at the top for the numbering */
  pltheight = pltheight-4.0*chrhty;

  /* We could express this as: */
  /* (1.0-12.3*chrhty-2.0*legendrows*chrhty)/nplots-4.0*chrhty */

  /* If this leaves room for 3 (???) numbers, then we apply the character height else we adjust the height */
  if (chrhty*3 > pltheight)
    chrhty = 1/(12.3+2*legendrows+7.0*nplots);

  chrhtx= chrhty/scalexy;

  gdsc -> numberheight = chrhtx/chrhtxref;
  gdsc -> symbolheight = 1.0;
  gdsc -> legendheight = 1.0;
  gdsc -> axdescheight = 1.0;

  /* legend */
  gdsc -> legendcols = legendcols;
  gdsc -> legendrows = legendrows;

  /* Horizontal and vertical borders are 2 characters */
  gdsc -> verbord = 2.0;
  gdsc -> horbord = 2.0;

  /* Stop of the axis numbering is 2 characters, only for the top of the graph */
  gdsc -> leftnum = 0.0;
  gdsc -> rightnum = 0.0;
  gdsc -> botnum = 0.0;
  gdsc -> topnum = 1.0;

  /* Margins */
  gdsc -> leftmargin = 7.0;
  gdsc -> rightmargin = 7.0;
  gdsc -> bottommargin = 1.8;
  gdsc -> topmargin = 2.5;

  /* Linewidths */
  gdsc -> graphlw = 1;
  gdsc -> boxlw = 1;

  /* reserve and return an initialised array for the axis labelling */
/*   if (!(gdsc -> logarcsx = (int *) malloc(gdsc -> nplots*sizeof(int)))){ */
/*     free(gdsc); */
/*     return 0; */
/*   } */
/*   for (i = 0; i < gdsc -> nplots; ++i) */
/*     gdsc -> logarcsx[i] = 0; */

  /* reserve and return an initialised array for the axis labelling */
/*   if (!(gdsc -> logarcsy = (int *) malloc(gdsc -> nplots*sizeof(int)))){ */
/*     free(gdsc -> logarcsx); */
/*     free(gdsc); */
/*     return 0; */
/*   } */
/*   for (i = 0; i < gdsc -> nplots; ++i) */
/*     gdsc -> logarcsy[i] = 0; */
  gdsc -> logarcsx = 0;
  gdsc -> logarcsy = 0;
  gdsc -> interptype_lines = PGP_I_LINEAR;
  gdsc -> interp_numlines = 500;

  return gdsc;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Make a polar diagram on an open device */
void pgp_polar(int ndots, float *xarray, float *yarray, int nlines, float *xlarray, float *ylarray, int nradii, float *radiis, float *lineangle, int linewidth, int circlewidth, int symbolwidth, int colour)
{
  int i;
  float rmax = 0;
  float twopointsx[2];
  float twopointsy[2];
  float xv1,yv1;
  float *radii;
  rmax = 0;

  /* Allocte memory for radii */
  if (!(radii = (float *) malloc(nradii*sizeof(float))))
    return;

  /* First find out about the extent of the graph */
  for (i = 0; i < ndots; ++i) {
    if (rmax < fabs(xarray[i]))
      rmax = fabs(xarray[i]);
    if (rmax < fabs(yarray[i]))
      rmax = fabs(yarray[i]);
  }

  for (i = 0; i < nlines; ++i) {
    if (rmax < fabs(xlarray[i]))
      rmax = fabs(xlarray[i]);
    if (rmax < fabs(ylarray[i]))
      rmax = fabs(ylarray[i]);
  }

  /* We make a local copy of the radii */

  for (i = 0; i < nradii; ++i) {
    radii[i] = fabs(radiis[i]);
    if (rmax < radii[i])
      rmax = radii[i];
  }

  /* Now scale rmax a bit */
  rmax = rmax*1.02;

  /* Define a drawing area, must be quadratic */
  cpgenv(-rmax, rmax, -rmax, rmax, 1, -2);

  /* Check the linewidths */
  if (linewidth < 1)
    linewidth = 1;
  if (linewidth > 201)
    linewidth = 201;
  if (circlewidth < 1)
    circlewidth = 1;
  if (circlewidth > 201)
    circlewidth = 201;
  if (symbolwidth < 1)
    symbolwidth = 1;
  if (symbolwidth > 201)
    symbolwidth = 201;

  /* Set the colour index foreground */
  cpgsci(1);

  /* Draw the circles */
  if (nradii > 0) {
    cpgslw(circlewidth);

    /* Make the circles open */
    cpgsfs(2);

    /* Chose a dashed style */
    cpgsls(4);

    for (i = 0; i < (nradii-1); ++i)
      cpgcirc(0,0, radii[i]);

    /* Draw the last circle, solid line with double width, solid */
    cpgsls(1);

    if (circlewidth > 100)
      circlewidth = 201;
    else
      circlewidth = 2*circlewidth;

    cpgslw(circlewidth);

    cpgcirc(0,0, radii[nradii-1]);

    /* Now do the Briggs style thingies */
    twopointsx[0] = 0;
    twopointsx[1] = 0;
    twopointsy[0] = radii[nradii-1];
    twopointsy[1] = radii[nradii-1]-0.05*radii[nradii-1];
    cpgline(2, twopointsx, twopointsy);

    twopointsx[0] = 0;
    twopointsx[1] = 0;
    twopointsy[0] = -radii[nradii-1];
    twopointsy[1] = -radii[nradii-1]+0.05*radii[nradii-1];
    cpgline(2, twopointsx, twopointsy);

    twopointsy[0] = 0;
    twopointsy[1] = 0;
    twopointsx[0] = radii[nradii-1];
    twopointsx[1] = radii[nradii-1]-0.05*radii[nradii-1];
    cpgline(2, twopointsx, twopointsy);

    twopointsy[0] = 0;
    twopointsy[1] = 0;
    twopointsx[0] = -radii[nradii-1];
    twopointsx[1] = -radii[nradii-1]+0.05*radii[nradii-1];
    cpgline(2, twopointsx, twopointsy);

    /* Check whether to plot a point */
    for (i = 0; i < (nradii); ++i) {
      if (radii[i] == 0){
	xv1 = yv1 = 0;
	cpgpt(1, &xv1, &yv1, -1);
      }
    }
    }
    /* Draw the reference thing */
    if ((lineangle)) {
      cpgslw(linewidth);
      twopointsy[0] = 0;
      twopointsy[1] = cos(DEGTORAD*(*lineangle))*radii[nradii-1];
      twopointsx[0] = 0;
      twopointsx[1] = -sin(DEGTORAD*(*lineangle))*radii[nradii-1];
    cpgline(2, twopointsx, twopointsy);
    }

  /* Plot the points */

  /* Set the linewidth */
  cpgslw(symbolwidth);

  /* Set the colour */
  cpgsci(colour);

 /* Plot points, number of points, vectors for x and y, symbol */
  if (ndots > 0)
  cpgpt(ndots, xarray, yarray, -1);

  /* Make a first dot */
  symbolwidth *= 2;
  if (symbolwidth > 201)
    symbolwidth = 210;
  cpgslw(symbolwidth);

  cpgpt(1, xarray, yarray, -1);

  /* Draw the lines */
  /* Set the linewidth */
  if (nlines > 1) {
  cpgslw(linewidth);
  cpgline(nlines, xlarray, ylarray);
  }

  /* Finished */
  return;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Opens a box in an opened window with a few specifications */
int pgp_openbox(pgp_gdsc *gdsc, int nplot, float xmin, float xmax, float ymin, float ymax, char *leftdeschi, char *leftdesclo, char *rightdeschi, char *rightdesclo, char *bottomdeschi, char *bottomdesclo, char *topdeschi, char *topdesclo, float lrzero, float lrscale, float btzero, float btscale, int logarcsx, int logarcsy)
{
  float plotheight;
  float chrhtxref;    /* Horizontal character width in x-map units, 1 is whole page */
  float chrhtyref;    /* Vertical character height in y-map units, 1 is whole page */
  /* float scalexy; */    /* scale from x size of numbers to y size of numbers (or inverse scale of x to y size of area */
  float dummy1, dummy2;

  int slw; /* Standard linewidth */

  float xlo; /* Box in window coordinates */
  float xhi;
  float ylo;
  float yhi;

  float xulo; /* Box in user coordinates */
  float xuhi;
  float yulo;
  float yuhi;

  float xlo2; /* Another box in window coordinates */
  float xhi2;
  float ylo2;
  float yhi2;

  float xulo2; /* Another box in user coordinates */
  float xuhi2;
  float yulo2;
  float yuhi2;

  float ax, ay, bx, by; /* linear coefficients from window to user coordinates */

  void (* boxingx)(const char *xopt, float xtick, int nxsub, const char *yopt, float ytick, int nysub);
  void (* boxingy)(const char *xopt, float xtick, int nxsub, const char *yopt, float ytick, int nysub);
  char steringxt[8],steringxb[8],steringyl[8],steringyr[8],steribgxt[8],steribgxb[8],steribgyl[8],steribgyr[8];

  /* First calculate the symbol height and width in map units
		    (is the same physical size, but in terms of window
		    size they are different) */

  /* Check if there is anything to draw */
  if (nplot > gdsc -> nplots || nplot < 1)
    return 0;

  /* Check the scaling */
  if ((logarcsx < 0) || (logarcsx > 3))
    return 0;

  /* per default the function to draw the box is cpgbox */
  boxingx = cpgbox;
  boxingy = cpgbox;

  /* Default numbering style */
  sprintf(steringyr,"MV");
  sprintf(steringyl,"NV");
  sprintf(steringxt,"M");
  sprintf(steringxb,"N");

  sprintf(steribgyr,"CTS");
  sprintf(steribgyl,"BTS");
  sprintf(steribgxt,"CTS");
  sprintf(steribgxb,"BTS");

  /* Default scaling */
  gdsc -> logarcsx = 0;
  gdsc -> logarcsy = 0;

  /* If a different scaling was requested, redefine labelling */
  if (logarcsx == 1) {
    gdsc -> logarcsx = 1;
/*     if (xmin <= 0) */
/*       return 0; */
/*     if (xmax <= 0) */
/*       return 0; */
/*     if ((xmin*btscale+btzero) <= 0) */
/*       return 0; */
/*     if ((xmax*btscale+btzero) <= 0) */
/*       return 0; */

    /* Do a trick */
    dummy1 = pgp_logf10f(xmin*btscale+btzero);
    dummy2 = pgp_logf10f(xmax*btscale+btzero);

    xmin = pgp_logf10f(xmin);
    xmax = pgp_logf10f(xmax);

    if (xmin == xmax) {
      btzero = dummy1-xmin;
      btscale = 0;
    }
    else {
      btscale = (dummy2-dummy1)/(xmax-xmin);
      btzero = dummy2-btscale*xmax;
    }
    sprintf(steringxt,"ML2");
    sprintf(steringxb,"NL2");
    sprintf(steribgxt,"CTSL");
    sprintf(steribgxb,"BTSL");
  }

  if (logarcsy == 1) {
    gdsc -> logarcsy = 1;
/*     if (ymin <= 0) */
/*       return 0; */
/*     if (ymax <= 0) */
/*       return 0; */
/*     if ((ymin*btscale+btzero) <= 0) */
/*       return 0; */
/*     if ((ymax*btscale+btzero) <= 0) */
/*       return 0; */

    /* Do a trick */
    dummy1 = pgp_logf10f(ymin*lrscale+lrzero);
    dummy2 = pgp_logf10f(ymax*lrscale+lrzero);

    ymin = pgp_logf10f(ymin);
    ymax = pgp_logf10f(ymax);

    if (ymin == ymax) {
      lrzero = dummy1-ymin;
      lrscale = 0;
    }
    else {
      lrscale = (dummy2-dummy1)/(ymax-ymin);
      lrzero = dummy2-lrscale*ymax;
    }
    sprintf(steringyr,"MVL2");
    sprintf(steringyl,"NVL2");
    sprintf(steribgyr,"CTSL");
    sprintf(steribgyl,"BTSL");
  }

  /* hms */
  if (logarcsx == 2) {
    gdsc -> logarcsx = 2;
    /* convert to seconds of time, input is interpreted as deg */
    btzero = btzero/240.0;
    btscale = btscale/240.0;
    xmin = xmin*240.0;
    xmax = xmax*240.0;
    boxingx = cpgtbox;
    sprintf(steringxt,"M");
    sprintf(steringxb,"NZH");
    sprintf(steribgxt,"CTS");
    sprintf(steribgxb,"BTSZXYH");
  }
  if (logarcsy == 2) {
    gdsc -> logarcsy = 2;

    /* convert to seconds of time, input is interpreted as deg */
    ymin = ymin*240.0;
    ymax = ymax*240.0;

    /* I'm really not sure why this is not the case */
/*      lrzero = lrzero/240.0; */
     lrscale = lrscale/240.0;
    boxingy = cpgtbox;
    sprintf(steringyr,"MV");
    sprintf(steringyl,"NVZXYH");
    sprintf(steribgyr,"CTS");
    sprintf(steribgyl,"BTSZXYH");
  }

  /* dms */
  if (logarcsx == 3) {
    gdsc -> logarcsx = 3;

    /* convert to seconds of time, input is interpreted as deg */
    xmin = xmin*3600.0;
    xmax = xmax*3600.0;

    /* I'm really not sure why this is not the case */
/*     btzero = btzero/3600.0; */
    btscale = btscale/3600.0;
    boxingx = cpgtbox;
    sprintf(steringxt,"M");
    sprintf(steringxb,"NZD");
    sprintf(steribgxt,"CTS");
    sprintf(steribgxb,"BTSZYD");
  }
  if (logarcsy == 3) {
    gdsc -> logarcsy = 3;

    /* convert to seconds of time, input is interpreted as deg */
    ymin = ymin*3600.0;
    ymax = ymax*3600.0;
/*     lrzero = lrzero/3600.0; */
    lrscale = lrscale/3600.0;
    boxingy = cpgtbox;
    sprintf(steringyr,"MV");
    sprintf(steringyl,"NVZYD");
    sprintf(steribgyr,"CTS");
    sprintf(steribgyl,"BTSZYD");
  }

  /* Ensure that minimum and maximum are different */
  if(xmin == xmax) {
    if (xmin != 0.0) {
      xmin = xmin*0.9;
      xmax = xmax*1.1;
    }
    else {
      xmin = -1.0;
      xmax = 1.0;
    }
  }

  if (logarcsx > 1) {
    if ((xmin-xmax) < 1.0) {
      xmin = xmin-1.0;
      xmax = xmax+1.0;
    }
  }

  if(ymin == ymax) {
    if (ymin != 0.0) {
      ymin = ymin*0.9;
      ymax = ymax*1.1;
    }
    else {
      ymin = -1.0;
      ymax = 1.0;
    }
  }

  if (logarcsy > 1) {
    if ((ymin-ymax) < 1.0) {
      ymin = ymin-1.0;
      ymax = ymax+1.0;
    }
  }

  /**********/
  /* scales */
  /**********/

  /* Standard linewidth */
  cpgqlw(&slw);

  /* Find the reference */
  cpgsch(gdsc -> numberheight);
  cpgqcs(0, &chrhtxref, &chrhtyref);

  /* The physical size in x is scalexy times the physical size in y */
  /* scalexy = chrhtyref/chrhtxref; */

  /* plot height */
  plotheight = (1.0 - 2.0*(gdsc -> legendrows+1)*gdsc -> legendheight*chrhtyref - 6*gdsc -> axdescheight*chrhtyref-6*chrhtyref)/gdsc -> nplots;

  /* Position of the plot, whole box */
  xlo = gdsc -> leftmargin*chrhtxref+3*gdsc -> axdescheight*chrhtxref;

  /* changed from 3*gdsc -> axdescheight to 4*gdsc -> axdescheight */
  xhi = 1.0-4*gdsc -> axdescheight*chrhtxref-gdsc -> rightmargin*chrhtxref;
  ylo = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot);
  yhi = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot-1);

  /* Check if the specification is wrong and return 0 if it is */
  if (xlo >= xhi)
    return 0;
  if (ylo >= yhi)
    return 0;

  /* Define the linear dependencies of the window scaling and the user scaling */
  ax = (xmax-xmin)/((xhi-gdsc -> horbord*gdsc -> symbolheight*chrhtxref) - (xlo+gdsc -> horbord*gdsc -> symbolheight*chrhtxref));
  ay =  (ymax-ymin)/((yhi-gdsc -> verbord*gdsc -> symbolheight*chrhtyref) - (ylo+gdsc -> verbord*gdsc -> symbolheight*chrhtyref));
  bx = xmin-ax*(xlo+gdsc -> horbord*gdsc -> symbolheight*chrhtxref);
  by = ymin-ay*(ylo+gdsc -> verbord*gdsc -> symbolheight*chrhtyref);

  /*************/
  /* Numbering */
  /*************/

  /* Set the symbol size (again) */
  cpgsch(gdsc -> numberheight);

  /* Draw the numbers axis by axis */

  /* First define the plot region */
  xlo2 = gdsc -> leftmargin*chrhtxref+3*gdsc -> axdescheight*chrhtxref+gdsc -> leftnum*chrhtxref;
  xhi2 = 1.0-gdsc -> rightmargin*chrhtxref-4*gdsc -> axdescheight*chrhtxref-gdsc -> rightnum*chrhtxref;
  ylo2 = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot)+gdsc -> botnum*chrhtyref;
  yhi2 = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot-1)-gdsc -> topnum*chrhtyref;

  xulo2 = ax*xlo2+bx;
  xuhi2 = ax*xhi2+bx;
  yulo2 = ay*ylo2+by;
  yuhi2 = ay*yhi2+by;

  /* Top and bottom numbering */
  if (nplot == 1) {
    cpgsvp(xlo2, xhi2, ylo, yhi);
    cpgswin(xulo2*btscale+btzero, xuhi2*btscale+btzero, yulo2*lrscale+lrzero, yuhi2*lrscale+lrzero);

    if ((topdeschi) || topdesclo)
      (*boxingx)(steringxt, 0.0, 0, "", 0.0, 0);
  }

  if (nplot == gdsc -> nplots) {
    cpgsvp(xlo2, xhi2, ylo, yhi);

    /* BUGFIX: changed from (...,yulo, yuhi) to (...,yulo2, yuhi2) */
    cpgswin(xulo2, xuhi2, yulo2, yuhi2);

    if ((bottomdeschi) || bottomdesclo)
      (*boxingx)(steringxb, 0.0, 0, "", 0.0, 0);
  }

  /* Left and right numbering */
  cpgsvp(xlo, xhi, ylo2, yhi2);
  cpgswin(xulo2*btscale+btzero, xuhi2*btscale+btzero, yulo2*lrscale+lrzero, yuhi2*lrscale+lrzero);

  if ((rightdeschi) || (rightdesclo))
    (*boxingy)("", 0.0, 0, steringyr, 0.0, 0);
  cpgsvp(xlo, xhi, ylo2, yhi2);
  cpgswin(xulo2, xuhi2, yulo2, yuhi2);

  if ((leftdeschi) || (leftdesclo))
    (*boxingy)("", 0.0, 0, steringyl, 0.0, 0);

  /***************/
  /* Box drawing */
  /***************/

  /* Set the line width */
  cpgslw(gdsc -> boxlw);

  /* Define an area of operation PGSVP (XMIN, XMAX, YMIN, YMAX), 0 to 1,  y from bottom to top */
  cpgsvp(xlo, xhi, ylo, yhi);

  xulo = ax*xlo+bx;
  xuhi = ax*xhi+bx;
  yulo = ay*ylo+by;
  yuhi = ay*yhi+by;

  cpgswin(xulo, xuhi, yulo, yuhi);

/* draw the box with the tickmarks of the lower and left panel */
  (*boxingx)(steribgxb, 0.0, 0, "", 0.0, 0);
  (*boxingy)("", 0.0, 0, steribgyl, 0.0, 0);

  /* Now change the user coordinates */
  cpgswin(xulo*btscale+btzero, xuhi*btscale+btzero, yulo*lrscale+lrzero, yuhi*lrscale+lrzero);

  /* draw the box with the tickmarks of the upper and right panel */
  (*boxingx)(steribgxt, 0.0, 0, "", 0.0, 0);
  (*boxingy)("", 0.0, 0, steribgyr, 0.0, 0);

  /* Reset the linewidth */
  cpgslw(slw);

  /*************/
  /* Labelling */
  /*************/

  /* Set character height (again) to the number value */
  cpgsch(gdsc -> numberheight);

  /* Label top */
  if (nplot == 1)
    putalabel('t', topdeschi, topdesclo, gdsc -> axdescheight, gdsc -> topmargin);

  /* Label bottom */
  if (nplot == gdsc -> nplots)
    putalabel('b', bottomdeschi, bottomdesclo, gdsc -> axdescheight, gdsc -> bottommargin);

  /* Label left and right axis */
  putalabel('l', leftdeschi, leftdesclo, gdsc -> axdescheight, gdsc -> leftmargin);
  putalabel('r', rightdeschi, rightdesclo, gdsc -> axdescheight, gdsc -> rightmargin);

  /* At the end we put the viewport and the window at the right values */
  switch (logarcsx) {
  case 2:
    xmin = xmin/240.0;
    xmax = xmax/240.0;
    break;
  case 3:
    xmin = xmin/3600.0;
    xmax = xmax/3600.0;
    break;
  }
  switch (logarcsy) {
  case 2:
    ymin = ymin/240.0;
    ymax = ymax/240.0;
    break;
  case 3:
    ymin = ymin/3600.0;
    ymax = ymax/3600.0;
    break;
  }

  /* Redefine */
  xlo = gdsc -> leftmargin*chrhtxref+3*gdsc -> axdescheight*chrhtxref;
  xhi = 1.0-4*gdsc -> axdescheight*chrhtxref-gdsc -> rightmargin*chrhtxref;
  ylo = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot);
  yhi = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(nplot-1);

  ax = (xmax-xmin)/((xhi-gdsc -> horbord*gdsc -> symbolheight*chrhtxref) - (xlo+gdsc -> horbord*gdsc -> symbolheight*chrhtxref));
  ay =  (ymax-ymin)/((yhi-gdsc -> verbord*gdsc -> symbolheight*chrhtyref) - (ylo+gdsc -> verbord*gdsc -> symbolheight*chrhtyref));
  bx = xmin-ax*(xlo+gdsc -> horbord*gdsc -> symbolheight*chrhtxref);
  by = ymin-ay*(ylo+gdsc -> verbord*gdsc -> symbolheight*chrhtyref);

  xulo = ax*xlo+bx;
  xuhi = ax*xhi+bx;
  yulo = ay*ylo+by;
  yuhi = ay*yhi+by;

  /* Do the redefinition */
  cpgsvp(xlo, xhi, ylo, yhi);
  cpgswin(xulo, xuhi, yulo, yuhi);

  return 1;
}
/* ------------------------------------------------------------ */

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts a label at a viewgraph as a fraction */

int putalabel(char where, char *overfrac, char *underfrac, float charheight, float dist)
{
  float oheight;    /* Character height at call time */
  int fraclength;
  char *fracbar;
  float ovpos;
  float unpos;
  float fracpos;

  char format[3];

  float frlen1;
  float frlen2;
  float xbox;
  float ybox;
  float chrhtxref;
  float chrhtyref;

  /* Check if there is anything to put */
  if ((overfrac == NULL) && (underfrac == NULL))
    return 1;

  /***********/
  /* scaling */
  /***********/

  /* Find the reference character height */
  cpgqch(&oheight);

  /* Set the character height to the required value */
  cpgsch(oheight*charheight);

  /* If the strings are empty we return a success */
  if (!(*overfrac)) {
    if (!(*underfrac))
      return 1;
    /* If only one of them is empty, we make the one not empty overfrac */
    else {
      fracbar = overfrac;
      overfrac = underfrac;
      underfrac = fracbar;
    }
  }

  /* A format string is required */
  switch (where) {
  case 'l':
    sprintf(format,"L");

    /* Calculate the maximum possible distance from the edge of the graph */
    ovpos = dist/charheight+2.0;
    fracpos = dist/charheight+1.4;
    unpos = dist/charheight;
    break;
  case 'r':
    sprintf(format,"R");

    /* Calculate the maximum possible distance from the edge of the graph */
    ovpos = dist/charheight+1;
    fracpos = dist/charheight+1.6;
    unpos = dist/charheight+3;
    break;

  case 't':
    sprintf(format,"T");
    ovpos = dist/charheight+2.0;
    fracpos = dist/charheight/charheight+1.4;
    unpos = dist/charheight/charheight+0.0;
    break;

  case 'b':
    sprintf(format,"B");
    ovpos = dist/charheight+1.0;
    fracpos = dist/charheight+1.6;
    unpos = dist/charheight+3.0;
    break;

  default:
    return 0;
  }

  /* To get the length of the fraction bar, first find a bounding box at an arbitrary position */
  cpglen(2, overfrac, &xbox, &ybox);
  frlen1 = ybox;
  cpglen(2, underfrac, &xbox, &ybox);
  frlen2 = ybox;
  frlen1= (frlen1 > frlen2)?frlen1:frlen2;

  /* Now frlen1 is the length of the string in world coordinates of the x axis, we get the character width of this axis in the same units */
  cpgqcs(2, &chrhtxref, &chrhtyref);

  /* Then, the length of the fraction bar is the uprounded fraction */
  fraclength = ((int) (frlen1/chrhtyref))+3;

  if (!(fracbar = (char *) malloc((fraclength+1)*sizeof(char))))
    return 0;

  fracbar[fraclength] = '\0';

  while (fraclength > 0) {
    fracbar[fraclength-1] = '_';
    --fraclength;
  }

  /************/
  /* plotting */
  /************/

  /* If there is a fraction bar to plot do this */
  if (*underfrac) {
  /* First plot the character string above the fraction bar */
  cpgmtxt(format, ovpos, 0.5, 0.5, overfrac);

  /* Now draw the fraction bar */
  cpgmtxt(format, fracpos, 0.5, 0.5, fracbar);

  /* At last put the argument below the fraction bar */
  cpgmtxt(format, unpos, 0.5, 0.5, underfrac);
  }

  /* No fraction bar */
  else {
    cpgmtxt(format, unpos, 0.5, 0.5, overfrac);
  }

  /* Deallocate */
  free(fracbar);

  /* Set the character height to the original value */
  cpgsch(oheight);

  return 1;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts a string to the specified legend position */
int pgp_legend(pgp_gdsc *gdsc, int nhor, int nver, char *string)
{
  float chrhtxref;    /* Horizontal character width in x-map units, 1 is whole page */
  float chrhtyref;    /* Vertical character height in y-map units, 1 is whole page */
  /* float charwidth; */
  float posix;
  float posiy;
  float plotheight;
  /* float boxwidth; */
  float xlo; /* Box in window coordinates */
  float xhi;
  float ylo;
  float yhi;

  /* First check if the intended legend position is available */
  if (nver < 1 || nver > gdsc -> legendrows || nhor < 1 || nhor > gdsc -> legendcols) {
    return 0;
  }

  /* Then open a box that nearly spans the whole surface at the bottom position of the graph */

  /* Find the reference */
  cpgsch(gdsc -> numberheight);
  cpgqcs(0, &chrhtxref, &chrhtyref);

  /* plot height */
  plotheight = (1.0 - 2.0*(gdsc -> legendrows+1)*gdsc -> legendheight*chrhtyref - 6*gdsc -> axdescheight*chrhtyref-6*chrhtyref)/gdsc -> nplots;
  ylo = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(gdsc -> nplots);
  yhi = 1.0-4*gdsc -> axdescheight*chrhtyref-gdsc -> topmargin*chrhtyref-plotheight*(gdsc -> nplots-1);

  /* At the left hand the box has the distance of one symbol to the border of the plotting surface */
  cpgsch(gdsc -> numberheight*gdsc -> legendheight);
  cpgqcs(0, &chrhtxref, &chrhtyref);

  xlo = chrhtxref;
  xhi = 1.0-chrhtxref;

  /* Check if the specification is wrong and return 0 if it is */
  if (xlo >= xhi)
    return 0;
  if (ylo >= yhi)
    return 0;

  /* Define the surface */
  cpgsvp(xlo, xhi, ylo, yhi);

  /* This is the length of a box element in units of the x-axis, spacing is one character between rows */
  /* boxwidth = (1-(gdsc -> legendcols-1)*charwidth)/gdsc -> legendcols; */

  /* This is the position of the string in x along the axis */
  posix = (nhor-0.5)/gdsc -> legendcols;

  /* Now in y */
  posiy = gdsc -> bottommargin/gdsc -> legendheight+3*gdsc -> axdescheight/gdsc -> legendheight+2*nver;

  /* Put it on the graph */
  cpgmtxt("B", posiy, posix, 0.5, string);

  /* Reset the character height */
  cpgsch(gdsc -> numberheight);

  return 1;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts markers on a predefined viewport */
int pgp_marker(pgp_gdsc *gdsc, int npoints, float *xpos, float *ypos, int colour, int empty, int symbol, float sizer)
{
  int colourbef;
  float chrhtxref, chrhtyref;
  int intlw;
  float *xposd = NULL;
  float *yposd = NULL;
  int i;

  /* Allocate space if necessary */
  if (gdsc -> logarcsx == 1) {
    if (!(xposd = (float *) malloc(npoints*sizeof(float))))
      return 1;
    for (i=0; i < npoints; ++i) {
	xposd[i] = pgp_logf10f(xpos[i]);
    }
    xpos = xposd;
  }
  if (gdsc -> logarcsy == 1) {
    if (!(yposd = (float *) malloc(npoints*sizeof(float)))) {
      if ((xposd))
	free(xposd);
      return 2;
    }
    for (i=0; i < npoints; ++i) {
	yposd[i] = pgp_logf10f(ypos[i]);
    }
    ypos = yposd;
  }

  /* Set the character height */
  if (symbol == -1)
    cpgsch(gdsc -> numberheight*gdsc -> symbolheight*sizer);
  else if (symbol == 3)
    cpgsch(gdsc -> numberheight*gdsc -> symbolheight*2*sizer*1.2);
  else
    cpgsch(gdsc -> numberheight*gdsc -> symbolheight*2*sizer);

  /* Inquire the character height in inches */
  cpgqcs(1, &chrhtxref, &chrhtyref);

  /* Set the linewidth as an equivalent */
  if (symbol == -1) {
    intlw = chrhtxref*200.0>201.0?201: (int) (chrhtxref*200.0);
  }
  else
    intlw = 1;
  cpgslw(intlw);

  /* Store the colour index */
  cpgqci(&colourbef);

  /* Set the colour index */
  cpgsci(colour);

 /* Plot points, number of points, vectors for x and y, symbol */
  cpgpt(npoints, xpos, ypos, symbol);

  /* If required, plot the same list on top of that with background colour and with a radius that is smaller */
  if ((empty)) {
    cpgslw(((intlw-4*gdsc -> graphlw) > 0)?(intlw-4*gdsc -> graphlw):0);
    cpgsci(0);
    cpgpt(npoints, xpos, ypos, symbol);
  }

  /* Reset the linewidth */
  cpgslw(gdsc -> graphlw);

  /* Reset the colour index */
  cpgsci(colourbef);

  /* Reset the character height */
  cpgsch(gdsc -> numberheight);

  /* Reset the linewidth */
  cpgslw(gdsc -> boxlw);

  /* Deallocate */

  if (gdsc -> logarcsy == 1) {
    free(yposd);
  }

  if (gdsc -> logarcsx == 1) {
    free(xposd);
  }

  return 0;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts lines on a predefined viewport */
int pgp_lines(pgp_gdsc *gdsc, int npoints, float *xpos, float *ypos, int colour)
{
  int i;
  int colourbef;
  float *xposd = NULL;
  float *yposd = NULL;

  if (gdsc -> interptype_lines != PGP_I_LINEAR) {
    if (pgp_get_interparray(xpos, ypos, npoints, gdsc -> interptype_lines, gdsc -> interp_numlines+1, &xposd, &yposd))
      goto error;
    npoints = gdsc -> interp_numlines+1;
    xpos = xposd;
    ypos = yposd;
  }

  /* Allocate space if necessary */
  if (gdsc -> logarcsx == 1) {
    if (!xposd) {
      if (!(xposd = (float *) malloc(npoints*sizeof(float))))
	goto error;
    }
    for (i=0; i < npoints; ++i) {
	xposd[i] = pgp_logf10f(xpos[i]);
    }
    xpos = xposd;
  }
  if (gdsc -> logarcsy == 1) {
    if (!yposd) {
      if (!(yposd = (float *) malloc(npoints*sizeof(float))))
	goto error;
    }
    for (i=0; i < npoints; ++i) {
	yposd[i] = pgp_logf10f(ypos[i]);
    }
    ypos = yposd;
  }

  /* Store the colour index */
  cpgqci(&colourbef);

  /* Set the colour index */
  cpgsci(colour);

  /* Set the linewidth */
  cpgslw(gdsc -> graphlw);

 /* Plot lines, number of points, vectors for x and y, symbol */
  cpgline(npoints, xpos, ypos);

  /* Reset the linewidth */
  cpgslw(gdsc -> graphlw);

  /* Reset the colour index */
  cpgsci(colourbef);

  /* Reset the linewidth */
  cpgslw(gdsc -> boxlw);

  /* Deallocate */
  if ((xposd))
    free(xposd);
  if ((yposd))
    free(yposd);
  return 0;

 error:
  if ((xposd))
    free(xposd);
  if ((yposd))
    free(yposd);
  return 1;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts lines on a predefined viewport */
int pgp_bars(pgp_gdsc *gdsc, int npoints, float *xpos, float *ypos, float width, int colour)
{
  int colourbef,i;
  float *xposd = NULL;
  float *yposd = NULL;

  /* Allocate space if necessary */
  if (gdsc -> logarcsx == 1) {
    if (!(xposd = (float *) malloc(npoints*sizeof(float))))
      return 1;
    for (i=0; i < npoints; ++i) {
      xposd[i] = pgp_logf10f(xpos[i]);
    }
    xpos = xposd;
  }
  if (gdsc -> logarcsy == 1) {
    if (!(yposd = (float *) malloc(npoints*sizeof(float)))) {
      if ((xposd))
	free(xposd);
      return 2;
    }
    for (i=0; i < npoints; ++i) {
	yposd[i] = pgp_logf10f(ypos[i]);
    }
    ypos = yposd;
  }

  /* Store the colour index */
  cpgqci(&colourbef);

  /* Set the colour index */
  cpgsci(colour);

  /* Set the linewidth */
  cpgslw(gdsc -> graphlw);

 /* Plot lines, number of points, vectors for x and y, symbol */
  for (i = 0; i < npoints; ++i)
  cpgerr1(5, xpos[i], ypos[i], width/2, 0.0);

  /* Reset the linewidth */
  cpgslw(gdsc -> graphlw);

  /* Reset the colour index */
  cpgsci(colourbef);

  /* Reset the linewidth */
  cpgslw(gdsc -> boxlw);

  /* Deallocate */
  if (gdsc -> logarcsy == 1) {
    free(yposd);
  }

  if (gdsc -> logarcsx == 1) {
    free(xposd);
  }

  return 0;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Puts lines on a predefined viewport */
int pgp_errby(pgp_gdsc *gdsc, int npoints, float *xpos, float *ypos, float *err, int colour)
{
  int colourbef, i;
  float *xposd = NULL;
  float *yposdl = NULL;
  float *yposdr = NULL;

  /* Allocate space if necessary */
  if (gdsc -> logarcsx == 1) {
    if (!(xposd = (float *) malloc(npoints*sizeof(float))))
      return 1;
    for (i=0; i < npoints; ++i) {
	xposd[i] = pgp_logf10f(xpos[i]);
    }
    xpos = xposd;
  }

  if (gdsc -> logarcsy == 1) {
    if (!(yposdl = (float *) malloc(npoints*sizeof(float)))) {
      if ((xposd))
	free(xposd);
      return 1;
    }
    if (!(yposdr = (float *) malloc(npoints*sizeof(float)))) {
      if ((xposd))
	free(xposd);
      free(yposdl);
      return 1;
    }
    for (i=0; i < npoints; ++i) {
	yposdr[i] = pgp_logf10f(ypos[i]+err[i]);
      }
    for (i=0; i < npoints; ++i) {
	yposdl[i] = pgp_logf10f(ypos[i]-err[i]);
    }
  }

  /* Store the colour index */
  cpgqci(&colourbef);

  /* Set the colour index */
  cpgsci(colour);

  /* Set the linewidth */
  cpgslw(gdsc -> graphlw);

 /* Plot lines, number of points, vectors for x and y, symbol */
  if (gdsc -> logarcsy == 1)
    cpgerry(npoints, xpos, yposdl, yposdr, 2.0);
  else
    cpgerrb(6, npoints, xpos, ypos, err, 2.0);

  /* Reset the linewidth */
  cpgslw(gdsc -> graphlw);

  /* Reset the colour index */
  cpgsci(colourbef);

  /* Reset the linewidth */
  cpgslw(gdsc -> boxlw);

  /* Clean up */
  if (gdsc -> logarcsx == 1)
    free(xposd);
  if (gdsc -> logarcsy == 1) {
    free(yposdl);
    free(yposdr);
  }

  return 0;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* End the plotting and close the device */
void pgp_end()
{
  cpgclos();
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Returns the logarithm on basis 10 and LOGFAULT if x <= 0 */
float pgp_logf10f(float x)
{
  if (x > 0)
    return log10f(x);
  else
    return LOGFAULT;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/* Returns an array of interpolated values. */

static int pgp_get_interparray(float *x, float *y, int npoints, int interptype, int noutpoints, float **xout, float **yout)
{
  int i, pgp_get_interparray = 0;
  double *xdouble = NULL, *ydouble = NULL; /* input is float, but gsl requires double */
  gsl_interp *gsl_interpv = NULL; /* gsl structs */
  gsl_interp_accel *gsl_interp_accelv = NULL; /* gsl structs */
  const gsl_interp_type * intytype = NULL; /* interpolation type, internal to gsl */
  float dx;

  /* At least two data points */
  if (npoints < 2) {
    pgp_get_interparray |= 1;
    goto error;
  }

  /* Wipe the input, gone ... */
  if (*xout) {
    free(*xout);
    *xout = NULL;
  }
  if (*yout) {
    free(*yout);
    *yout = NULL;
  }

  /* Allocate arrays */
  if (!(xdouble = (double *) malloc(npoints*sizeof(double)))) {
    pgp_get_interparray |= 2;
    goto error;
  }
  if (!(ydouble = (double *) malloc(npoints*sizeof(double)))) {
    pgp_get_interparray |= 2;
    goto error;
  }
  for (i = 0; i < npoints; ++i) {
    xdouble[i] = x[i];
    ydouble[i] = y[i];
  }
  if (!(*xout = (float *) malloc(noutpoints*sizeof(float)))) {
    pgp_get_interparray |= 2;
    goto error;
  }
  if (!(*yout = (float *) malloc(noutpoints*sizeof(float)))) {
    pgp_get_interparray |= 2;
    goto error;
  }

  switch(interptype) {
    case PGP_I_LINEAR:
      intytype = gsl_interp_linear;
      break;

    case PGP_I_CSPLINE:
      if (npoints > 2)
	intytype = gsl_interp_cspline;
      else {
	printf("Must have at least 3 radii for spline, using linear\n");
	intytype = gsl_interp_linear;
      }
      break;

    case PGP_I_AKIMA:
      if (npoints > 4) {
	intytype = gsl_interp_akima;
      }
      else {
	printf("Must have at least 5 radii for Akima\n");
	if (npoints > 2) {
	  printf("Using natural cubic spline\n");
	  intytype = gsl_interp_cspline;
	}
	else {
	  printf("Using linear\n");
	  intytype = gsl_interp_linear;
	}
      }
      break;

    /* No real default, but we set linear */
    default:
      intytype = gsl_interp_linear;
      break;
    }

  /* Attempt to allocate all gsl structures */
  if (!(gsl_interpv = gsl_interp_alloc(intytype, npoints))) {
    pgp_get_interparray |= 2;
    goto error;
  }
  if (!(gsl_interp_accelv = gsl_interp_accel_alloc())) {
    pgp_get_interparray |= 2;
    goto error;
  }

  gsl_interp_init(gsl_interpv, xdouble, ydouble, npoints);

  (*xout)[0] = x[0];
  (*yout)[0] = y[0];

  dx = (x[npoints-1]-x[0])/((float)(noutpoints-1));

  /* Bugfix: rounding error could lead to failure of interpolation */
  --noutpoints;

  for (i = 1; i < noutpoints; ++i) {
    (*xout)[i] = (*xout)[i-1]+dx;
    (*yout)[i] = gsl_interp_eval(gsl_interpv, xdouble, ydouble, (*xout)[i], gsl_interp_accelv);
  }
  (*xout)[i] = x[npoints-1];
  (*yout)[i] = y[npoints-1];

  /* last points should be identical regardless of rounding errors */
    free(xdouble);
    free(ydouble);
    gsl_interp_free(gsl_interpv);
    gsl_interp_accel_free(gsl_interp_accelv);

  return pgp_get_interparray;

 error:
  if (xdouble)
    free(xdouble);
  if (ydouble)
    free(ydouble);
  if (*xout) {
    free(*xout);
    *xout = NULL;
  }
  if (*yout) {
    free(*yout);
    *yout = NULL;
  }
  if (gsl_interpv)
    gsl_interp_free(gsl_interpv);
  if (gsl_interp_accelv)
    gsl_interp_accel_free(gsl_interp_accelv);

  return pgp_get_interparray;
}

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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   $Log: pgp.c,v $
   Revision 1.14  2011/05/25 22:25:26  jozsa
   Left work

   Revision 1.13  2011/05/10 00:30:16  jozsa
   Left work

   Revision 1.12  2009/05/26 07:56:41  jozsa
   Left work

   Revision 1.11  2007/08/22 15:58:43  gjozsa
   Left work

   Revision 1.10  2006/11/22 14:16:21  gjozsa
   Bugfix concerning RASH and horizontal/vertical lines

   Revision 1.9  2006/11/03 10:49:02  gjozsa
   Introduced logarithmic scaling, hms dms

   Revision 1.8  2006/10/05 12:39:58  gjozsa
   nasty bugfix

   Revision 1.7  2006/04/03 11:47:57  gjozsa
   Left work

   Revision 1.6  2006/02/08 16:13:17  gjozsa
   Extended graphics routines

   Revision 1.5  2005/10/12 14:51:25  gjozsa
   First point gets thick in the polar plot

   Revision 1.4  2005/10/12 09:53:09  gjozsa
   Included polar plot facility

   Revision 1.3  2005/06/24 12:00:29  gjozsa
   Left work

   Revision 1.2  2005/05/24 10:42:24  gjozsa
   Left work

   Revision 1.1  2005/05/17 12:58:15  gjozsa
   Added to cvs control

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

max_stars_repo_name: kernsuite-debian/tirific

max_issues_repo_name: gigjozsa/tirific

max_forks_repo_name: kernsuite-debian/tirific


text:
/* linalg/bidiag.c
 *
 * Copyright (C) 2001, 2007 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Factorise a matrix A into
 *
 * A = U B V^T
 *
 * where U and V are orthogonal and B is upper bidiagonal.
 *
 * On exit, B is stored in the diagonal and first superdiagonal of A.
 *
 * U is stored as a packed set of Householder transformations in the
 * lower triangular part of the input matrix below the diagonal.
 *
 * V is stored as a packed set of Householder transformations in the
 * upper triangular part of the input matrix above the first
 * superdiagonal.
 *
 * The full matrix for U can be obtained as the product
 *
 *       U = U_1 U_2 .. U_N
 *
 * where
 *
 *       U_i = (I - tau_i * u_i * u_i')
 *
 * and where u_i is a Householder vector
 *
 *       u_i = [0, .. , 0, 1, A(i+1,i), A(i+3,i), .. , A(M,i)]
 *
 * The full matrix for V can be obtained as the product
 *
 *       V = V_1 V_2 .. V_(N-2)
 *
 * where
 *
 *       V_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [0, .. , 0, 1, A(i,i+2), A(i,i+3), .. , A(i,N)]
 *
 * See Golub & Van Loan, "Matrix Computations" (3rd ed), Algorithm 5.4.2
 *
 * Note: this description uses 1-based indices. The code below uses
 * 0-based indices
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

int
gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V)
{
  if (A->size1 < A->size2)
    {
      GSL_ERROR ("bidiagonal decomposition requires M>=N", GSL_EBADLEN);
    }
  else if (tau_U->size  != A->size2)
    {
      GSL_ERROR ("size of tau_U must be N", GSL_EBADLEN);
    }
  else if (tau_V->size + 1 != A->size2)
    {
      GSL_ERROR ("size of tau_V must be (N - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t M = A->size1;
      const size_t N = A->size2;
      gsl_vector * tmp = gsl_vector_alloc(M);
      size_t j;

      for (j = 0 ; j < N; j++)
        {
          /* apply Householder transformation to current column */
          gsl_vector_view v = gsl_matrix_subcolumn(A, j, j, M - j);
          double tau_j = gsl_linalg_householder_transform (&v.vector);

          /* apply the transformation to the remaining columns */
          if (j + 1 < N)
            {
              gsl_matrix_view m = gsl_matrix_submatrix (A, j, j + 1, M - j, N - j - 1);
              gsl_vector_view work = gsl_vector_subvector(tau_U, j, N - j - 1);
              gsl_linalg_householder_left (tau_j, &v.vector, &m.matrix, &work.vector);
            }

          gsl_vector_set (tau_U, j, tau_j);

          /* apply Householder transformation to current row */
          if (j + 1 < N)
            {
              v = gsl_matrix_subrow (A, j, j + 1, N - j - 1);
              tau_j = gsl_linalg_householder_transform (&v.vector);

              /* apply the transformation to the remaining rows */
              if (j + 1 < M)
                {
                  gsl_matrix_view m = gsl_matrix_submatrix (A, j + 1, j + 1, M - j - 1, N - j - 1);
                  gsl_vector_view work = gsl_vector_subvector(tmp, 0, M - j - 1);
                  gsl_linalg_householder_right (tau_j, &v.vector, &m.matrix, &work.vector);
                }

              gsl_vector_set (tau_V, j, tau_j);
            }
        }

      gsl_vector_free(tmp);

      return GSL_SUCCESS;
    }
}

/* Form the orthogonal matrices U, V, diagonal d and superdiagonal sd
   from the packed bidiagonal matrix A */

int
gsl_linalg_bidiag_unpack (const gsl_matrix * A,
                          const gsl_vector * tau_U,
                          gsl_matrix * U,
                          const gsl_vector * tau_V,
                          gsl_matrix * V,
                          gsl_vector * diag,
                          gsl_vector * superdiag)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  const size_t K = GSL_MIN(M, N);

  if (M < N)
    {
      GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
    }
  else if (tau_U->size != K)
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (tau_V->size + 1 != K)
    {
      GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
    }
  else if (U->size1 != M || U->size2 != N)
    {
      GSL_ERROR ("size of U must be M x N", GSL_EBADLEN);
    }
  else if (V->size1 != N || V->size2 != N)
    {
      GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
    }
  else if (diag->size != K)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (superdiag->size + 1 != K)
    {
      GSL_ERROR ("size of subdiagonal must be (diagonal size - 1)", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      /* Copy diagonal into diag */

      for (i = 0; i < N; i++)
        {
          double Aii = gsl_matrix_get (A, i, i);
          gsl_vector_set (diag, i, Aii);
        }

      /* Copy superdiagonal into superdiag */

      for (i = 0; i < N - 1; i++)
        {
          double Aij = gsl_matrix_get (A, i, i+1);
          gsl_vector_set (superdiag, i, Aij);
        }

      /* Initialize V to the identity */

      gsl_matrix_set_identity (V);

      for (i = N - 1; i-- > 0;)
        {
          /* Householder row transformation to accumulate V */
          gsl_vector_const_view h = gsl_matrix_const_subrow (A, i, i + 1, N - i - 1);
          double ti = gsl_vector_get (tau_V, i);
          gsl_matrix_view m = gsl_matrix_submatrix (V, i + 1, i + 1, N- i - 1, N - i - 1);
          gsl_vector_view work = gsl_matrix_subrow(U, 0, 0, N - i - 1);

          gsl_linalg_householder_left (ti, &h.vector, &m.matrix, &work.vector);
        }

      /* Initialize U to the identity */

      gsl_matrix_set_identity (U);

      for (j = N; j-- > 0;)
        {
          /* Householder column transformation to accumulate U */
          gsl_vector_const_view h = gsl_matrix_const_subcolumn (A, j, j, M - j);
          double tj = gsl_vector_get (tau_U, j);
          gsl_matrix_view m = gsl_matrix_submatrix (U, j, j, M - j, N - j);

          gsl_linalg_householder_hm (tj, &h.vector, &m.matrix);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
                           gsl_vector * tau_U,
                           gsl_vector * tau_V,
                           gsl_matrix * V)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  const size_t K = GSL_MIN(M, N);

  if (M < N)
    {
      GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
    }
  else if (tau_U->size != K)
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (tau_V->size + 1 != K)
    {
      GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
    }
  else if (V->size1 != N || V->size2 != N)
    {
      GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      /* Initialize V to the identity */

      gsl_matrix_set_identity (V);

      for (i = N - 1; i-- > 0;)
        {
          /* Householder row transformation to accumulate V */
          gsl_vector_const_view r = gsl_matrix_const_row (A, i);
          gsl_vector_const_view h =
            gsl_vector_const_subvector (&r.vector, i + 1, N - (i+1));

          double ti = gsl_vector_get (tau_V, i);

          gsl_matrix_view m =
            gsl_matrix_submatrix (V, i + 1, i + 1, N-(i+1), N-(i+1));

          gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
        }

      /* Copy superdiagonal into tau_v */

      for (i = 0; i < N - 1; i++)
        {
          double Aij = gsl_matrix_get (A, i, i+1);
          gsl_vector_set (tau_V, i, Aij);
        }

      /* Allow U to be unpacked into the same memory as A, copy
         diagonal into tau_U */

      for (j = N; j-- > 0;)
        {
          /* Householder column transformation to accumulate U */
          double tj = gsl_vector_get (tau_U, j);
          double Ajj = gsl_matrix_get (A, j, j);
          gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M-j, N-j);

          gsl_vector_set (tau_U, j, Ajj);
          gsl_linalg_householder_hm1 (tj, &m.matrix);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
                            gsl_vector * diag,
                            gsl_vector * superdiag)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  const size_t K = GSL_MIN(M, N);

  if (diag->size != K)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (superdiag->size + 1 != K)
    {
      GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* Copy diagonal into diag */

      for (i = 0; i < K; i++)
        {
          double Aii = gsl_matrix_get (A, i, i);
          gsl_vector_set (diag, i, Aii);
        }

      /* Copy superdiagonal into superdiag */

      for (i = 0; i < K - 1; i++)
        {
          double Aij = gsl_matrix_get (A, i, i+1);
          gsl_vector_set (superdiag, i, Aij);
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


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

#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>

#define N (4096*4*16)

double g (double *x, int dim);

int
main (void)
{
  double sum;
  double s = 12.;
  int dim = 4;
  double x[4];

  gsl_rng *r;

  r = gsl_rng_alloc (gsl_rng_taus2);
  gsl_rng_set (r, 1UL);

  long i, j, nn;

  printf ("# Iter        Npoints   Volume     AbsErr\n");
  nn = N;
  for (j = 0; j < 16; j++)
    {
      sum = 0.;
      for (i = 0; i < nn; i++)
	{
	  // 4-dimensional random point
	  for (int k = 0; k < dim; k++)
	    {
	      x[k] = s * gsl_rng_uniform (r) - s / 2;
	    }
	  sum += g (x, dim);
	}
      double integ = (sum * pow (s, (double) dim)) / nn;
      double exact = M_PI * M_PI;

      printf ("%3ld  %16ld   %f   %f\n", j, nn, integ, fabs (exact - integ));
      nn *= 2;
    }

  gsl_rng_free (r);

  return 0;
}

double
g (double *x, int dim)
{
  double r2 = 0., q;

  for (int i = 0; i < dim; i++)
    {
      q = x[i];
      r2 += q * q;
    }

  return exp (-r2);
}


text:
// Created by BlurringShadow at 2021-02-27-下午 10:24

#pragma once

#include <algorithm>
#include <chrono>

#ifdef __cpp_lib_concepts
#include <concepts>
#include <ranges>
#endif

#include <execution>
#include <optional>
#include <string>
#include <type_traits>
#include <utility>

#include <nameof.hpp>
#include <boost/type_traits.hpp>
#include <entt/entt.hpp>
#include <fmt/ostream.h>
#include <glm/ext.hpp>
#include <gsl/gsl>
#include <range/v3/all.hpp>
#include <rxcpp/rx.hpp>
#include <tl/expected.hpp>

using namespace std::literals;

namespace polycumic::utility
{
#ifndef __cpp_lib_concepts
#define SFINAE(...) std::enable_if_t<__VA_ARGS__>* = nullptr
#endif
    using tl::expected;

    template<typename T>
    struct auto_cast
    {
        T&& t;

        template<typename U>
#ifdef __cpp_lib_concepts
        requires requires { static_cast<std::decay_t<U>>(std::forward<T>(t)); }
#endif
        [[nodiscard]] constexpr operator
        U() noexcept(std::is_nothrow_constructible_v<T, std::decay_t<U>>)
        {
            return static_cast<std::decay_t<U>>(std::forward<T>(t));
        }
    };

    template<typename T>
    auto_cast(T&& t) -> auto_cast<T>;

    template<typename T, typename U, typename Comp>
    constexpr T& set_if(T& left, U&& right, Comp comp)
    {
        if(comp(right, left)) left = std::forward<U>(right);
        return left;
    }

    template<typename T, typename U>
    constexpr T& set_if_greater(T& left, U&& right)
    {
        return set_if(left, std::forward<U>(right), std::greater<>{});
    }

    template<typename T, typename U>
    constexpr T& set_if_lesser(T& left, U&& right)
    {
        return set_if(left, std::forward<U>(right), std::less<>{});
    }

    template<typename T, typename Compare>
    [[nodiscard]] constexpr bool is_between(
        const T& v,
        const boost::type_identity_t<T>& min,
        const boost::type_identity_t<T>& max,
        Compare cmp
    ) { return std::addressof(std::clamp(v, min, max, cmp)) == std::addressof(v); }

    template<typename T>
    [[nodiscard]] constexpr bool is_between(
        const T& v,
        const boost::type_identity_t<T>& min,
        const boost::type_identity_t<T>& max
    ) { return is_between(v, min, max, std::less<>{}); }

#define AUTO_MEMBER(name, ini) decltype(ini) name = ini

#define CONST_AUTO_MEMBER(name, ini) const decltype(ini) name = ini

#define CONST_AUTO_REF_MEMBER(name, ini) const decltype(ini)& name = ini

    namespace details
    {
        inline static std::random_device random_device;

        inline static std::mt19937 random_engine{std::random_device{}()};
    }

    constexpr auto& get_random_device() { return details::random_device; }

    constexpr auto& get_random_engine() { return details::random_engine; }

    CPP_template(typename T)(requires std::is_enum_v<T>)
    constexpr auto to_underlying(const T v) { return static_cast<std::underlying_type_t<T>>(v); }
}

ext: h

max_stars_repo_name: BlurringShadow/Polycumic

max_issues_repo_name: BlurringShadow/Polycumic

max_forks_repo_name: BlurringShadow/Polycumic


text:
#ifndef COMMON_HEADER_H
#define COMMON_HEADER_H
//Include Standard Headers
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdarg.h>
#include <stdint.h>
#include <stdbool.h>
//Math Libraries

//#include <math.h>
#include <limits.h>
//#include <gsl/gsl_sf_bessel.h>
//RNG Libraries:
#include <linux/random.h>

//POSIX libraries:
#include <unistd.h>
#include <pthread.h>

//Ncurses
//#include <ncurses.h>
#endif

max_stars_repo_path: Third Edition/Chapter 2/common_header.h

max_issues_repo_path: Third Edition/Chapter 2/common_header.h

max_forks_repo_path: Third Edition/Chapter 2/common_header.h


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <stddef.h>
#include <unistd.h>
#include <math.h>
#include "bigfile-mpi.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void usage() {
    fprintf(stderr, "usage: bigfile-sample-mpi [-r ratio] [-N Nfile] [-f newfilepath] filepath block newblock\n");
    exit(1);

}
#define DONE_TAG 1293
#define ERROR_TAG 1295
#define DIE_TAG 1290
#define WORK_TAG 1291

MPI_Datatype MPI_TYPE_WORK;
BigFile bf = {0};
BigFile bfnew = {0};
BigBlock bb = {0};
BigBlock bbnew = {0};
int verbose = 0;
int Nfile = -1;
size_t CHUNKSIZE = 1 * 1024 * 1024;
int ThisTask, NTask;
char * newfilepath = NULL;
void slave(void);
void server(void);

double ratio = 1.0;
struct work {
    int64_t offset;
    int64_t seed;
    int64_t chunksize;
    int64_t offsetnew;
    int64_t nsel;
};

static size_t filesize();

int main(int argc, char * argv[]) {

    MPI_Init(&argc, &argv);

    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
    MPI_Comm_size(MPI_COMM_WORLD, &NTask);

    MPI_Type_contiguous(sizeof(struct work), MPI_BYTE, &MPI_TYPE_WORK);
    MPI_Type_commit(&MPI_TYPE_WORK);

    int ch;
    while(-1 != (ch = getopt(argc, argv, "n:N:vf:r:"))) {
        switch(ch) {
            case 'r':
                ratio = atof(optarg);
                break;
            case 'N':
            case 'n':
                Nfile = atoi(optarg);
                break;
            case 'f':
                newfilepath = optarg;
                break;
            case 'v':
                verbose = 1;
                break;
            default:
                usage();
        }
    }
    if(argc - optind + 1 != 4) {
        usage();
    }
    argv += optind - 1;
    if(0 != big_file_mpi_open(&bf, argv[1], MPI_COMM_WORLD)) {
        fprintf(stderr, "failed to open: %s\n", big_file_get_error_message());
        exit(1);
    }
    if(0 != big_file_mpi_open_block(&bf, &bb, argv[2], MPI_COMM_WORLD)) {
        fprintf(stderr, "failed to open: %s\n", big_file_get_error_message());
        exit(1);
    }
    if(Nfile == -1 || bb.Nfile == 0) {
        Nfile = bb.Nfile;
    }
    if(newfilepath == NULL) {
        newfilepath = argv[1];
    }
    if(0 != big_file_mpi_create(&bfnew, newfilepath, MPI_COMM_WORLD)) {
        fprintf(stderr, "failed to open: %s\n", big_file_get_error_message());
        exit(1);
    }
    size_t newsize = filesize();
    if(0 != big_file_mpi_create_block(&bfnew, &bbnew, argv[3], bb.dtype, bb.nmemb, Nfile, newsize, MPI_COMM_WORLD)) {
        fprintf(stderr, "failed to create temp: %s\n", big_file_get_error_message());
        exit(1);
    }

    /* copy attrs */
    size_t nattr;
    BigAttr * attrs = big_block_list_attrs(&bb, &nattr);
    int i;
    for(i = 0; i < nattr; i ++) {
        BigAttr * attr = &attrs[i];
        big_block_set_attr(&bbnew, attr->name, attr->data, attr->dtype, attr->nmemb);
    }

    if(bb.nmemb > 0 && bb.size > 0) {
    /* copy data */
        if(ThisTask == 0) {
            server();
        } else {
            slave();
        }
    }
    if(0 != big_block_mpi_close(&bbnew, MPI_COMM_WORLD)) {
        fprintf(stderr, "failed to close new: %s\n", big_file_get_error_message());
        exit(1);
    }
    big_block_mpi_close(&bb, MPI_COMM_WORLD);
    big_file_mpi_close(&bf, MPI_COMM_WORLD);
    big_file_mpi_close(&bfnew, MPI_COMM_WORLD);
    return 0;
}
static size_t filesize() {
    gsl_rng * rng = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rng, 1984);
    int64_t offset = 0;
    int64_t offsetnew = 0;
    struct work work;
    for(offset = 0; offset < bb.size; ) {
        int64_t chunksize = CHUNKSIZE;

        /* never read beyond my end (read_simple caps at EOF) */
        if(offset + chunksize >= bb.size) {
            /* this is the last chunk */
            chunksize = bb.size - offset;
        }
        work.offset = offset;
        work.chunksize = chunksize;
        work.seed = gsl_rng_get(rng);
        work.offsetnew = offsetnew;
        if(ratio == 1.0) {
            work.nsel = chunksize;
        } else {
            work.nsel = gsl_ran_poisson(rng, chunksize * ratio);
        }

        offset += chunksize;
        offsetnew += work.nsel;
    }
    return offsetnew;
}
void server() {
    int64_t offset = 0;
    int64_t offsetnew = 0;
    struct work work;
    gsl_rng * rng = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rng, 1984);
    for(offset = 0; offset < bb.size; ) {
        int64_t chunksize = CHUNKSIZE;
        MPI_Status status;
        int result = 0;
        MPI_Recv(&result, 1, MPI_INT, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD,
                &status);
        if(status.MPI_TAG == ERROR_TAG) {
            break;
        }

        /* never read beyond my end (read_simple caps at EOF) */
        if(offset + chunksize >= bb.size) {
            /* this is the last chunk */
            chunksize = bb.size - offset;
        }
        work.offset = offset;
        work.chunksize = chunksize;
        work.seed = gsl_rng_get(rng);
        work.offsetnew = offsetnew;
        if(ratio == 1.0) {
            work.nsel = chunksize;
        } else {
            work.nsel = gsl_ran_poisson(rng, chunksize * ratio);
        }
        MPI_Send(&work, 1, MPI_TYPE_WORK, status.MPI_SOURCE, WORK_TAG, MPI_COMM_WORLD);

        offset += chunksize;
        offsetnew += work.nsel;
        if(verbose) {
            fprintf(stderr, "%td / %td done (%0.4g%%)\r", offset, bb.size, (100. / bb.size) * offset);
        }
    }
    int i;
    for(i = 1; i < NTask; i ++) {
        struct work work;
        MPI_Send(&work, 1, MPI_TYPE_WORK, i, DIE_TAG, MPI_COMM_WORLD);
    }

}
void slave() {
    gsl_rng * rng = gsl_rng_alloc(gsl_rng_mt19937);
    int result = 0;
    MPI_Send(&result, 1, MPI_INT, 0, DONE_TAG, MPI_COMM_WORLD);
    while(1) {
        struct work work;
        MPI_Status status;
        MPI_Recv(&work, 1, MPI_TYPE_WORK, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);

        if(status.MPI_TAG == DIE_TAG) {
            break;
        }
        gsl_rng_set(rng, work.seed);

        int64_t offset = work.offset;
        int64_t chunksize = work.chunksize;
        int64_t offsetnew = work.offsetnew;
        int64_t nsel = work.nsel;
        BigArray array;
        BigBlockPtr ptrnew;
        BigArray arraynew;

        size_t dims[2];
        void * buffer = malloc(dtype_itemsize(bb.dtype) * bb.nmemb * nsel);
        dims[0] = nsel;
        dims[1] = bb.nmemb;
        big_array_init(&arraynew, buffer, bb.dtype, 2, dims, NULL);

        ptrdiff_t i;
        size_t step = dtype_itemsize(bb.dtype) * bb.nmemb;
        size_t leftover = chunksize;
        char * p = buffer;
        char * q;
        if(0 != big_block_read_simple(&bb, offset, chunksize, &array, NULL)) {
            fprintf(stderr, "failed to read original: %s\n", big_file_get_error_message());
            result = -1;
            goto bad;
        }
        q = array.data;

//        printf("%ld %ld\n", nsel, leftover);
        for(i = 0; i < chunksize; i ++) {
            int64_t r = gsl_rng_uniform_int(rng, leftover);
            if(r < nsel) {
                memcpy(p, q, step);
                p += step;
                nsel --;
            }
            if(nsel == 0) break;
            leftover --;
            q += step;
        }
        if(nsel != 0) abort();
        free(array.data);
        if(0 != big_block_seek(&bbnew, &ptrnew, offsetnew)) {
            fprintf(stderr, "failed to seek new: %s\n", big_file_get_error_message());
            result = -1;
            free(arraynew.data);
            goto bad;
        }

        if(0 != big_block_write(&bbnew, &ptrnew, &arraynew)) {
            fprintf(stderr, "failed to write new: %s\n", big_file_get_error_message());
            result = -1;
            free(arraynew.data);
            goto bad;
        }

        free(arraynew.data);
        MPI_Send(&result, 1, MPI_INT, 0, DONE_TAG, MPI_COMM_WORLD);
        continue;
    bad:
        MPI_Send(&result, 1, MPI_INT, 0, ERROR_TAG, MPI_COMM_WORLD);
        continue;
    }
    return;
}

max_stars_repo_name: rainwoodman/bigfile

max_issues_repo_name: rainwoodman/bigfile

max_forks_repo_name: rainwoodman/bigfile


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

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

#include "utils.h"
#include "norm.h"

/*
 * * Program that reads reference spectra and measure spectra data
 * and linearly unmixes the data using method of least squares fit.
 * Assumes readings is set to 400 and file path is preset to:
 *
 * syntax: ./linearUnmixer MUfilename ref1filename ref2filename ref3filename ...
 *  returns: completion status, prints composition of ref1, ref2, ... , refn
 *
 *  * ****************************************************************
 *  * Author       Dept.           Date            Notes
 *  * ****************************************************************
 *  * Jake Z       Bio. Eng.       Mar 12 2020     Initial version
 *  * Jake Z            "          Apr 03 2020     It works now!
 *  TODO: investigate the discrepancy between what you set $readings and the amount
 *  of space you allocate to the matrices as you forgot to account for the initial
 *  data clearing of the title headers in the data.
 * */

const int readings = 1000;	//num of data points in each reading
const int skip = 50;		//num of data points to discard from beginning

int main(int argc, char *argv[])
{
	if ( argc <= 2 )
	{
		fprintf( stderr, "error, incorrect usage.\n correct usage: ./unmix unknown ref1 ref2 ... refn\n");
		return 3;
	}

	//allocate memory for vars
	int numRefs = argc - 2;
	char *filePath = "./assets/";	// create prompt/arg for desired assets

	gsl_vector *muVector = gsl_vector_alloc(readings);
	gsl_vector *fVector = gsl_vector_alloc(readings);
	gsl_vector *x = gsl_vector_alloc(numRefs);

	gsl_matrix *refMatrix = gsl_matrix_alloc(readings, numRefs);
	gsl_matrix *C = gsl_matrix_alloc( numRefs, numRefs );
	gsl_matrix *Cinverse = gsl_matrix_alloc( numRefs, numRefs );

	int s;
	gsl_permutation *p = gsl_permutation_alloc(numRefs);

	//initialize vars to 0

	gsl_vector_set_zero(muVector);
	gsl_vector_set_zero(fVector);
	gsl_vector_set_zero(x);
	gsl_matrix_set_zero(refMatrix);
	gsl_matrix_set_zero(C);
	gsl_permutation_init(p);	//INITIALISED THE PERMUTATION

	readFileVector( muVector, filePath, argv[1] );	//read vals of unknown arg
	readFileMatrix( refMatrix, filePath, argc, argv );	//read the refs into matrix

	normalizeVector( muVector );

	//create fn that reads the files

	FILE *csvfile = fopen( "./muVector.csv", "wt" );
	gsl_vector_fprintf( csvfile, muVector, "%f" );
	fclose(csvfile);

	FILE *csvfile2 = fopen( "./refVectors.txt", "wt" );
	gsl_matrix_fprintf( csvfile2, refMatrix, "%f" );
	fclose(csvfile2);

	//now we have variable muVector with vals and refMatrix with vals.

	gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1, refMatrix, refMatrix, 0, C);       //matrix C is now A^T A
	gsl_blas_dgemv( CblasTrans, 1, refMatrix, muVector, 0, fVector);        //fVector is 'b'
	gsl_linalg_LU_decomp( C, p, &s);
	gsl_linalg_LU_invert( C, p, Cinverse );

	//TODO maybe implement the other way you can solve for basis using inverse on both sides, (A^T A )^-1 A^T A x = (a^t a)^-1 (a^t b) = x

	gsl_blas_dgemv( CblasNoTrans, 1, Cinverse, fVector, 0, x);

	gsl_vector_fprintf( stdout, x, "%f" );

	gsl_matrix_free(Cinverse);
	gsl_matrix_free(C);
	gsl_matrix_free(refMatrix);
	gsl_vector_free(muVector);
	gsl_vector_free(fVector);
	gsl_vector_free(x);

	return 0;
}

max_stars_repo_name: jakeinater/spectralILU

max_issues_repo_name: jakeinater/spectralILU

max_forks_repo_name: jakeinater/spectralILU


text:
/*
 * kjg_fpca.c
 *
 *  Created on: Apr 28, 2014
 *      Author: Kevin
 */

#include <stdio.h>

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

#include "kjg_fpca.h"
#include "kjg_gsl.h"

#include "admutils.h"
#include "gval.h"

size_t KJG_FPCA_ROWS = 256;

void
kjg_fpca (size_t K, size_t L, size_t I, double *eval, double *evec)
{
  if (K >= L)
    exit (1);
  if (I == 0)
    exit (1);

  size_t m = get_ncols ();
  size_t n = get_nrows ();

  // PART A - compute Q such that X ~ Q * (Q^T) * X
  gsl_matrix *G1 = gsl_matrix_alloc (n, L);
  gsl_matrix *G2 = gsl_matrix_alloc (n, L);
  gsl_matrix *Q = gsl_matrix_alloc (m, (I + 1) * L);
  gsl_matrix *Gswap;

  gsl_rng *r = kjg_gsl_rng_init ();
  kjg_gsl_ran_ugaussian_matrix (r, G1);
  gsl_rng_free (r);

  size_t i;
  for (i = 0; i < I; i++)
    {
      gsl_matrix_view Qi = gsl_matrix_submatrix (Q, 0, i * L, m, L);

      // do the multiplication
      kjg_fpca_XTXA (G1, &Qi.matrix, G2);

      // scale to prevent G2 from blowing up
      gsl_matrix_scale (G2, 1.0 / m);

      Gswap = G2;
      G2 = G1;
      G1 = Gswap;
    }

  gsl_matrix_view Qi = gsl_matrix_submatrix (Q, 0, I * L, m, L);
  kjg_fpca_XA (G1, &Qi.matrix);

  {
    gsl_matrix *V = gsl_matrix_alloc (Q->size2, Q->size2);
    gsl_vector *S = gsl_vector_alloc (Q->size2);

    kjg_gsl_SVD (Q, V, S);

    gsl_matrix_free (V);
    gsl_vector_free (S);
  }

  // kjg_gsl_matrix_QR(Q); // QR decomposition is less accurate than SVD

  gsl_matrix_free (G1);
  gsl_matrix_free (G2);

  // PART B - compute B matrix, take SVD and return
  gsl_matrix *B = gsl_matrix_alloc (n, (I + 1) * L);
  kjg_fpca_XTB (Q, B);

  gsl_matrix *Utilda = gsl_matrix_alloc ((I + 1) * L, (I + 1) * L);
  gsl_vector *Stilda = gsl_vector_alloc ((I + 1) * L);

  kjg_gsl_SVD (B, Utilda, Stilda);

  gsl_matrix_view Vk = gsl_matrix_submatrix (B, 0, 0, n, K);
  gsl_matrix_view evec_view = gsl_matrix_view_array (evec, n, K);
  gsl_matrix_memcpy (&evec_view.matrix, &Vk.matrix);

  gsl_vector_view Sk = gsl_vector_subvector (Stilda, 0, K);
  gsl_vector_view eval_view = gsl_vector_view_array (eval, K);
  gsl_vector_mul (&Sk.vector, &Sk.vector);
  gsl_vector_scale (&Sk.vector, 1.0 / m);
  gsl_vector_memcpy (&eval_view.vector, &Sk.vector);

  gsl_matrix_free (Q);
  gsl_matrix_free (B);
  gsl_matrix_free (Utilda);
  gsl_vector_free (Stilda);
}

void
kjg_fpca_XTXA (const gsl_matrix * A1, gsl_matrix * B, gsl_matrix * A2)
{
  size_t m = get_ncols ();
  size_t n = get_nrows ();

  size_t i, r;			// row index
  double *Y = malloc (sizeof (double) * n * KJG_FPCA_ROWS);	// normalized

  gsl_matrix_view Bi, Xi;

  gsl_matrix_set_zero (A2);

  for (i = 0; i < m; i += KJG_FPCA_ROWS)
    {
      r = kjg_geno_get_normalized_rows (i, KJG_FPCA_ROWS, Y);
      Xi = gsl_matrix_view_array (Y, r, n);
      Bi = gsl_matrix_submatrix (B, i, 0, r, B->size2);
      gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1, &Xi.matrix, A1, 0,
		      &Bi.matrix);
      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1, &Xi.matrix, &Bi.matrix,
		      1, A2);
    }

  free (Y);
}

void
kjg_fpca_XA (const gsl_matrix * A, gsl_matrix * B)
{
  size_t n = get_nrows ();
  size_t m = get_ncols ();

  size_t i, r;
  double *Y = malloc (sizeof (double) * n * KJG_FPCA_ROWS);

  gsl_matrix_view Hmat, Xmat;

  gsl_matrix_set_zero (B);

  for (i = 0; i < m; i += KJG_FPCA_ROWS)
    {
      r = kjg_geno_get_normalized_rows (i, KJG_FPCA_ROWS, Y);
      Xmat = gsl_matrix_view_array (Y, r, n);
      Hmat = gsl_matrix_submatrix (B, i, 0, r, B->size2);
      gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1, &Xmat.matrix, A, 0,
		      &Hmat.matrix);
    }

  free (Y);
}

void
kjg_fpca_XTB (const gsl_matrix * B, gsl_matrix * A)
{
  size_t n = get_nrows ();
  size_t m = get_ncols ();

  size_t i, r;
  double *Y = malloc (sizeof (double) * n * KJG_FPCA_ROWS);
  gsl_matrix_view Xmat;

  gsl_matrix_set_zero (A);

  for (i = 0; i < m; i += KJG_FPCA_ROWS)
    {
      r = kjg_geno_get_normalized_rows (i, KJG_FPCA_ROWS, Y);
      Xmat = gsl_matrix_view_array (Y, r, n);
      gsl_matrix_const_view Hmat = gsl_matrix_const_submatrix (B, i, 0, r,
							       B->size2);
      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1, &Xmat.matrix, &Hmat.matrix,
		      1, A);
    }

  free (Y);
}

max_stars_repo_name: nevrome/EIG

max_issues_repo_name: nevrome/EIG

max_forks_repo_name: nevrome/EIG


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

#ifndef __GSL_SPLINE2D_H__
#define __GSL_SPLINE2D_H__

#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

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

__BEGIN_DECLS

/*
 * A 2D interpolation object which stores the arrays defining the function.
 * In all other respects, this is just like a gsl_interp2d object.
 */
typedef struct
{
  gsl_interp2d interp_object; /* low-level interpolation object */
  double * xarr;              /* x data array */
  double * yarr;              /* y data array */
  double * zarr;              /* z data array */
} gsl_spline2d;

gsl_spline2d * gsl_spline2d_alloc(const gsl_interp2d_type * T, size_t xsize, size_t ysize);

int gsl_spline2d_init(gsl_spline2d * interp, const double xa[],
                      const double ya[], const double za[],
                      size_t xsize, size_t ysize);

void gsl_spline2d_free(gsl_spline2d * interp);

double gsl_spline2d_eval(const gsl_spline2d * interp, const double x,
                         const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

int gsl_spline2d_eval_e(const gsl_spline2d * interp, const double x,
                        const double y, gsl_interp_accel* xa, gsl_interp_accel* ya,
                        double * z);

double gsl_spline2d_eval_deriv_x(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

int gsl_spline2d_eval_deriv_x_e(const gsl_spline2d * interp, const double x,
                                const double y, gsl_interp_accel* xa,
                                gsl_interp_accel* ya, double * z);

double gsl_spline2d_eval_deriv_y(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya);

int gsl_spline2d_eval_deriv_y_e(const gsl_spline2d * interp, const double x,
                                const double y, gsl_interp_accel* xa,
                                gsl_interp_accel* ya, double * z);

double gsl_spline2d_eval_deriv_xx(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

int gsl_spline2d_eval_deriv_xx_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

double gsl_spline2d_eval_deriv_yy(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

int gsl_spline2d_eval_deriv_yy_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

double gsl_spline2d_eval_deriv_xy(const gsl_spline2d * interp, const double x,
                                  const double y, gsl_interp_accel* xa, gsl_interp_accel* ya);

int gsl_spline2d_eval_deriv_xy_e(const gsl_spline2d * interp, const double x,
                                 const double y, gsl_interp_accel* xa,
                                 gsl_interp_accel* ya, double * z);

size_t gsl_spline2d_min_size(const gsl_spline2d * interp);

const char * gsl_spline2d_name(const gsl_spline2d * interp);

int gsl_spline2d_set(const gsl_spline2d * interp, double zarr[],
                     const size_t i, const size_t j, const double z);
double gsl_spline2d_get(const gsl_spline2d * interp, const double zarr[],
                        const size_t i, const size_t j);

__END_DECLS

#endif /* __GSL_SPLINE2D_H__ */

max_forks_repo_name: ielomariala/Hex-Game


text:
#include <stdio.h>
#include <stdlib.h>
#include <argp.h>
#include <string.h>
#include <assert.h>
#include <stdint.h>
#include "utils.h"

#ifdef USE_MKL
#pragma message "Using the MKL for BLAS."
#include <mkl.h>
#else
#pragma message "Not using the MKL for BLAS."
#include <cblas.h>
#endif

#define MAX_LINES 2000
#define NB_RUNS 1
#define MAX_NAME_SIZE 256

#define MAX_OFFSET 1000000
#define WRITE_MEMORY_SIZE 100000000

FILE* active_file = NULL;
char* basename = "calibration";
char* dir_name  = ".";

static double *matrix_A;
static double *matrix_B;
static double *matrix_C;

unsigned long long base_time;

typedef struct {
  char* sizefile;
  int loop;
  char* resultfile;
} my_args;

static char doc[] = "Runs BLAS benchmarks to compute values used for SMPI calibration";

static struct argp_option options[] = {
  {"sizeFile", 's', "SIZEFILE", 0, "filename of the size list"},
  {"resultfile", 'o', "RESULTFILE", 0, "filename of the results"},
  {"loop", 'l', "NB_LOOPS", 0, "number of (non-recorded) loops to run before and after the execution"},
  { 0 }
};

static int parse_options (int key, char *arg, struct argp_state *state)
{
  my_args *arguments = state->input;
  switch (key){
  case 'l':
    arguments->loop = atoi(arg);
    break;
  case 's':
    arguments->sizefile = arg;
    break;
  case 'o':
    arguments->resultfile = arg;
    break;

  default: return ARGP_ERR_UNKNOWN;
  }
  return 0;
}

struct argp argp = { options, parse_options, NULL, doc };
my_args arguments;

/*
 * Return 0b000...01...1, with n ones at the end.
 */
uint64_t __get_mask(unsigned n) {
    if(n == 0)
        return 0;
    assert(n >= 0 && n < 64);
    uint64_t mask = 1;
    return (mask << n) - 1;
}

/*
 * Return 0b000...01...10...0.
 * The bits at positions [start, stop] are equal to 1, the others to 0 (start is the lower order).
 */
uint64_t get_mask(unsigned start, unsigned stop) {
    assert(0 <= start && start <= stop && stop < 64);
    uint64_t ones = __get_mask(stop);
    uint64_t zeroes = ~__get_mask(start);
    return ones & zeroes;
}

/*
 * Apply the given mask to a double.
 */
double apply_mask(double x, uint64_t mask) {
    uint64_t *tmp = (uint64_t*)&x;
    (*tmp) |= mask;
    return *((double*)tmp);
}

void __print_bits(uint64_t n, unsigned i) {
    if(i == 64)
        return;
    __print_bits(n / 2, i+1);
    if(n % 2)
        printf("1");
    else
        printf("0");
}

void print_bits(uint64_t n) {
    printf("0b");
    __print_bits(n, 0);
    printf("\n");
}

void print_bits_f(double x) {
    uint64_t *tmp = (uint64_t*)&x;
    print_bits(*tmp);
}

double *allocate_matrix(size_t size1, size_t size2) {
    size_t alloc_size = (size1*size2 + MAX_OFFSET) * sizeof(double);
    double *result = (double*) malloc(alloc_size);
    if(!result) {
      perror("malloc");
      exit(errno);
    }
#ifdef MASK_SIZE
    uint64_t mask = get_mask(0, MASK_SIZE);
#pragma message "Using a mask for the values of the matrix."
#else
    // Note that this particular mask does not do anything: apply_mask(x, get_mask(0, 0)) == x
    // One need to change the stop and/or start indices to really apply a mask.
    uint64_t mask = get_mask(0, 0);
#endif
    printf("Using the mask: ");
    print_bits(mask);
    double x = 3.1415926535897932384626433;
    assert(x == apply_mask(x, get_mask(0, 0)));
    assert(x != apply_mask(x, get_mask(0, 1)));
    for(int i = 0; i < size1*size2; i++) {
        result[i] = apply_mask((double)rand()/(double)(RAND_MAX), mask);
    }
    return result;
}

void write_memory(void) {
    static char *buff = NULL;
    if(!buff) {
        buff = malloc(WRITE_MEMORY_SIZE);
        if(!buff) {
            perror("malloc");
            exit(errno);
        }
    }
    memset(buff, rand()%256, WRITE_MEMORY_SIZE);
}

FILE *open_file(const char* filename){
    FILE *file = fopen(filename, "w");
    if(!file) {
        perror("open_file");
        exit(errno);
    }
    return file;
}

void get_dgemm(FILE *file, int *sizes, int nb_it, unsigned long long base_time, int write_file) {
  unsigned long long start_time, total_time;
  double alpha = 1., beta=1.;
  for(int i=0; i<nb_it; i++) {
    start_time=get_time();
    size_t offset = rand() % MAX_OFFSET;
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, sizes[0], sizes[1], sizes[2], alpha, matrix_A+offset, sizes[3],
            matrix_B+offset, sizes[4], beta, matrix_C+offset, sizes[5]);
    total_time=get_time()-start_time;
    if(write_file)
      print_in_file(file, "dgemm", sizes, 6, start_time-base_time, total_time);
  }
}

static const char *names[] = {"dgemm", NULL};
static const void (*functions[])(FILE*, int*, int, unsigned long long, int) = {get_dgemm};

void test_op(FILE *result_file, experiment_t *exp, int nb_runs, unsigned long long base_time, int write_file) {
  functions[exp->op_id](result_file, exp->sizes, nb_runs, base_time, write_file);
}

int get_max_size(experiment_t *exp, int nb_exp, int idx) {
  int max = -1;
  for(int i = 0; i < nb_exp; i++) {
    int size = exp[i].sizes[idx];
    if(max < size)
      max = size;
  }
  return max;
}

int main(int argc, char** argv){

  srand(42);
  bzero (&arguments, sizeof(my_args));

  if (argp_parse (&argp, argc, argv, 0, 0, &arguments) == ARGP_KEY_ERROR){
    fprintf(stderr,"error during the parsing of parameters\n");
    return 1;
  }

  int nb_loops = arguments.loop;

  if(arguments.sizefile==NULL){
    fprintf(stderr, "Please provide a name for a file containing a list of sizes\n");
    return -1;
  }
  if(arguments.resultfile==NULL){
    fprintf(stderr, "Please provide a name for a result file\n");
    return -1;
  }

  FILE *result_file = NULL;
  if(arguments.resultfile)
    result_file = open_file(arguments.resultfile);

  int nb_exp, largest_size;

  experiment_t *experiments = parse_experiment_file(names, arguments.sizefile, &nb_exp, &largest_size, -1, 300000, 6);

  printf("nb_exp=%d, largest_size=%d\n", nb_exp, largest_size);

  size_t max_size = (size_t)largest_size*(size_t)largest_size;
  int max_M = get_max_size(experiments, nb_exp, 0);
  int max_N = get_max_size(experiments, nb_exp, 1);
  int max_K = get_max_size(experiments, nb_exp, 2);
  printf("max(M)=%d | max(N)=%d | max(K)=%d\n", max_M, max_N, max_K);
  printf("Alloc size: was %.2e bytes\n", (double)(max_size)*sizeof(double)*3);
  double alloc_size = (double)max_M*(double)max_K;
  alloc_size += (double)max_K*(double)max_N;
  alloc_size += (double)max_M*(double)max_N;
  printf("Alloc size: now %.2e bytes\n", alloc_size*sizeof(double));

  matrix_A = allocate_matrix(max_M, max_K);
  matrix_B = allocate_matrix(max_K, max_N);
  matrix_C = allocate_matrix(max_M, max_N);

  //Init time
  base_time=get_time();
  for(int i = 0; i < nb_loops; i++) {
    for(int j = 0; j < nb_exp; j++) {
      test_op(result_file, &experiments[j], NB_RUNS, base_time, 0);
    }
  }
  for(int j = 0; j < nb_exp; j++) {
    test_op(result_file, &experiments[j], NB_RUNS, base_time, 1);
  }
  for(int i = 0; i < nb_loops; i++) {
    for(int j = 0; j < nb_exp; j++) {
      test_op(result_file, &experiments[j], NB_RUNS, base_time, 0);
    }
  }

  fflush(result_file);
  fclose(result_file);
  free(experiments);
  free(matrix_A);
  free(matrix_B);
  free(matrix_C);
  return 0;
}

max_stars_repo_name: Ezibenroc/platform-calibration

max_issues_repo_name: Ezibenroc/platform-calibration

max_forks_repo_name: Ezibenroc/platform-calibration


text:
#if !defined(UTIL_H)
#define UTIL_H

#include <gsl/gsl_matrix.h>

// Flag to set if the program should exit on GSL error
const static bool EXIT_ON_GSL_ERROR = false;

struct stability_params {
  double lambda[8];
  double l;
  double u;
  double m;
  double v;
};

struct integration_params {
  double M;
  double m1;
  double m2;
  double gamma;
};

struct integration_params2 {
  double M;
  double x1;
  double m;
  double gamma;
};

struct integration_params_tt {
  double M;
  double mW;
  double mb;
  double mt;
  double x1;
  double gtop;
  int h;
};

struct integration_params_tb {
  double M;
  double mW;
  double mb;
  double mt;
  double x1;
  double gtop;
  double md;
  double Z_u;
  double Z_d;
};

int sign(double x);

double stability_fcn(double x, void *params);

bool stability_minimum(stability_params &p);

void print_gsl_matrix(gsl_matrix *mat,int m,int n);

double fint(double x, void* par);

double gint(double x, void* par);

double Lint_s(double x, void *param);

double Lint_p(double x, void *param);

double Lint_c(double x, void *param);

double hvv_fcn(double y, void *params);

double hvv_fcn1(double y, void *params);

double hvv_fcn2(double y, void *params);

double htt_fcn1(double y, void *params);

double htt_fcn2(double y, void *params);

double htb_fcn1(double y, void *params);

double htb_fcn2(double y, void *params);

double cubic(double x, double X[4], double Y[4]);

double hvh_fcn(double x, void *params);

double  L(double x, double y, double z);
double  DHp(double ui, double uj, double xi, double xj, double sqL);
double  DHm(double ui, double uj, double xi, double xj, double sqL);
double  BHp(double ui, double uj, double xi, double xj, double sqL);
double sqrtlambda(double a1, double a2, double a3);

#endif


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

#pragma once

#ifndef rectpack_89b58399_0f97_4cb9_9160_0226504762b1_h
#define rectpack_89b58399_0f97_4cb9_9160_0226504762b1_h

#include <ariel/config.h>
#include <gslib/bintree.h>
#include <gslib/std.h>
#include <gslib/math.h>
#include <ariel/type.h>

__ariel_begin__

struct rp_rect
{
    float               x, y;                           /* left, top */
    float               width, height;

public:
    rp_rect() { x = y = width = height = 0.f; }
    rp_rect(float l, float t, float w, float h) { set_rect(l, t, w, h); }
    void set_rect(float l, float t, float w, float h) { x = l, y = t, width = w, height = h; }
    float left() const { return x; }
    float top() const { return y; }
    float right() const { return x + width; }
    float bottom() const { return y + height; }
    void move_to(float destx, float desty) { x = destx, y = desty; }
    void offset(float cx, float cy) { x += cx, y += cy; }
    float area() const { return width * height; }
};

struct rp_node:
    public rp_rect
{
    bool                transposed;
    void*               bind_ptr;

public:
    rp_node()
    {
        transposed = false;
        bind_ptr = nullptr;
    }
    virtual ~rp_node() {}
    virtual bool is_dynamic() const { return false; }
    virtual void set_transposed(bool b) { transposed = b; }
    virtual string to_string() const { return string(); }   /* for debug only */
    void set_bind_ptr(void* p) { bind_ptr = p; }
    bool is_empty() const { return !bind_ptr; }
};

struct rp_dynamic_node:
    public rp_node
{
    rp_rect             blank;
    bool                has_blank;                      /* tag of blank.                | only valid on a non-leaf node. */
    bool                is_vert_arrange;                /* vertical arrangement?        | only valid on a non-leaf node. */

public:
    rp_dynamic_node() { has_blank = false; }
    virtual bool is_dynamic() const override { return true; }
    virtual void set_transposed(bool b) override;
    virtual string to_string() const override;
    float get_area() const;
    void set_no_blank();
    void set_blank(const rp_rect& rc);
    void set_from_local(rp_dynamic_node& p);
    void set_dimensions(float w, float h);
    void move_to_pos(float x, float y);
    void flip_leaf();
    void flip_non_leaf();
};

typedef _bintreenode_wrapper<rp_node> rp_wrapper;
typedef _bintree_allocator<rp_wrapper> rp_allocator;
typedef bintree<rp_node, rp_wrapper, rp_allocator> rp_tree;
typedef rp_tree::iterator rp_iterator;
typedef rp_tree::const_iterator rp_const_iterator;

struct rp_input
{
    float               width, height;
    void*               binding;
};

typedef list<rp_input> rp_input_list;

extern void rp_global_location(rp_const_iterator i, rp_rect& rc);
extern bool rp_is_node_transposed(rp_const_iterator i);

class rect_packer
{
public:
    enum packing_strategy
    {
        ps_unknown,
        ps_compactly,
        ps_dynamically,
    };

public:
    rect_packer();
    ~rect_packer();
    const rp_node& get_root_node() const { return *_tree.const_root(); }
    bool is_empty() const { return !_tree.is_valid(); }
    float get_width() const { return _tree.const_root()->width; }
    float get_height() const { return _tree.const_root()->height; }
    int get_pack_times() const { return _pack_times; }
    void pack_automatically(rp_input_list& inputs);
    void tracing() const;
    void trace_total() const;
    void trace_blank() const;
    void trace_tree() const;
    bool checkups() const;

public:
    template<class _lamb>
    void for_each(rp_iterator i, _lamb fn)
    {
        assert(i);
        _tree.preorder_traversal([&fn](rp_wrapper* w) {
            assert(w);
            rp_iterator i(w);
            if(i.is_leaf() && !i->is_empty()) {
                rp_rect rc;
                rp_global_location(i, rc);
                fn(i->bind_ptr, rc, i->transposed);
            }
        });
    }
    template<class _lamb>
    void for_each(_lamb fn)
    {
        if(!_tree.is_valid())
            return;
        for_each(_tree.get_root(), fn);
    }

protected:
    rp_tree             _tree;
    packing_strategy    _strategy;
    int                 _pack_times;

protected:
    void initialize_compactly(float w, float h);
    bool pack_compactly(const rp_input_list& inputs);
    bool add_rect_compactly(float w, float h, void* binding);
    void trace_compactly() const;
    void initialize_dynamically();
    void add_rect_dynamically(float w, float h, void* binding);
    void add_rect_dynamically(rp_tree& rc);
    bool add_rect_dynamically(rp_tree& rc, rp_tree& xchg);
    void proc_fill_up(rp_iterator p, rp_tree& t);
    void proc_combine(rp_iterator p, rp_tree& t);
    void proc_combine_root(rp_iterator p, rp_tree& t);
    void proc_combine_non_root(rp_iterator p, rp_iterator i, rp_tree& t);
    void proc_replace(rp_iterator p, rp_tree& t, rp_tree& xchg);
    void trace_dynamically() const;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
// matrix/snowboy-blas.h

// Copyright 2016  KITT.AI (author: Guoguo Chen)

#ifndef SNOWBOY_MATRIX_SNOWBOY_BLAS_H_
#define SNOWBOY_MATRIX_SNOWBOY_BLAS_H_

#ifdef HAVE_ATLAS
  extern "C" {
    #include <cblas.h>
    #include <clapack.h>
  }
#elif defined(HAVE_CLAPACK)
  #ifdef __APPLE__
    #ifndef __has_extension
    #define __has_extension(x) 0
    #endif
    #define vImage_Utilities_h
    #define vImage_CVUtilities_h
    #include <Accelerate/Accelerate.h>
    typedef __CLPK_integer          integer;
    typedef __CLPK_logical          logical;
//    typedef __CLPK_real             real;
    typedef __CLPK_doublereal       doublereal;
    typedef __CLPK_complex          complex;
    typedef __CLPK_doublecomplex    doublecomplex;
    typedef __CLPK_ftnlen           ftnlen;
  #else
    extern "C" {
      #include <cblas.h>
      #include <f2c.h>
      #include <clapack.h>

      // Removes dangerous macros from f2c.h.
      #undef abs
      #undef dabs
      #undef min
      #undef max
      #undef dmin
      #undef dmax
      #undef bit_test
      #undef bit_clear
      #undef bit_set
    }
  #endif
  // Defines SnowboyBlasInt type that will be used in SVD.
  typedef integer SnowboyBlasInt;
#elif defined(HAVE_OPENBLAS)
  #include "cblas.h"
#ifndef NO_FORTRAN
  #include "lapacke.h"
#endif
  #undef I
  #undef complex
  // get rid of macros from f2c.h -- these are dangerous.
  #undef abs
  #undef dabs
  #undef min
  #undef max
  #undef dmin
  #undef dmax
  #undef bit_test
  #undef bit_clear
  #undef bit_set
#else
  #error "You have to define HAVE_CLAPACK or HAVE_ATLAS or HAVE_OPENBLAS."
#endif

#endif  // SNOWBOY_MATRIX_SNOWBOY_BLAS_H_

max_stars_repo_name: songmeixu/bit-matrix

max_issues_repo_name: songmeixu/bit-matrix

max_forks_repo_name: songmeixu/bit-matrix