text:
/*
**
** Lit:
**   Y.Yu, S.T. Acton:
**   Speckle reducing anisotropic diffusion.
**   IEEE Trans Image Proc., Vol. 11, No. 11, Nov 2002
**
**  G.Lohmann, MPI-CBS, 2008
*/

#include <viaio/Vlib.h>
#include <viaio/VImage.h>

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

#include <gsl/gsl_math.h>

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

/*
** repair border
*/
void
Border(VImage src)
{
  int b=0,r,c,nrows,ncols;
  double u;

  nrows = VImageNRows(src);
  ncols = VImageNColumns(src);

  for (r=0; r<VImageNRows(src); r++) {
    u = VGetPixel(src,b,r,1);
    VSetPixel(src,b,r,0,u);
    u = VGetPixel(src,b,r,ncols-2);
    VSetPixel(src,b,r,ncols-1,u);
  }

  for (c=0; c<VImageNColumns(src); c++) {
    u = VGetPixel(src,b,1,c);
    VSetPixel(src,b,0,c,u);
    u = VGetPixel(src,b,nrows-2,c);
    VSetPixel(src,b,nrows-1,c,u);
  }
}

VImage
VSRAD(VImage src,VImage dest,VShort numiter,VShort type,VFloat rho)
{
  static VImage tmp1=NULL,tmp2=NULL,ctmp=NULL;
  int nbands,nrows,ncols;
  int b,r,c,iter;
  double u,v,w,x,y;
  double d,da,dd,dt,t,q,q0;
  double cx,c0,c1,c2,c3;
  double tiny=1.0e-6;

  nbands = VImageNBands(src);
  nrows  = VImageNRows(src);
  ncols  = VImageNColumns(src);

  if (tmp1 == NULL) {
    tmp1 = VCreateImage(1,nrows,ncols,VFloatRepn);
    tmp2 = VCreateImage(1,nrows,ncols,VFloatRepn);
    ctmp = VCreateImage(1,nrows,ncols,VFloatRepn);
  }

  dt = 0.05;

  dest = VCopyImage(src,dest,VAllBands);

  for (b=0; b<nbands; b++) {

    t = 0;
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	VPixel(tmp1,0,r,c,VFloat) = VGetPixel(src,b,r,c);
      }
    }

    VFillImage(ctmp,VAllBands,0);
    VFillImage(tmp2,VAllBands,0);

    q0 = 1;
    for (iter=0; iter < numiter; iter++) {

      for (r=1; r<nrows-1; r++) {
	for (c=1; c<ncols-1; c++) {

	  u  = VPixel(tmp1,0,r,c,VFloat);
	  if (ABS(u) < tiny) continue;

	  v  = VPixel(tmp1,0,r+1,c,VFloat);
	  w  = VPixel(tmp1,0,r,c+1,VFloat);
	  x  = VPixel(tmp1,0,r-1,c,VFloat);
	  y  = VPixel(tmp1,0,r,c-1,VFloat);

	  da = (SQR(u-v) + SQR(u-w) + SQR(u-x) + SQR(u-y))/(u*u);
	  dd = (v+w+x+y - 4.0*u)/(u*u);
	  q = (0.5*da - 0.0625*dd) / (SQR(1.0 + 0.25*dd));
	  if (q > 0) {
	    q = sqrt(q);
	    u = (q*q - q0*q0) / ((q0*q0)*(1.0 + q0*q0));
	  }
	  else
	    u = 0;

	  cx = 0;
	  if (type == 0)
	    cx = 1.0 / (1.0 + u);
	  else
	    cx = exp(-u);

	  if (gsl_isnan(cx) || gsl_isinf(cx)) {
	    cx = 0;
	    VWarning(" cx, insnan, isinf");
	  }

	  VPixel(ctmp,0,r,c,VFloat) = cx;
	}
      }
      Border(ctmp);

      for (r=1; r<nrows-1; r++) {
	for (c=1; c<ncols-1; c++) {

	  u  = VPixel(tmp1,0,r,c,VFloat);
	  if (ABS(u) < tiny) continue;

	  v  = VPixel(tmp1,0,r+1,c,VFloat);
	  w  = VPixel(tmp1,0,r,c+1,VFloat);
	  x  = VPixel(tmp1,0,r-1,c,VFloat);
	  y  = VPixel(tmp1,0,r,c-1,VFloat);

	  c0  = VPixel(ctmp,0,r+1,c,VFloat);
	  c1  = VPixel(ctmp,0,r,c+1,VFloat);
	  c2  = VPixel(ctmp,0,r-1,c,VFloat);
	  c3  = VPixel(ctmp,0,r,c-1,VFloat);

          /*  ???
          cx =  VPixel(ctmp,0,r,c,VFloat);
	  d = c0*(v-u) + c1*(w-u) + cx*(x-u) + cx*(y-u);
          */

          d = c0*(v-u) + c1*(w-u) + c2*(x-u) + c3*(y-u);

	  u = u + 0.25*dt * d;
	  if (gsl_isnan(u) || gsl_isinf(u)) {
	    u = 0;
	    VWarning(" insnan, isinf");
	  }

	  if (u > VPixelMaxValue (tmp2)) u = VPixelMaxValue (tmp2);
	  if (u < VPixelMinValue (tmp2)) u = VPixelMinValue (tmp2);
	  VPixel(tmp2,0,r,c,VFloat) = u;

	}
      }
      q0 *= exp(-rho*t);

      Border(tmp2);
      tmp1 = VCopyImagePixels(tmp2,tmp1,VAllBands);
      t += dt;
    }

    /*
    ** output
    */
    for (r=1; r<nrows-1; r++) {
      for (c=1; c<ncols-1; c++) {
	v = VPixel(tmp2,0,r,c,VFloat);
	if (gsl_isnan(v) || gsl_isinf(v)) {
	  v = 0;
	  VWarning(" insnan, isinf");
	}

	if (v > VPixelMaxValue (dest)) v = VPixelMaxValue (dest);
	if (v < VPixelMinValue (dest)) v = VPixelMinValue (dest);
	VSetPixel(dest,b,r,c,(VDouble) v);
      }
    }
  }

  /*
  VDestroyImage(tmp1);
  VDestroyImage(tmp2);
  */
  return dest;
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include "matrix.h"

void Matrix_Multiply(Matrix a, Matrix b, Matrix c) {
    /* for dimentions:
     * a[m*k],b[k*n],c[m*n]
     * We use:
     * RowMajor, NoTrans, NoTrans, m, n, k, 1, A, k, B, n, 0, C, n);
     */

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, a.rows, b.cols, a.cols,
                1, a.values, a.cols, b.values, b.cols, 0, c.values, b.cols);
}

void Matrix_Print(Matrix m) {
    printf("=====\n");
    printf("rows: %lld, columns: %lld\n", m.rows, m.cols);
    for (int i = 0; i < m.rows; i++) {
        for (int j = 0; j < m.cols; j++) {
            printf("%lf ", m.values[i * m.cols + j]);
        }
        printf("\n");
    }
}

void Matrix_Add(Matrix a, Matrix b, Matrix c) {
    for (int i = 0; i < a.rows * a.cols; i++) {
        c.values[i] = a.values[i] + b.values[i];
    }
}

void Matrix_Scale(Matrix a, double n) {
    cblas_dscal(a.rows * a.cols, n, a.values, 1);
}

void Matrix_Test(void) {
    double *A, *B, *C;
    int m, n, k, i;
    double alpha, beta;
    m = 2000, k = 200, n = 1000;
    printf(
            " Initializing data for matrix multiplication C=A*B for matrix \n"
                    " A(%ix%i) and matrix B(%ix%i)\n",
            m, k, k, n);
    alpha = 1.0;
    beta = 0.0;

    A = (double *) malloc(m * k * sizeof(double));
    B = (double *) malloc(k * n * sizeof(double));
    C = (double *) malloc(m * n * sizeof(double));
    if (A == NULL || B == NULL || C == NULL) {
        printf("\n ERROR: Can't allocate memory for matrices. Aborting... \n\n");
        free(A);
        free(B);
        free(C);
        return;
    }

    printf(" Intializing matrix data \n");
    for (i = 0; i < (m * k); i++) {
        A[i] = (double) (i + 1);
    }

    for (i = 0; i < (k * n); i++) {
        B[i] = (double) (-i - 1);
    }

    for (i = 0; i < (m * n); i++) {
        C[i] = 0.0;
    }

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k,
                B, n, beta, C, n);
    printf("\n Computations completed.\n");
    free(A);
    free(B);
    free(C);
}

max_stars_repo_name: RedisLabsModules/redis-ml

max_issues_repo_name: rmikio/redis-ml

max_forks_repo_name: rmikio/redis-ml


text:
/*
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include "spc_FITScards.h"
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_PET.h"
#include "inout_aper.h"
#include "trace_conf.h"
#include "aper_conf.h"
#include "disp_conf.h"
#include "drizzle_utils.h"
#include "crossdisp_utils.h"
#include "spc_cfg.h"
#include "inima_utils.h"
#include "spce_output.h"

#define AXE_IMAGE_PATH "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

int
main (int argc, char *argv[])
{
  char *opt;
  char aper_file[MAXCHAR];
  char aper_file_path[MAXCHAR];

  char list_file[MAXCHAR];
  char list_file_path[MAXCHAR];

  char conf_file[MAXCHAR];
  char conf_file_path[MAXCHAR];

  char grism_file[MAXCHAR];
  char grism_file_path[MAXCHAR];

  char PET_file[MAXCHAR];
  char PET_file_path[MAXCHAR];

  char SEC_file[MAXCHAR];
  char SEC_file_path[MAXCHAR];

  char outputroot[MAXCHAR];
  char outputroot_path[MAXCHAR];

  char output_path[MAXCHAR];

  aperture_conf *conf;

  object **oblist;
  observation *obs;

  ap_pixel *pri_PET;
  ap_pixel *sec_PET;

  d_point     pixel;
  d_point     refwave_pos;
  d_point     outref;
  d_point     minxy_PET;
  px_point    pixmax;
  dispstruct  *disp, *outdisp;
  trace_func  *trace;
  gsl_matrix  *drizzcoeffs;
  //gsl_vector  *lambdaref = gsl_vector_alloc(2);

  FILE *in_file;

  objectobs  **allobjects;

  int nobjects, nobjects2;
  int boxwidth, boxheight,trlength;
  double relx, rely, objwidth, orient;
  //double min_px, min_py;
  double drizzle_width;
  double cdref, cdscale, cdcorr, cdmeanscale;
  double sprefreso, spreso, spmeanreso, spcorr;
  double sky_cps;
  double exptime;

  double *gaga;

  int index, i, for_grism, in_index;
  char label[MAXCHAR];

  fitsfile *PET_ptr, *SEC_ptr, *DPP_ptr;
  int f_status = 0;

  int aperID, beamID, objindex;
  FITScards *cards;

  drzstamp_dim dimension;
  drzprep *drzprep_stamps;

  axe_inputs *input_list;

  //int rectified = 0;
  //int drizzled  = 0;
  int bckmode   = 0;
  int backpet   = 1;
  //int usemode   = 0;
  int quant_cont= 0;
  int opt_extr  = 0;

  if ((argc < 3) || (opt = get_online_option ("help", argc, argv)))
    {
      fprintf (stdout,
	       "aXe_DRZPREP Version %s:\n"
	       "           aXe task to produce a set of Drizzle PrePare (DPP) files\n"
	       "           for a set of images given in an image list. A DPP-file is\n"
	       "           a multi extension fits file with an pixel stamp image and a\n"
	       "           contamination stamp image for each aperture in the\n"
	       "           grism image. Uses the PET file to derive the pixel/contamination\n"
	       "           values for the stamp images and the aperture files\n"
	       "           (AF's) to define a common geometry for the individual objects,\n"
	       "           which are usually located on several images/PET's. The task also\n"
	       "           derives and stores keywords such that aXe_DRIZZLE can create\n"
	       "           coadded images for all objects in the set of DPP-files.\n"
	       "\n"
	       "           The image list is looked for in the current directory.\n"
	       "           The images listed in the image list is looked for in\n"
	       "           $AXE_IMAGE_PATH. The OAF/BAF- and PET-files computed from\n"
	       "           those images is looked for in $AXE_OUTPUT_PATH\n"
	       "           The DPP's are written to $AXE_OUTPUT_PATH\n"
	       "\n"
	       "Usage:\n"
	       "      aXe_DRZPREP [image list filename] [aXe_config1,aXe_config2,..] [options]\n"
	       "\n"
	       "Options:\n"
	       "           -bck  - generating a background DPP from a background PET\n"
	       "\n",RELEASE);
      exit (1);
    }

  fprintf (stdout, "aXe_DRZPREP: Starting...\n\n");

  index = 0;
  strcpy (list_file, argv[++index]);
  strcpy (list_file_path, list_file);

  strcpy (conf_file, argv[++index]);

  // build the aXe inputs structure
  input_list = get_axe_inputs(list_file_path, conf_file);
  //print_axe_inputs(input_list);

  /* Determine if we are using the special bck mode */
  /* In this mode, file names are handled differently */
  if ((opt = get_online_option("bck", argc, argv)))
    bckmode = 1;
  if ((opt = get_online_option("backpet", argc, argv)))
    backpet = 1;
  if ((opt = get_online_option("opt_extr", argc, argv)))
    opt_extr = 1;

  // give feedback onto the screen:
  // report on input and output
  // and also on specific parameter
  fprintf (stdout, "aXe_DRZPREP: Input Image List:           %s\n",
	   list_file_path);
  fprintf (stdout, "aXe_DRZPREP: Configuration file names:   %s\n",
	   conf_file);
  if (opt_extr)
    fprintf (stdout, "aXe_DRZPREP: Storing optimal weighting extentions.\n");
  fprintf(stdout, "\n\n");

  nobjects = 0;
  nobjects2 = 0;

  allobjects = malloc_objectobs();

  fprintf (stdout, "aXe_DRZPREP: Checking input files ...");
  for (in_index=0; in_index < input_list->nitems; in_index++)
    {
      // copy the config file name and the grism file
      // name from the struct to loacal variables
      strcpy(conf_file, input_list->axe_items[in_index].config_file);
      strcpy(grism_file, input_list->axe_items[in_index].grism_file);

      fprintf(stdout, "(%i): %s, %s\n", in_index, conf_file, grism_file);

      // construct the full path to the current configuration fille
      build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);
      for (i=0; i<= DRZMAX; i++)
	{
	  for_grism = check_for_grism (conf_file_path,i);
	  if (!for_grism)
	    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			 "aXe_DRZPREP: Beam %c in configuration file %s has no grism dispersion and can not be drizzled.\n",
			 BEAM(i), conf_file_path);
	}

      conf = get_aperture_descriptor (conf_file_path);

      // compose the full name to the grism image
      build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);

      // check whether the grism file exists
      in_file = fopen(grism_file_path, "r");
      if (!in_file)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_DRZPREP: File %s does not exist!\n", grism_file_path);
      fclose (in_file);

      get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

      // find the full name to the PET file
      if (bckmode)
	{
	  /* Build object PET file name */
	  replace_file_extension (grism_file, PET_file, ".fits",
				  ".BCK.PET.fits", conf->science_numext);
	      build_path (AXE_OUTPUT_PATH, PET_file, PET_file_path);
	}
      else
	{
	  /* Build object PET file name */
	  replace_file_extension (grism_file, PET_file, ".fits",
				  ".PET.fits", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, PET_file, PET_file_path);
	}

      //  Open the OPET file for reading
      fits_open_file (&PET_ptr, PET_file_path, READONLY, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_DRZPREP: Could not open file: %s\n",
		       PET_file_path);
	}

      if (!bckmode)
	// check whether there is quantitative contamination
	quant_cont = check_quantitative_contamination(PET_ptr);
      else
	quant_cont=1;

      // check whether optimal extraction can be done
      if (opt_extr && !quant_cont)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_DRZPREP: PET file: %s was assembled without quantitative contamination.\n",
		     "Optimal extraction is not possible.", PET_file_path);

      // close the PET file
      fits_close_file (PET_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_DRZPREP: Could not" " close PET file: %s\n", PET_file_path);
	}

      if (bckmode && opt_extr)
	{
	  /* Build object PET file name */
	  replace_file_extension (grism_file, PET_file, ".fits",
				  ".PET.fits", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, PET_file, PET_file_path);

	  //  Open the OPET file for reading
	  fits_open_file (&PET_ptr, PET_file_path, READONLY, &f_status);
	  if (f_status)
	    {
	      ffrprt (stderr, f_status);
	      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			   "aXe_DRZPREP: Could not open file: %s\n",
			   PET_file_path);
	    }

	  // check whether there is quantitative contamination
	  quant_cont = check_quantitative_contamination(PET_ptr);

	  if (!quant_cont)
	    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			 "aXe_DRZPREP: PET file: %s was assembled without quantitative contamination.\n",
			 "Optimal extraction is not possible.", PET_file_path);

	  // close the PET file
	  fits_close_file (PET_ptr, &f_status);
	  if (f_status)
	    {
	      ffrprt (stderr, f_status);
	      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			   "aXe_DRZPREP: Could not" " close PET file: %s\n", PET_file_path);
	    }
	}

    }
  fprintf (stdout, "Done.\n\n");

  for (i=0; i < input_list->nitems; i++)
    {
      // copy the config file name and the grism file
      // name from the struct to loacal variables
      strcpy(conf_file, input_list->axe_items[i].config_file);
      strcpy(grism_file, input_list->axe_items[i].grism_file);

      build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);
      conf = get_aperture_descriptor (conf_file_path);

      build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);
      get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

      if (bckmode){
	/* Build aperture file name */
	replace_file_extension (grism_file, aper_file, ".fits",
				".OAF", conf->science_numext);
	build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);
      }
      else{
	if (backpet)
	  replace_file_extension (grism_file, aper_file, ".fits",
				  ".OAF", conf->science_numext);
	else
	  replace_file_extension (grism_file, aper_file, ".fits",
				  ".OAF", conf->science_numext);
	build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);
      }

      fprintf(stdout,
	      "aXe_DRZPREP: Input image file list:    %s\n", list_file_path);
      fprintf(stdout,
	      "aXe_DRZPREP: Input grism image:        %s\n", grism_file_path);
      fprintf(stdout,
	      "aXe_DRZPREP: Input aperture file:      %s\n", aper_file_path);
      fprintf(stdout,
	      "aXe_DRZPREP: Input configuration file: %s\n", conf_file_path);

      obs = load_dummy_observation ();
      /* Loading the object list */
      fprintf (stdout, "aXe_DRZPREP: Loading object aperture list ... ");
      oblist = file_to_object_list_seq (aper_file_path, obs);

      if (oblist != NULL) {
	pixmax = get_npixels (grism_file_path, conf->science_numext);
	fprintf (stdout,"%d objects loaded.\n\n",object_list_size(oblist));
	nobjects2 = add_observation(grism_file_path, conf_file_path,
				    allobjects, nobjects2, oblist,
				    object_list_size(oblist), pixmax,
				    conf->science_numext);
	free_oblist(oblist);
      }
    }

  fprintf (stdout, "aXe_DRZPREP: %i objects in all observations.\n", nobjects2);
  fprintf (stdout, "aXe_DRZPREP: Starting to compute the mean values...");
  f_status =  make_refpoints(conf_file_path, grism_file_path,
			     pixmax, allobjects, nobjects2);
  fprintf (stdout, "Done.\n\n");

  for (in_index=0; in_index < input_list->nitems; in_index++)
    {
      // copy the config file name and the grism file
      // name from the struct to loacal variables
      strcpy(conf_file, input_list->axe_items[in_index].config_file);
      strcpy(grism_file, input_list->axe_items[in_index].grism_file);

      build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);
      conf = get_aperture_descriptor (conf_file_path);

      build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);

      //
      // Start of working on one image
      //
      //
      //

      /* Determine where the various extensions are in the FITS file */
      get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);
      sky_cps = (double)get_float_from_keyword(grism_file_path, conf->science_numext, "SKY_CPS");
      if (isnan(sky_cps))
	sky_cps = 0.0;
      //
      // try to get the descriptor 'sky_cps' from the 'sci'-extension
      //
      exptime = (double)get_float_from_keyword(grism_file_path, conf->science_numext, conf->exptimekey);
      if (isnan(exptime))
	exptime = (double)get_float_from_keyword(grism_file_path, 1, conf->exptimekey);
      if (isnan(exptime))
	exptime = 1.0;

      if (bckmode)
	{
	  /* Build aperture file name */
	  //	  replace_file_extension (grism_file, aper_file, ".fits",
	  //				  ".BAF", conf->science_numext);
	  replace_file_extension (grism_file, aper_file, ".fits",
				  ".OAF", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);

	  /* Build object PET file name */
	  replace_file_extension (grism_file, PET_file, ".fits",
				  ".BCK.PET.fits", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, PET_file, PET_file_path);

	  if (opt_extr)
	    {
	      /* Build second PET file name */
	      replace_file_extension (grism_file, SEC_file, ".fits",
				      ".PET.fits", conf->science_numext);
	      build_path (AXE_OUTPUT_PATH, SEC_file, SEC_file_path);
	    }
	}
      else
	{
	  /* Build aperture file name */
	  //	  replace_file_extension (grism_file, aper_file, ".fits",
	  //				  ".BAF", conf->science_numext);
	  replace_file_extension (grism_file, aper_file, ".fits",
				  ".OAF", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);

	  /* Build object PET file name */
	  replace_file_extension (grism_file, PET_file, ".fits",
				  ".PET.fits", conf->science_numext);
	  build_path (AXE_OUTPUT_PATH, PET_file, PET_file_path);
	}

      if ( (opt = get_online_option ("outputroot", argc, argv)) ){
	strcpy (outputroot, opt);
      }
      else{
	replace_file_extension (PET_file, outputroot, ".PET.fits", "",-1);
      }

      build_path (AXE_OUTPUT_PATH, outputroot, outputroot_path);
      sprintf(output_path,"%s.DPP.fits",outputroot_path);

      fprintf (stdout, "aXe_DRZPREP: Input PET file name:       %s\n",
	       PET_file_path);
      fprintf (stdout, "aXe_DRZPREP: Input Aperture file name:  %s\n",
	       aper_file_path);
      fprintf (stdout, "aXe_DRZPREP: Name of DPP file :         %s\n",
	       output_path);

      obs = load_dummy_observation ();
      /* Loading the object list */
      fprintf (stdout, "aXe_DRZPREP: Loading object aperture list...");
      oblist = file_to_object_list_seq (aper_file_path, obs);
      fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));

      //  Open the OPET file for reading
      fits_open_file (&PET_ptr, PET_file_path, READONLY, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_DRZPREP: Could not open file: %s\n",
		       PET_file_path);
	}

      if (!bckmode)
	// check whether there is quantitative contamination
	quant_cont = check_quantitative_contamination(PET_ptr);
      else
	quant_cont=1;

      if (bckmode && opt_extr)
	{
	  //  Open the OPET file for reading
	  fits_open_file (&SEC_ptr, SEC_file_path, READONLY, &f_status);
	  if (f_status)
	    {
	      ffrprt (stderr, f_status);
	      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			   "aXe_DRZPREP: Could not open file: %s\n",
			   SEC_file_path);
	    }
	}

      DPP_ptr = create_FITSimage_opened (output_path, 1);
      /* Copy the header info from the grism image */
      {
	FITScards *cards;
	if (bckmode && opt_extr)
	  cards = get_FITS_cards_opened (SEC_ptr);
	else
	  cards = get_FITS_cards_opened (PET_ptr);
	put_FITS_cards_opened(DPP_ptr,cards);
	free_FITScards(cards);
      }

      i = 0;
      if (oblist!=NULL)
	{
	  while (1)
	    {

	      /* Get the PET for this object */
	      pri_PET = get_ALL_from_next_in_PET(PET_ptr, &aperID, &beamID);
	      if ((aperID==-1) && (beamID==-1)) break;

	      if (bckmode && opt_extr)
		sec_PET = get_ALL_from_next_in_PET(SEC_ptr, &aperID, &beamID);
	      else
		sec_PET = NULL;

	      if (beamID >= 0 && beamID <= DRZMAX){
		//	      if (beamID == 0){
		fprintf (stdout, "aXe_DRZPREP: BEAM_%d%c.", aperID, BEAM(beamID));
		objindex =  find_object_in_object_list(oblist,aperID);

		sprintf (label, "%s.%d%c.ps/CPS", outputroot_path,
			 oblist[objindex]->ID, BEAM (oblist[objindex]->beams[beamID].ID));

		pixel.x = oblist[objindex]->beams[beamID].refpoint.x - conf->refx;
		pixel.y = oblist[objindex]->beams[beamID].refpoint.y - conf->refy;
		disp = get_dispstruct_at_pos(conf_file_path, 1,
					     oblist[objindex]->beams[beamID].ID,pixel);
		trace = oblist[objindex]->beams[beamID].spec_trace;

		outref = get_mean_refpoint(allobjects, nobjects2, oblist[objindex]->ID,
					   &boxwidth, &boxheight, &relx, &rely, &objwidth,
					   &orient, &cdref, &cdscale, &cdmeanscale,
					   &sprefreso, &spreso, &spmeanreso, &trlength);

		minxy_PET = get_minxy_from_PET(pri_PET);
		relx = oblist[objindex]->beams[beamID].refpoint.x - minxy_PET.x;
		rely = oblist[objindex]->beams[beamID].refpoint.y - minxy_PET.y;

		//*************************************************
		// patch to correct the reference point in case
		// that the the trace descritpion
		// does have a non negligeable first order term!
		gaga = oblist[objindex]->beams[beamID].spec_trace->data;
		rely = rely + gaga[1];
		//**************************************************

		outdisp = get_dispstruct_at_pos(conf_file_path, 1,
						  oblist[objindex]->beams[beamID].ID,outref);

		cdcorr = cdscale / cdref;
		spcorr = spreso  / sprefreso;

		// determine the trace length
		// NOTE: putting in this line is wrong, but opens to make
		//       drizzle results as in aXe-1.7
		//trlength = get_beam_trace_length(oblist[objindex]->beams[beamID]);

		drizzcoeffs = get_drizzle_coeffs( disp, trace, boxwidth, boxheight,
						  trlength, relx, rely, conf, orient,
						  outdisp, cdcorr, sprefreso, spmeanreso);

		drizzle_width = cdcorr * get_drizzle_width(oblist[objindex],beamID,trace);
		objwidth      = cdmeanscale / cdref * objwidth;

		refwave_pos = get_refwave_position( disp,  trace, pixel, conf);
		refwave_pos.x = refwave_pos.x - minxy_PET.x;
		refwave_pos.y = refwave_pos.y - minxy_PET.y;

		trlength = gsl_matrix_get(drizzcoeffs, 0,10);
		dimension =  get_drzprep_dim(pri_PET, oblist[objindex]->beams[beamID].width,
					     boxwidth, boxheight);
		{
		  px_point tmp_in;
		  px_point tmp_out;

		  d_point new_pos;

		  tmp_in.x = dimension.xsize;
		  tmp_in.y = dimension.ysize;

		  tmp_out.x = (int)trlength;
		  tmp_out.y = 2*(int)ceil(objwidth) + 10;
		  //		    gsl_matrix_fprintf (stdout, drizzcoeffs, "%f");

		  gsl_matrix_set(drizzcoeffs, 0,10,0.0);
		  new_pos = get_drz_position_free(refwave_pos, drizzcoeffs, tmp_in, tmp_out);
		}

		drzprep_stamps = stamp_img_drzprep(opt_extr, pri_PET, sec_PET,
						   oblist[objindex]->beams[beamID].width,
						   -1000000.0, quant_cont, dimension,
						   drizzcoeffs, exptime, sky_cps,
						   (double)conf->rdnoise, bckmode);

		// store the count-extension
		sprintf (label, "BEAM_%d%c",oblist[objindex]->ID,
			 BEAM (oblist[objindex]->beams[beamID].ID));
		gsl_to_FITSimage_opened (drzprep_stamps->counts, DPP_ptr ,0,label);
		cards = drzinfo_to_FITScards(oblist[objindex],beamID,
					     outref, conf, drizzcoeffs,
					     trlength, relx, rely,objwidth,
					     refwave_pos, sky_cps,
					     drizzle_width, cdref, spcorr);
		put_FITS_cards_opened(DPP_ptr,cards);

		// store the error-extension
		sprintf (label, "ERR_%d%c",oblist[objindex]->ID,
			 BEAM (oblist[objindex]->beams[beamID].ID));
		gsl_to_FITSimage_opened (drzprep_stamps->error, DPP_ptr ,0,label);
		put_FITS_cards_opened(DPP_ptr,cards);

		// store the contamination-extension
		sprintf (label, "CONT_%d%c",oblist[objindex]->ID,
			 BEAM (oblist[objindex]->beams[beamID].ID));
		gsl_to_FITSimage_opened (drzprep_stamps->cont, DPP_ptr ,0,label);
		put_FITS_cards_opened(DPP_ptr,cards);

		//		  if (opt_extr)
		if (drzprep_stamps->model)
		  {
		    // store the model-extension
		    sprintf (label, "MOD_%d%c",oblist[objindex]->ID,
			     BEAM (oblist[objindex]->beams[beamID].ID));
		    gsl_to_FITSimage_opened (drzprep_stamps->model, DPP_ptr ,0,label);
		    put_FITS_cards_opened(DPP_ptr,cards);
		  }

		if (drzprep_stamps->vari)
		  {
		    // store the variance extension
		    sprintf (label, "VAR_%d%c",oblist[objindex]->ID,
			     BEAM (oblist[objindex]->beams[beamID].ID));
		    gsl_to_FITSimage_opened (drzprep_stamps->vari, DPP_ptr ,0,label);
		    put_FITS_cards_opened(DPP_ptr,cards);
		  }

		// release memory
		free_drzprep(drzprep_stamps);
		free_FITScards(cards);
		free_dispstruct(disp);
		free_dispstruct(outdisp);
		gsl_matrix_free(drizzcoeffs);

		fprintf (stdout, " Done.\n");
		i++;
	      }

	      if (pri_PET!=NULL) free(pri_PET);
	      if (sec_PET!=NULL) free(sec_PET);
	    }
	}
      if (oblist!=NULL) free_oblist (oblist);

      fits_close_file (DPP_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_DRZPREP: Could not" " close DPP file: %s\n", output_path);
	}

      fits_close_file (PET_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_DRZPREP: Could not" " close PET file: %s\n", PET_file_path);
	}

      if (bckmode && opt_extr)
	{
	  fits_close_file (SEC_ptr, &f_status);
	  if (f_status)
	    {
	      ffrprt (stderr, f_status);
	      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			   "aXe_DRZPREP: Could not" " close PET file: %s\n", SEC_file_path);
	    }
	}

      fprintf (stdout, "aXe_DRZPREP: %s Done...\n\n", output_path);
      //
      // End of working on one image
      //
      //
      //

    }

  free_axe_inputs(input_list);
  free_objectobs(allobjects);
  exit (0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#include <string>

#include <gsl/gsl_complex.h>

#include "Expression.h"
#include "UnitConversionExpression/Units.h"

class BaseUnitExpression: public Expression {

    UnitType type;
    std::string abbr;
    double val;

    BaseUnitExpression(const std::string& unit, double val);
    BaseUnitExpression(UnitType type, const std::string& unit, double val);

    public:
        static expression construct(const std::string& unit, double val = 1);
        static expression construct(UnitType type, const std::string& unit, double val);

        std::string repr() const override;
        int id() const override;

        EXPRESSION_OVERRIDES

        friend expression operator*(const BaseUnitExpression& unit1, const BaseUnitExpression& unit2);
        friend expression operator*(const BaseUnitExpression& unit1, const expression expr);
        friend expression operator/(const BaseUnitExpression& unit1, const BaseUnitExpression& unit2);
        friend expression operator/(const BaseUnitExpression& unit1, const expression expr);
        friend expression operator^(const BaseUnitExpression& unit1, const expression expr);

        friend expression convert(const BaseUnitExpression& from, const BaseUnitExpression& to);
};

class ConvertedUnitExpression: public Expression {

    UnitType type;
    std::string abbr;
    double val;

    ConvertedUnitExpression(UnitType type, const std::string& unit, double val);

    public:
        static expression construct(UnitType type, const std::string& unit, double val);

        std::string repr() const override;
        int id() const override;

        EXPRESSION_OVERRIDES
};

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
#ifndef PARALLEL_PROCESSING_H_2FVRLMCH
#define PARALLEL_PROCESSING_H_2FVRLMCH

#include <chrono>
#include <gsl/gsl>
#include <iomanip>
#include <ios>
#include <iostream>
#include <sens_loc/util/console.h>
#include <sens_loc/util/progress_bar_observer.h>
#include <taskflow/taskflow.hpp>
#include <type_traits>

namespace sens_loc::apps {

/// Helper function that processes a range of files based on index.
/// The boolean function \c f is applied to each function. Error handling
/// and reporting is done if \c f returns \c false.
///
/// \tparam BoolFunction Apply this functor for each index.
/// \param start,end inclusive range of integers for the files
/// \param f functor that is applied for each index
template <typename BoolFunction>
bool parallel_indexed_file_processing(int          start,
                                      int          end,
                                      BoolFunction f) noexcept {
    static_assert(std::is_nothrow_invocable_r_v<bool, BoolFunction, int>,
                  "Functor needs to be noexcept callable and return bool!");

    try {
        if (start > end)
            std::swap(start, end);

        int total_tasks = end - start + 1;

        tf::Executor executor;
        executor.make_observer<util::progress_bar_observer>(total_tasks);
        tf::Taskflow tf;

        bool batch_success = true;
        int  fails         = 0;

        tf.parallel_for(
            start, end + 1, 1, [&batch_success, &fails, &f](int idx) {
                const bool success = f(idx);
                if (!success) {
                    auto s = synced();
                    fails++;
                    std::cerr << util::err{};
                    std::cerr << "Could not process index \""
                              << rang::style::bold << idx << "\""
                              << rang::style::reset << "!" << std::endl;
                    batch_success = false;
                }
            });

        const auto before = std::chrono::steady_clock::now();
        executor.run(tf).wait();
        const auto after = std::chrono::steady_clock::now();
        const auto dur_deci_seconds =
            std::chrono::duration_cast<std::chrono::duration<long, std::centi>>(
                after - before);
        std::cout << std::endl;

        Ensures(fails >= 0);

        using namespace std::chrono;
        {
            auto s = synced();
            std::cerr << util::info{};
            std::cerr << "Processing " << rang::style::bold
                      << std::abs(end - start) + 1 - fails << rang::style::reset
                      << " images took " << rang::style::bold << std::fixed
                      << std::setprecision(2)
                      << (dur_deci_seconds.count() / 100.) << rang::style::reset
                      << " seconds!\n";
        }

        if (fails > 0) {
            auto s = synced();
            std::cerr << util::warn{} << "Encountered " << rang::style::bold
                      << fails << rang::style::reset << " problematic files!\n";
        }

        return batch_success;
    } catch (...) {
        auto s = synced();
        std::cerr << util::err{} << "System error in batch processing!\n";
        return false;
    }
}

}  // namespace sens_loc::apps

#endif /* end of include guard: PARALLEL_PROCESSING_H_2FVRLMCH */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


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

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#pragma GCC diagnostic ignored "-Weffc++"

#include "fabric_transport.h"
#include "mcas_client_config.h"
#include "protocol.h"
#include "protocol_ostream.h"

#include <api/fabric_itf.h>
#include <api/mcas_itf.h>
#include <common/exceptions.h>
#include <common/utils.h>
#include <common/byte_buffer.h>
#include <common/string_view.h>
#include <gsl/pointers>
#include <sys/mman.h>
#include <sys/uio.h>
#include <unistd.h>
#include <gnutls/gnutls.h>
#include <gnutls/x509.h>
#include <gnutls/crypto.h>

#include <boost/numeric/conversion/cast.hpp>
#include <map>
#include <set>
#include <tuple>

/* Enable this to introduce locks to prevent re-entry by multiple
   threads.  The client is not re-entrant because of the state machine
   and multi-packet operations
*/
#define THREAD_SAFE_CLIENT

#define CAFILE "/etc/ssl/certs/ca-bundle.trust.crt"

#ifdef THREAD_SAFE_CLIENT
#define API_LOCK() std::lock_guard<std::mutex> g(_api_lock);
#else
#define API_LOCK()
#endif

namespace mcas
{

namespace client
{

struct TLS_transport
{
  static unsigned debug_level();
  static ssize_t gnutls_pull_func(gnutls_transport_ptr_t, void*, size_t);
  static ssize_t gnutls_vec_push_func(gnutls_transport_ptr_t, const giovec_t * , int );
  static int gnutls_pull_timeout_func(gnutls_transport_ptr_t, unsigned int);
};

struct async_buffer_set_t;
struct iob_free;

static const char *env_scbe = ::getenv("SHORT_CIRCUIT_BACKEND");

/* Adaptor point for other transports */
using Connection_base = mcas::client::Fabric_transport;

/**
 * Client side connection handler
 *
 */
class Connection_handler : public Connection_base {
  using iob_ptr        = std::unique_ptr<client::Fabric_transport::buffer_t, iob_free>;
  using locate_element = protocol::Message_IO_response::locate_element;

public:
  friend struct TLS_transport;

  using memory_region_t = typename Transport::memory_region_t;
  template <typename T>
    using basic_string_view = common::basic_string_view<T>;
  using byte = common::byte;

  /**
   * Constructor
   *
   * @param debug_level
   * @param connection
   * @param patience Time to wait (in seconds) for single fabric post to complete
   * @param other Additional configuration string
   *
   * @return
   */
  Connection_handler(const unsigned debug_level,
                     Connection_base::Transport *connection,
                     const unsigned patience,
                     const std::string other);

  ~Connection_handler();

private:
  enum State {
              INITIALIZE,
              HANDSHAKE_SEND,
              HANDSHAKE_GET_RESPONSE,
              SHUTDOWN,
              STOPPED,
              READY,
  };

  State _state = State::INITIALIZE;

  template <typename MT>
    status_t invoke_ado_common(
      const iob_ptr & iobs
      , const MT *msg
      , std::vector<component::IMCAS::ADO_response>& out_response
      , unsigned int flags
    );

public:
  using pool_t = uint64_t;

  void bootstrap()
  {
    set_state(INITIALIZE);
    while (tick() > 0)
      ;
  }

  void shutdown()
  {
    set_state(SHUTDOWN);
    while (tick() > 0) sleep(1);
  }

  pool_t open_pool(const std::string name,
                   const unsigned int flags,
                   const addr_t base);

  pool_t create_pool(const std::string  name,
                     const size_t       size,
                     const unsigned int flags,
                     const uint64_t     expected_obj_count,
                     const addr_t       base);

  status_t close_pool(const pool_t pool);

  status_t delete_pool(const std::string &name);

  status_t delete_pool(const pool_t pool);

  status_t configure_pool(const component::IKVStore::pool_t pool, const std::string &json);

  status_t put(const pool_t       pool,
               const std::string  key,
               const void *       value,
               const size_t       value_len,
               const unsigned int flags);

  status_t put(const pool_t       pool,
               const void *       key,
               const size_t       key_len,
               const void *       value,
               const size_t       value_len,
               const unsigned int flags);

  status_t put_direct(pool_t                               pool,
                      const void *                         key,
                      size_t                               key_len,
                      const void *                         value,
                      size_t                               value_len,
                      component::Registrar_memory_direct * rmd,
                      component::IKVStore::memory_handle_t handle,
                      unsigned int                         flags);

  status_t async_put(const pool_t                      pool,
                     const void *                      key,
                     const size_t                      key_len,
                     const void *                      value,
                     size_t                            value_len,
                     component::IMCAS::async_handle_t &out_handle,
                     unsigned int                      flags);

  status_t async_put_direct(const component::IMCAS::pool_t       pool,
                            const void *                         key,
                            size_t                               key_len,
                            const void *                         value,
                            size_t                               value_len,
                            component::IMCAS::async_handle_t &   out_handle,
                            component::Registrar_memory_direct * rmd,
                            component::IKVStore::memory_handle_t handle = component::IMCAS::MEMORY_HANDLE_NONE,
                            unsigned int                         flags  = component::IMCAS::FLAGS_NONE);

  status_t async_get_direct(const component::IMCAS::pool_t       pool,
                            const void *                         key,
                            size_t                               key_len,
                            void *                               value,
                            size_t &                             value_len,
                            component::IMCAS::async_handle_t &   out_handle,
                            component::Registrar_memory_direct * rmd,
                            component::IKVStore::memory_handle_t handle = component::IMCAS::MEMORY_HANDLE_NONE,
                            unsigned int                         flags  = component::IMCAS::FLAGS_NONE);

  status_t check_async_completion(component::IMCAS::async_handle_t &handle);

  status_t get(const pool_t pool, const std::string &key, std::string &value);

  status_t get(const pool_t pool, const std::string &key, void *&value, size_t &value_len);

  status_t get_direct(const pool_t                         pool,
                      const void *                         key,
                      size_t                               key_len,
                      void *                               value,
                      size_t &                             out_value_len,
                      component::Registrar_memory_direct * rmd,
                      component::IKVStore::memory_handle_t handle = component::IKVStore::HANDLE_NONE);

  status_t get_direct_offset(pool_t                              pool,
                             std::size_t                         offset,
                             std::size_t &                       length,
                             void *                              buffer,
                             component::Registrar_memory_direct *rmd,
                             component::IMCAS::memory_handle_t   handle);

  status_t async_get_direct_offset(pool_t                              pool,
                                   std::size_t                         offset,
                                   std::size_t &                       length,
                                   void *                              buffer,
                                   component::IMCAS::async_handle_t &  out_handle,
                                   component::Registrar_memory_direct *rmd,
                                   component::IMCAS::memory_handle_t   handle);

  status_t put_direct_offset(pool_t                              pool,
                             std::size_t                         offset,
                             std::size_t &                       length,
                             const void *                        buffer,
                             component::Registrar_memory_direct *rmd,
                             component::IMCAS::memory_handle_t   handle);

  status_t async_put_direct_offset(pool_t                              pool,
                                   std::size_t                         offset,
                                   std::size_t &                       length,
                                   const void *                        buffer,
                                   component::IMCAS::async_handle_t &  out_handle,
                                   component::Registrar_memory_direct *rmd,
                                   component::IMCAS::memory_handle_t   handle);

  status_t erase(const pool_t pool, const std::string &key);

  status_t async_erase(const component::IMCAS::pool_t    pool,
                       const std::string &               key,
                       component::IMCAS::async_handle_t &out_handle);

  uint64_t key_hash(const void *key, const size_t key_len);

  uint64_t auth_id() const
  {
    /* temporary */
    auto     env = getenv("MCAS_AUTH_ID");
    uint64_t auth_id;
    if (env) {
      auto id = std::strtoll(env, nullptr, 10);
      auth_id = boost::numeric_cast<uint64_t>(id);
    }
    else {
      auth_id = boost::numeric_cast<uint64_t>(getuid());
    }

    return auth_id;
  }

  size_t count(const pool_t pool);

  status_t get_attribute(const component::IKVStore::pool_t    pool,
                         const component::IKVStore::Attribute attr,
                         std::vector<uint64_t> &              out_attr,
                         const std::string *                  key);

  status_t get_statistics(component::IMCAS::Shard_stats &out_stats);

  status_t find(const component::IKVStore::pool_t pool,
                const std::string &               key_expression,
                const offset_t                    offset,
                offset_t &                        out_matched_offset,
                std::string &                     out_matched_key);

  status_t invoke_ado(const component::IMCAS::pool_t               pool,
                      basic_string_view<byte>                      key,
                      basic_string_view<byte>                      request,
                      const unsigned int                           flags,
                      std::vector<component::IMCAS::ADO_response> &out_response,
                      const size_t                                 value_size);

  status_t invoke_ado_async(const component::IMCAS::pool_t               pool,
                            basic_string_view<byte>                      key,
                            basic_string_view<byte>                      request,
                            const component::IMCAS::ado_flags_t          flags,
                            std::vector<component::IMCAS::ADO_response> &out_response,
                            component::IMCAS::async_handle_t &           out_async_handle,
                            const size_t                                 value_size);

  status_t invoke_put_ado(const component::IMCAS::pool_t               pool,
                          basic_string_view<byte>                      key,
                          basic_string_view<byte>                      request,
                          basic_string_view<byte>                      value,
                          size_t                                       root_len,
                          const unsigned int                           flags,
                          std::vector<component::IMCAS::ADO_response> &out_response);

  status_t invoke_put_ado_async(const component::IMCAS::pool_t                  pool,
                                const basic_string_view<byte>                   key,
                                const basic_string_view<byte>                   request,
                                const basic_string_view<byte>                   value,
                                const size_t                                    root_len,
                                const component::IMCAS::ado_flags_t             flags,
                                std::vector<component::IMCAS::ADO_response>&    out_response,
                                component::IMCAS::async_handle_t&               out_async_handle);

  bool check_message_size(size_t size) const { return size > _max_message_size; }

  status_t receive_and_process_ado_response(
    const iob_ptr & iobr
    , std::vector<component::IMCAS::ADO_response> & out_response
  );

private:
  /**
   * FSM tick call
   *
   */
  int tick();

  /**
   * FSM state change
   *
   * @param s State to change to
   */
  inline void set_state(State s)
  {
    if (2 < debug_level()) {
      static const std::map<State, const char *> m{
                                                   {INITIALIZE, "INITIALIZE"},
                                                   {HANDSHAKE_SEND, "HANDSHAKE_SEND"},
                                                   {HANDSHAKE_GET_RESPONSE, "HANDSHAKE_GET_RESPONSE"},
                                                   {SHUTDOWN, "SHUTDOWN"},
                                                   {STOPPED, "STOPPED"},
                                                   {READY, "READY"},
      };
      PLOG("Client state %s -> %s", m.find(_state)->second, m.find(s)->second);
    }
    _state = s;
  } /* we could add transition checking later */

  void start_tls();

  template <typename MT>
  void msg_send_log(const MT *m, const void *context, const char *desc) { msg_send_log(1, m, context, desc); }

  template <typename MT>
  void msg_send_log(const unsigned level, const MT *msg, const void *context, const char *desc)
  {
    msg_log(level, msg, context, desc, "SEND");
  }

  template <typename MT>
  void msg_log(const unsigned level_, const MT *msg_, const void *context_, const char *desc_, const char *direction_)
  {
    if ( level_ < debug_level() )
      {
        std::ostringstream m;
        m << *msg_;
        PLOG("%s (%p) (%s) %s", direction_, context_, desc_, m.str().c_str());
      }
  }

  template <typename MT>
  void msg_recv_log(const MT *m, const void *context, const char *desc) { msg_recv_log(1, m, context, desc); }

  template <typename MT>
  void msg_recv_log(const unsigned level, const MT *msg, const void *context, const char *desc)
  {
    msg_log(level, msg, context, desc, "RECV");
  }

public:

  template <typename MT>
  gsl::not_null<const MT *>msg_recv(const buffer_t *iob, const char *desc)
  {
    /*
     * First, cast the response buffer to Message (checking version).
     * Second, cast the Message to a specific message type (checking message type).
     */
    const auto *const msg = mcas::protocol::message_cast(iob->base());
    const auto *const response_msg = msg->ptr_cast<MT>();
    msg_recv_log(response_msg, iob, desc);
    return response_msg;
  }

  template <typename MT>
  void sync_inject_send(buffer_t *iob, const MT *msg, std::size_t size, const char *desc)
  {
    msg_send_log(msg, iob, desc);
    Connection_base::sync_inject_send(iob, size);
  }

  template <typename MT>
  void sync_inject_send(buffer_t *iob, const MT *msg, const char *desc)
  {
    sync_inject_send(iob, msg, msg->msg_len(), desc);
  }

  template <typename MT>
  void sync_send(buffer_t *iob, const MT *msg, const char *desc)
  {
    msg_send_log(msg, iob, desc);
    Connection_base::sync_send(iob);
  }

  void sync_inject_send(buffer_t *iob, std::size_t size)
  {
    Connection_base::sync_inject_send(iob, size);
  }

private:
  /* unused */
#if 0
  void post_send(buffer_t *iob, const protocol::Message_IO_request *msg, buffer_external *iob_extra, const char *desc)
  {
    msg_send_log(msg, iob, desc);
    Connection_base::post_send(iob, iob_extra);
  }
#endif
  template <typename MT>
  void post_send(buffer_t *iob, const MT *msg, const char *desc)
  {
    msg_send_log(2, msg, iob, desc);
    Connection_base::post_send(iob);
  }

  template <typename MT>
  void post_send(const ::iovec *first, const ::iovec *last, void **descriptors, void *context, const MT *msg, const char *desc)
  {
    msg_send_log(msg, context, desc);
    Connection_base::post_send(first, last, descriptors, context);
  }

  /**
   * Put used when the value exceeds the size of the basic
   * IO buffer (e.g., 2MB).  This version of put will perform
   * a two-stage exchange, advance notice and then value
   *
   * @param pool Pool identifier
   * @param key Key
   * @param key_len Key length
   * @param value_len Value length
   * @param flags ignored?
   *
   * @return
   */
  std::tuple<uint64_t, uint64_t> put_locate(const pool_t   pool,
                                            const void *   key,
                                            const size_t   key_len,
                                            const size_t   value_len,
                                            const unsigned flags);

  component::IMCAS::async_handle_t put_locate_async(pool_t                              pool,
                                                    const void *                        key,
                                                    size_t                              key_len,
                                                    const void *                        value,
                                                    size_t                              value_len,
                                                    component::Registrar_memory_direct *rmd,
                                                    void *                              desc,
                                                    unsigned                            flags);

  std::tuple<uint64_t, uint64_t, std::size_t> get_locate(const pool_t   pool,
                                                         const void *   key,
                                                         const size_t   key_len,
                                                         const unsigned flags);

  component::IMCAS::async_handle_t get_locate_async(pool_t                              pool,
                                                    const void *                        key,
                                                    size_t                              key_len,
                                                    void *                              value,
                                                    size_t &                            value_len,
                                                    component::Registrar_memory_direct *rmd,
                                                    void *                              desc,
                                                    unsigned                            flags);
public:
  iob_ptr make_iob_ptr(buffer_t::completion_t);
  iob_ptr make_iob_ptr_recv();
  iob_ptr make_iob_ptr_send();
  iob_ptr make_iob_ptr_write();
  iob_ptr make_iob_ptr_read();

private:
  static void send_complete(void *, buffer_t *iob);
  static void recv_complete(void *, buffer_t *iob);
  static void write_complete(void *, buffer_t *iob);
  static void read_complete(void *, buffer_t *iob);

  std::tuple<uint64_t, std::vector<locate_element>> locate(pool_t pool, std::size_t offset, std::size_t size);

  component::IMCAS::async_handle_t get_direct_offset_async(pool_t                              pool,
                                                           std::size_t                         offset,
                                                           void *                              buffer,
                                                           std::size_t &                       length,
                                                           component::Registrar_memory_direct *rmd,
                                                           void *                              desc);

  component::IMCAS::async_handle_t put_direct_offset_async(pool_t                              pool,
                                                           std::size_t                         offset,
                                                           const void *                        buffer,
                                                           std::size_t &                       length,
                                                           component::Registrar_memory_direct *rmd,
                                                           void *                              desc);

private:
#ifdef THREAD_SAFE_CLIENT
  std::mutex _api_lock;
#endif

  bool     _exit;
  uint64_t _request_id;

public: /* for async "move_along" processing */
  uint64_t request_id() { return ++_request_id; }

private:
  size_t _max_message_size;
  size_t _max_inject_size;

  struct options_s {
    bool short_circuit_backend;
    unsigned tls   : 1;
    unsigned hmac : 1;

    options_s()
      : short_circuit_backend(env_scbe && env_scbe[0] == '1'), tls(0), hmac(0)
    {}
  };

  options_s _options;

  gnutls_certificate_credentials_t _xcred;
  gnutls_priority_t                _priority;
  gnutls_session_t                 _session;
  common::Byte_buffer              _tls_buffer;

};

}  // namespace client
}  // namespace mcas

#pragma GCC diagnostic pop
#endif


text:
//===--- Sudoku/precompiled_basic.h                                     ---===//
//
// Lists the standard system includes consumed by this library,
// or project specific include files that are used frequently, but
// are changed infrequently.
//
// - To be used as part of a precompiled header configuration for a project
//   consuming this library.
// - Where frequent changes to the Sudoku library code are expected.
// - Not required when the library itself is added to a precompiled header.
//
//===----------------------------------------------------------------------===//
#pragma once

// STL containers
#include <array>
#include <bitset>
#include <vector>
// STL
#include <algorithm>
#include <iterator>
#include <limits>
#include <memory>
#include <stdexcept>
#include <utility>

// Libraries
#include <gsl/gsl>

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


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

#include <config.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  GSL_ERROR (
"the IEEE interface for this platform is unsupported or could not be "
"determined at configure time\n", GSL_EUNSUP) ;
}


text:
#ifndef EXECUTION_HANDLER_INCLUDE
#define EXECUTION_HANDLER_INCLUDE

#include <map>

#include <gsl/gsl>

#include "executionContent.h"

namespace execHelper {
namespace test {
namespace baseUtils {
class ExecutionHandler {
  public:
    void add(const std::string& key, ExecutionContent&& content) noexcept;

    const ExecutionContent& at(const std::string& key) const noexcept;

  private:
    using ExecutionContentCollection = std::map<std::string, ExecutionContent>;

    /**
     * \brief   Used for handling an execution iteration
     */
    class ExecutionHandlerIterationRAII {
      public:
        explicit ExecutionHandlerIterationRAII(
            gsl::not_null<ExecutionContentCollection*> outputs);
        ~ExecutionHandlerIterationRAII();

        ExecutionHandlerIterationRAII(
            const ExecutionHandlerIterationRAII& other) = default;
        ExecutionHandlerIterationRAII(ExecutionHandlerIterationRAII&& other) =
            default;

        ExecutionHandlerIterationRAII&
        operator=(const ExecutionHandlerIterationRAII& other) = default;
        ExecutionHandlerIterationRAII&
        operator=(ExecutionHandlerIterationRAII&& other) noexcept =
            default; // NOLINT(misc-noexcept-move-constructor)

        const ExecutionContent& at(const std::string& key) const noexcept;

      private:
        gsl::not_null<ExecutionContentCollection*> m_outputs;
    };

    ExecutionContentCollection m_outputs;

  public:
    ExecutionHandlerIterationRAII startIteration() noexcept;
};
} // namespace baseUtils
} // namespace test
} // namespace execHelper

#endif /* EXECUTION_HANDLER_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


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

/* the code in this module performs Romberg integration */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#define ROMBERG_PRINT_ROW(i, r) \
  do { \
    size_t jj; \
    fprintf(stderr, "R[%zu] = ", i); \
    for (jj = 0; jj <= i; ++jj) \
      fprintf(stderr, "%.8e ", r[jj]); \
    fprintf(stderr, "\n"); \
  } while (0)

gsl_integration_romberg_workspace *
gsl_integration_romberg_alloc(const size_t n)
{
  gsl_integration_romberg_workspace *w;

  /* check inputs */
  if (n < 1)
    {
      GSL_ERROR_VAL ("workspace size n must be at least 1", GSL_EDOM, 0);
    }

  w = calloc(1, sizeof(gsl_integration_romberg_workspace));
  if (w == NULL)
    {
      GSL_ERROR_VAL ("unable to allocate workspace", GSL_ENOMEM, 0);
    }

  /* ceiling on n, since the number of points is 2^n + 1 */
  w->n = GSL_MIN(n, 30);

  w->work1 = malloc(w->n * sizeof(double));
  if (w->work1 == NULL)
    {
      gsl_integration_romberg_free(w);
      GSL_ERROR_VAL ("unable to allocate previous row", GSL_ENOMEM, 0);
    }

  w->work2 = malloc(w->n * sizeof(double));
  if (w->work2 == NULL)
    {
      gsl_integration_romberg_free(w);
      GSL_ERROR_VAL ("unable to allocate current row", GSL_ENOMEM, 0);
    }

  return w;
}

void
gsl_integration_romberg_free(gsl_integration_romberg_workspace * w)
{
  if (w->work1)
    free(w->work1);

  if (w->work2)
    free(w->work2);

  free(w);
}

int
gsl_integration_romberg(const gsl_function * f, const double a, const double b,
                        const double epsabs, const double epsrel, double * result,
                        size_t * neval, gsl_integration_romberg_workspace * w)
{
  if (epsabs < 0.0)
    {
      GSL_ERROR("epsabs must be non-negative", GSL_EDOM);
    }
  else if (epsrel < 0.0)
    {
      GSL_ERROR("epsrel must be non-negative", GSL_EDOM);
    }
  else
    {
      const size_t n = w->n;
      double *Rp = &(w->work1[0]); /* previous row */
      double *Rc = &(w->work2[0]); /* current row */
      double *Rtmp;
      double h = 0.5 * (b - a);    /* step size */
      size_t i;

      /* R(0,0) */
      Rp[0] = h * (GSL_FN_EVAL(f, a) + GSL_FN_EVAL(f, b));
      *neval = 2;

      /*ROMBERG_PRINT_ROW((size_t) 0, Rp);*/

      for (i = 1; i < n; ++i)
        {
          size_t j;
          double sum = 0.0;
          double err;
          double four_j;         /* 4^j */
          size_t two_i = 1 << i; /* 2^i */

          for (j = 1; j < two_i; j += 2)
            {
              sum += GSL_FN_EVAL(f, a + j * h);
              ++(*neval);
            }

          /* R(i,0) */
          Rc[0] = sum * h + 0.5 * Rp[0];

          four_j = 4.0;
          for (j = 1; j <= i; ++j)
            {
              Rc[j] = (four_j * Rc[j - 1] - Rp[j - 1]) / (four_j - 1.0);
              four_j *= 4.0;
            }

          /*ROMBERG_PRINT_ROW(i, Rc);*/

          /*
           * compute difference between current and previous result and
           * check for convergence
           */
          err = fabs(Rc[i] - Rp[i - 1]);
          if ((err < epsabs) || (err < epsrel * fabs(Rc[i])))
            {
              *result = Rc[i];
              return GSL_SUCCESS;
            }

          /* swap Rp and Rc */
          Rtmp = Rp;
          Rp = Rc;
          Rc = Rtmp;

          h *= 0.5;
        }

      /* did not converge - return best guess */
      *result = Rp[n - 1];

      return GSL_EMAXITER;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*solver/integration.c
*
* Function for numerical integration
* using quadrature and Monte Carlo
* methods.
*
* Author: Benjamin Vatter.
* email: benjaminvatterj@gmail.com
* date: 15 August 2015
*/

#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "integration.h"
#include "dbg.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>

typedef struct {
	int dims;
	int depth;
	double *lb;
	double *ub;
	double *real_x;
	double epsabs;
	double epsrel;
	size_t limit;
	gsl_integration_workspace * workspace;
	double * abserr;
	nquad_function *f;
} nquad_params;

void gsl_error (const char * reason, const char * file, int line, int gsl_errno)
{
  gsl_stream_printf ("NQUAD ERROR", file, line, reason);

  if (gsl_errno != GSL_EROUND) {
	  fflush (stderr);
	  fflush (stdout);
	  abort ();
  } else {
  	log_warn("NQUAD ERROR: rounding error detected. This result might be unreliable.");
  	fflush (stderr);
  	return;
  }
}

/*
SUbroutine for recursive multidimensional integration
*/
double nquad_f(double x, void * params)
{
	nquad_params *in_pars = (nquad_params *) params;
	double result;
	double abserr;
	int step = in_pars->depth;
	// Add value
	if (step>0) {
	 	in_pars->real_x[step-1] = x;
	}
	if (in_pars->dims == in_pars->depth){
		result = in_pars->f->function(in_pars->real_x, in_pars->dims, in_pars->f->params);
		return result;
	} else {
		// Consistency
		assert((step+1) <= in_pars->dims);
		// increase depth
		in_pars->depth += 1;
		gsl_integration_workspace *w = gsl_integration_workspace_alloc(in_pars->limit);
		gsl_function F = {.function = &nquad_f, .params = in_pars};
		gsl_integration_qags(&F, in_pars->lb[step], in_pars->ub[step], in_pars->epsabs, in_pars->epsrel, in_pars->limit, w,&result, &abserr);
		//gsl_integration_qag(&F, in_pars->lb[step], in_pars->ub[step], in_pars->epsabs, in_pars->epsrel, in_pars->limit, 2, w,&result, &abserr);
		in_pars->abserr[step] = abserr;
		gsl_integration_workspace_free(w);
		in_pars->depth -= 1;
		return result;
	}
}

/*
An extension of gsl_integration_qags
to n dimensions. The first parameter is
the number or dimensions. the signature of
f must be
f(double *x, void * params)

Parameters:
n - dimensions
f - function and parameters
lb - lower bounds to integrals
ub - upper bounds to integrals
epsabs - desired absolute error in each integral (only in nested-quadrature)
epsrel - desired relative error in each integral (only in nested-quadrature)
limit - max number of iterations
result - on out contains result
abserr - on out contains estimated error
*/
void nquad(int n, nquad_function *f,  double *lb, double *ub,
           double epsabs, double epsrel, size_t limit, double * result, double *abserr)
{
	assert(n>=1);
	nquad_params *new_params = malloc(sizeof(nquad_params));
	new_params->dims = n;
	new_params->depth = 0;
	new_params->lb = lb;
	new_params->ub = ub;
	new_params->epsabs = epsabs;
	new_params->epsrel = epsrel;
	new_params->limit = limit;
	new_params->abserr = malloc(sizeof(double) * n);
	new_params->f = f;
	new_params->real_x = malloc(sizeof(double) * n);

	// Set error handle
	gsl_error_handler_t *old_handler = gsl_set_error_handler(&gsl_error);

	// integrate
	result[0] = nquad_f(0.0, new_params);
	abserr[0] = 1.0;
	int i;
	double err = 0.0;
	for (i = 0; i<n; i++){
		err = (err < new_params->abserr[i])? new_params->abserr[i]: err;
	}
	*abserr = err;

	// release real x
	free(new_params->real_x);
	free(new_params->abserr);
	free(new_params);
	// recover error handle
	gsl_set_error_handler(old_handler);
}

/* Monte Carlo integration routine
* It has the same signature as the nquad function
* but it doesn't really use the epsrel value. It is
* done so to allow interchanging those two functions
* during runtime.
*
* This is based on the Monte Carlo plain routine found
* in GSL, written by Michael Booth
*/
void mcint(int n, nquad_function *f, double *lb, double *ub, double epsabs,
           double epsrel, size_t limit, double *result, double *abserr)
{
	double vol=1.0, m = 0, q = 0;
	double *x = malloc(sizeof(double)*n);
	size_t i, j;
	const gsl_rng_type *T;
	gsl_rng *r;
	gsl_rng_env_setup ();
	T = gsl_rng_default;
	r = gsl_rng_alloc (T);

	/* Sanity check */
	for(i=0; i<n; i++) {
		if(lb[i] >= ub[i]) {
			GSL_ERROR("Lower bound is higher than upper bound. Null feasible region.", GSL_EINVAL);
		}
		if(ub[i] - lb[i] > GSL_DBL_MAX) {
			GSL_ERROR("Range of integration is too large, please rescale",
                     GSL_EINVAL);
		}
		vol *= ub[i] - lb[i];
	}

	j=0;

	do
	{
		for (i=0; i < n; i++) {
			x[i] = lb[i] + gsl_rng_uniform_pos(r) * (ub[i] - lb[i]);
		}
		double fval = (f->function)(x, n, f->params);
		double d = fval - m;
		m += d / (n + 1.0);
		q += d * d * (n/(n +1.0));

		if (j < 2){ *abserr = GSL_POSINF;}
		else {
			*abserr = vol * sqrt(q / (limit * (limit - 1.0)));
		}
		j++;
	} while((*abserr > epsabs) && j < limit);

	*result = vol * m;
	if (*abserr > epsabs){
		log_warn("Monte Carlo integration max evaluation %d was insufficient"
		"to achieve the desired absolute error of %e", (int)limit, epsabs);
	}

	free(x);
	gsl_rng_free (r);
}


text:
#include "GeneralFunctions.h"
#include <math.h>
#include <stdio.h>
#include <iostream>
#include <time.h>
//
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include "gsl/gsl_cdf.h"
#include "gsl/gsl_randist.h"
#include "gsl/gsl_rng.h"

int AcceleratedGibbs (int maxK,int bias, int N, int D, int K, char *C,  int *R, double alpha, double s2B, double s2Y, gsl_matrix **Y, gsl_matrix *Z, int *nest, gsl_matrix *P, gsl_matrix *Pnon, gsl_matrix **lambda, gsl_matrix **lambdanon);
void SampleY (double missing, int N, int D, int K, char Cd,  int Rd, double fd, double mud, double wd, double s2Y,double s2u, double s2theta, gsl_matrix *X, gsl_matrix *Z, gsl_matrix *Yd,  gsl_matrix *Bd, gsl_vector *thetad, const gsl_rng *seed);
double Samples2Y (double missing, int N, int d, int K, char Cd,  int Rd, double fd, double mud, double wd, double s2u, double s2theta, gsl_matrix *X, gsl_matrix *Z, gsl_matrix *Yd,  gsl_matrix *Bd, gsl_vector *thetad, const gsl_rng *seed);
int IBPsampler_func (double missing, gsl_matrix *X, char *C, gsl_matrix *Z, gsl_matrix **B, gsl_vector **theta, int *R, double *f, double *mu,  double *w, int maxR, int bias, int N, int D, int K, double alpha, double s2B, double *s2Y,double s2u, int maxK,int Nsim);
int initialize_func (int N, int D, int maxK, double missing, gsl_matrix *X, char *C, gsl_matrix **B, gsl_vector **theta, int *R, double *f, double *mu,  double *w, double *s2Y);

max_stars_repo_name: ferjorosa/test-glfm

max_issues_repo_name: ferjorosa/test-glfm

max_forks_repo_name: ferjorosa/test-glfm


text:
// BSD 3-Clause License
//
// Copyright (c) 2018, Gregory Meyer
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// * Redistributions of source code must retain the above copyright
//   notice, this list of conditions and the following disclaimer.
//
// * Redistributions in binary form must reproduce the above copyright
//   notice, this list of conditions and the following disclaimer in
//   the documentation and/or other materials provided with the
//   distribution.
//
// * Neither the name of the copyright holder nor the names of its
//   contributors may be used to endorse or promote products derived
//   from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

#ifndef SWEEP_UTIL_TYPES_H
#define SWEEP_UTIL_TYPES_H

#include <cstdint>

#include <gsl/gsl_util>

namespace sweep::util {
inline namespace types {

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using i8 = std::int8_t;
using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;

using usize = std::size_t;
using isize = std::ptrdiff_t;

using f32 = float;
using f64 = double;

inline namespace literals {

constexpr u8 operator""_u8(unsigned long long x) noexcept {
    return gsl::narrow_cast<u8>(x);
}

constexpr u16 operator""_u16(unsigned long long x) noexcept {
    return gsl::narrow_cast<u16>(x);
}

constexpr u32 operator""_u32(unsigned long long x) noexcept {
    return gsl::narrow_cast<u32>(x);
}

constexpr u64 operator""_u64(unsigned long long x) noexcept {
    return gsl::narrow_cast<u64>(x);
}

constexpr i8 operator""_i8(unsigned long long x) noexcept {
    return gsl::narrow_cast<i8>(x);
}

constexpr i16 operator""_i16(unsigned long long x) noexcept {
    return gsl::narrow_cast<i16>(x);
}

constexpr i32 operator""_i32(unsigned long long x) noexcept {
    return gsl::narrow_cast<i32>(x);
}

constexpr i64 operator""_i64(unsigned long long x) noexcept {
    return gsl::narrow_cast<i64>(x);
}

constexpr usize operator""_usize(unsigned long long x) noexcept {
    return gsl::narrow_cast<usize>(x);
}

constexpr isize operator""_isize(unsigned long long x) noexcept {
    return gsl::narrow_cast<isize>(x);
}

constexpr f32 operator""_f32(long double x) noexcept {
    return gsl::narrow_cast<f32>(x);
}

constexpr f64 operator""_f64(long double x) noexcept {
    return gsl::narrow_cast<f64>(x);
}

} // inline namespace literals
} // inline namespace types
} // namespace sweep::util

#endif

max_stars_repo_name: Gregory-Meyer/minesweeper

max_issues_repo_name: Gregory-Meyer/minesweeper

max_forks_repo_name: Gregory-Meyer/minesweeper


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef NOMPI
#include <mpi.h>
#endif
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

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

/*! \file driftfac.c
 *  \brief compute loop-up tables for prefactors in cosmological integration
 */

static double logTimeBegin;
static double logTimeMax;

/*! This function computes look-up tables for factors needed in
 *  cosmological integrations. The (simple) integrations are carried out
 *  with the GSL library.  Separate factors are computed for the "drift",
 *  and the gravitational and hydrodynamical "kicks".  The lookup-table is
 *  used for reasons of speed.
 */
void init_drift_table(void)
{
#define WORKSIZE 100000
  int i;
  double result, abserr;
  gsl_function F;
  gsl_integration_workspace *workspace;

  logTimeBegin = log(All.TimeBegin);
  logTimeMax = log(All.TimeMax);

  workspace = gsl_integration_workspace_alloc(WORKSIZE);

  for(i = 0; i < DRIFT_TABLE_LENGTH; i++)
    {
      F.function = &drift_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), 0,
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      DriftTable[i] = result;

      F.function = &gravkick_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), 0,
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      GravKickTable[i] = result;

      F.function = &hydrokick_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), 0,
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      HydroKickTable[i] = result;
    }

  gsl_integration_workspace_free(workspace);
}

/*! This function integrates the cosmological prefactor for a drift step
 *  between time0 and time1. The value returned is * \f[ \int_{a_0}^{a_1}
 *  \frac{{\rm d}a}{H(a)} * \f]
 */
double get_drift_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * DriftTable[0];
  else
    df1 = DriftTable[i1 - 1] + (DriftTable[i1] - DriftTable[i1 - 1]) * (u1 - i1);

  u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * DriftTable[0];
  else
    df2 = DriftTable[i2 - 1] + (DriftTable[i2] - DriftTable[i2 - 1]) * (u2 - i2);

  return df2 - df1;
}

/*! This function integrates the cosmological prefactor for a kick step of
 *  the gravitational force.
 */
double get_gravkick_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * GravKickTable[0];
  else
    df1 = GravKickTable[i1 - 1] + (GravKickTable[i1] - GravKickTable[i1 - 1]) * (u1 - i1);

  u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * GravKickTable[0];
  else
    df2 = GravKickTable[i2 - 1] + (GravKickTable[i2] - GravKickTable[i2 - 1]) * (u2 - i2);

  return df2 - df1;
}

/*! This function integrates the cosmological prefactor for a kick step of
 *  the hydrodynamical force.
 */
double get_hydrokick_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * HydroKickTable[0];
  else
    df1 = HydroKickTable[i1 - 1] + (HydroKickTable[i1] - HydroKickTable[i1 - 1]) * (u1 - i1);

  u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * HydroKickTable[0];
  else
    df2 = HydroKickTable[i2 - 1] + (HydroKickTable[i2] - HydroKickTable[i2 - 1]) * (u2 - i2);

  return df2 - df1;
}

/*! Integration kernel for drift factor computation.
 */
double drift_integ(double a, void *param)
{
  double h;

  h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
  h = All.Hubble * sqrt(h);

  return 1 / (h * a * a * a);
}

/*! Integration kernel for gravitational kick factor computation.
 */
double gravkick_integ(double a, void *param)
{
  double h;

  h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
  h = All.Hubble * sqrt(h);

  return 1 / (h * a * a);
}

/*! Integration kernel for hydrodynamical kick factor computation.
 */
double hydrokick_integ(double a, void *param)
{
  double h;

  h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
  h = All.Hubble * sqrt(h);

  return 1 / (h * pow(a, 3 * GAMMA_MINUS1) * a);
}

double growthfactor_integ(double a, void *param)
{
  double s;

  s = All.Omega0 + (1 - All.Omega0 - All.OmegaLambda) * a + All.OmegaLambda * a * a * a;
  s = sqrt(s);

  return pow(sqrt(a) / s, 3);
}

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
#ifndef SIGSCANNERMEMORYDATA_HHHHH
#define SIGSCANNERMEMORYDATA_HHHHH

#include <vector>
#include <gsl/span>
#include <string>

namespace MPSig {
    class SigScannerMemoryData {
        gsl::span<const char> m_refData;
    public:
        explicit SigScannerMemoryData(gsl::span<const char> data);

        bool IsInRange(std::intptr_t addr) const;
        char* Deref(std::intptr_t addrWithOffset) const;
        std::pair<char*, bool> DerefTry(std::intptr_t addrWithOffset) const;

        gsl::span<const char> Get() const;
        std::intptr_t GetOffset() const;
    };
}

#endif


text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

#ifndef setup_h
#define setup_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include <petscksp.h>
#include <stdbool.h>
#include <string.h>

#if PETSC_VERSION_LT(3,14,0)
#  define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l,m,n) DMAddBoundary(a,b,c,e,h,i,j,k,f,g,m)
#elif PETSC_VERSION_LT(3,16,0)
#  define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l,m,n) DMAddBoundary(a,b,c,e,h,i,j,k,l,f,g,m)
#else
#  define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l,m,n) DMAddBoundary(a,b,c,d,f,g,h,i,j,k,l,m,n)
#endif

#ifndef PHYSICS_STRUCT
#define PHYSICS_STRUCT
typedef struct Physics_private *Physics;
struct Physics_private {
  CeedScalar   nu;      // Poisson's ratio
  CeedScalar   E;       // Young's Modulus
};
#endif

// -----------------------------------------------------------------------------
// Command Line Options
// -----------------------------------------------------------------------------
// Problem options
typedef enum {
  ELAS_LINEAR = 0, ELAS_SS_NH = 1, ELAS_FSInitial_NH1 = 2, ELAS_FSInitial_NH2 = 3,
  ELAS_FSCurrent_NH1 = 4, ELAS_FSCurrent_NH2 = 5
} problemType;
static const char *const problemTypes[] = {"Linear",
                                           "SS-NH",
                                           "FSInitial-NH1",
                                           "FSInitial-NH2",
                                           "FSCurrent-NH1",
                                           "FSCurrent-NH2",
                                           "problemType","ELAS_",0
                                          };
static const char *const problemTypesForDisp[] = {"Linear elasticity",
                                                  "Hyperelasticity small strain, Neo-Hookean",
                                                  "Hyperelasticity finite strain Initial config Neo-Hookean w/ dXref_dxinit, Grad(u) storage",
                                                  "Hyperelasticity finite strain Initial config Neo-Hookean w/ dXref_dxinit, Grad(u), C_inv, constant storage",
                                                  "Hyperelasticity finite strain Current config Neo-Hookean w/ dXref_dxinit, Grad(u) storage",
                                                  "Hyperelasticity finite strain Current config Neo-Hookean w/ dXref_dxcurr, tau, constant storage",
                                                 };

// Forcing function options
typedef enum {
  FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2
} forcingType;
static const char *const forcing_types[] = {"none",
                                            "constant",
                                            "mms",
                                            "forcingType","FORCE_",0
                                           };
static const char *const forcing_types_for_disp[] = {"None",
                                                     "Constant",
                                                     "Manufactured solution"
                                                    };

// Multigrid options
typedef enum {
  MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2
} multigridType;
static const char *const multigrid_types [] = {"logarithmic",
                                               "uniform",
                                               "none",
                                               "multigridType","MULTIGRID",0
                                              };
static const char *const multigrid_types_for_disp[] = {"P-multigrid, logarithmic coarsening",
                                                       "P-multigrind, uniform coarsening",
                                                       "No multigrid"
                                                      };

typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt,
                              PetscScalar *, void *);
// Note: These variables should be updated if additional boundary conditions
//         are added to boundary.c.
BCFunc BCMMS, BCZero, BCClamp;

// -----------------------------------------------------------------------------
// Structs
// -----------------------------------------------------------------------------
// Units
typedef struct Units_private *Units;
struct Units_private {
  // Fundamental units
  PetscScalar meter;
  PetscScalar kilogram;
  PetscScalar second;
  // Derived unit
  PetscScalar Pascal;
};

// Application context from user command line options
typedef struct AppCtx_private *AppCtx;
struct AppCtx_private {
  char          ceed_resource[PETSC_MAX_PATH_LEN];     // libCEED backend
  char          mesh_file[PETSC_MAX_PATH_LEN];         // exodusII mesh file
  char          output_dir[PETSC_MAX_PATH_LEN];
  PetscBool     test_mode;
  PetscBool     view_soln;
  PetscBool     view_final_soln;
  problemType   problem_choice;
  forcingType   forcing_choice;
  multigridType multigrid_choice;
  PetscScalar   nu_smoother;
  PetscInt      degree;
  PetscInt      q_extra;
  PetscInt      num_levels;
  PetscInt      *level_degrees;
  PetscInt      num_increments;                        // Number of steps
  PetscInt      bc_clamp_count;
  PetscInt      bc_clamp_faces[16];
  // [translation; 3] [rotation axis; 3] [rotation magnitude c_0, c_1]
  // The rotations are (c_0 + c_1 s) \pi, where s = x · axis
  PetscScalar   bc_clamp_max[16][8];
  PetscInt      bc_traction_count;
  PetscInt      bc_traction_faces[16];
  PetscScalar   bc_traction_vector[16][3];
  PetscScalar   forcing_vector[3];
};

// Problem specific data
// *INDENT-OFF*
typedef struct {
  CeedInt           q_data_size;
  CeedQFunctionUser setup_geo, apply, jacob, energy, diagnostic;
  const char        *setup_geo_loc, *apply_loc, *jacob_loc, *energy_loc,
                    *diagnostic_loc;
  CeedQuadMode      quad_mode;
} problemData;
// *INDENT-ON*

// Data specific to each problem option
extern problemData problem_options[6];

// Forcing function data
typedef struct {
  CeedQFunctionUser setup_forcing;
  const char        *setup_forcing_loc;
} forcingData;

extern forcingData forcing_options[3];

// Data for PETSc Matshell
typedef struct UserMult_private *UserMult;
struct UserMult_private {
  MPI_Comm        comm;
  DM              dm;
  Vec             X_loc, Y_loc, neumann_bcs;
  CeedVector      x_ceed, y_ceed;
  CeedOperator    op;
  CeedQFunction   qf;
  Ceed            ceed;
  PetscScalar     load_increment;
  CeedQFunctionContext ctx_phys, ctx_phys_smoother;
};

// Data for Jacobian setup routine
typedef struct FormJacobCtx_private *FormJacobCtx;
struct FormJacobCtx_private {
  UserMult     *jacob_ctx;
  PetscInt     num_levels;
  SNES         snes_coarse;
  Mat          *jacob_mat, jacob_mat_coarse;
  Vec          u_coarse;
};

// Data for PETSc Prolongation/Restriction Matshell
typedef struct UserMultProlongRestr_private *UserMultProlongRestr;
struct UserMultProlongRestr_private {
  MPI_Comm     comm;
  DM           dm_c, dm_f;
  Vec          loc_vec_c, loc_vec_f;
  CeedVector   ceed_vec_c, ceed_vec_f;
  CeedOperator op_prolong, op_restrict;
  Ceed         ceed;
};

// libCEED data struct for level
typedef struct CeedData_private *CeedData;
struct CeedData_private {
  Ceed                ceed;
  CeedBasis           basis_x, basis_u, basis_c_to_f, basis_energy,
                      basis_diagnostic;
  CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_qd_i,
                      elem_restr_gradu_i,
                      elem_restr_energy, elem_restr_diagnostic,
                      elem_restr_dXdx, elem_restr_tau,
                      elem_restr_C_inv, elem_restr_lam_log_J, elem_restr_qd_diagnostic_i;
  CeedQFunction       qf_apply, qf_jacob, qf_energy, qf_diagnostic;
  CeedOperator        op_apply, op_jacob, op_restrict, op_prolong, op_energy,
                      op_diagnostic;
  CeedVector          q_data, q_data_diagnostic, grad_u, x_ceed,
                      y_ceed, true_soln, dXdx, tau, C_inv, lam_log_J;
};

// Translate PetscMemType to CeedMemType
static inline CeedMemType MemTypeP2C(PetscMemType mem_type) {
  return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST;
}

// -----------------------------------------------------------------------------
// Process command line options
// -----------------------------------------------------------------------------
// Process general command line options
PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx);

// Process physics options
PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units);

// -----------------------------------------------------------------------------
// Setup DM
// -----------------------------------------------------------------------------
PetscErrorCode CreateBCLabel(DM dm, const char name[]);

// Create FE by degree
PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc,
                                     PetscBool is_simplex, const char prefix[],
                                     PetscInt order, PetscFE *fem);

// Read mesh and distribute DM in parallel
PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx app_ctx, DM *dm);

// Setup DM with FE space of appropriate degree
PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order,
                               PetscBool boundary, PetscInt num_comp_u);

// -----------------------------------------------------------------------------
// libCEED Functions
// -----------------------------------------------------------------------------
// Destroy libCEED objects
PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data);

// Utility function - essential BC dofs are encoded in closure indices as -(i+1)
PetscInt Involute(PetscInt i);

// Utility function to create local CEED restriction from DMPlex
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P,
    CeedInt height, DMLabel domain_label, CeedInt value,
    CeedElemRestriction *elem_restr);

// Utility function to get Ceed Restriction for each domain
PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height,
                                       DMLabel domain_label, PetscInt value, CeedInt P, CeedInt Q, CeedInt q_data_size,
                                       CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x,
                                       CeedElemRestriction *elem_restr_qd_i);

// Set up libCEED for a given degree
PetscErrorCode SetupLibceedFineLevel(DM dm, DM dm_energy, DM dm_diagnostic,
                                     Ceed ceed, AppCtx app_ctx,
                                     CeedQFunctionContext phys_ctx,
                                     CeedData *data, PetscInt fine_level,
                                     PetscInt num_comp_u, PetscInt U_g_size,
                                     PetscInt u_loc_size, CeedVector force_ceed,
                                     CeedVector neumann_ceed);

// Set up libCEED multigrid level for a given degree
PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx,
                                 CeedData *data, PetscInt level,
                                 PetscInt num_comp_u, PetscInt U_g_size,
                                 PetscInt u_loc_size, CeedVector fine_mult);

// Setup context data for Jacobian evaluation
PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V,
                                Vec V_loc, CeedData ceed_data, Ceed ceed,
                                CeedQFunctionContext ctx_phys,
                                CeedQFunctionContext ctx_phys_smoother,
                                UserMult jacobian_ctx);

// Setup context data for prolongation and restriction operators
PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c,
                                       DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f,
                                       CeedData ceed_data_c, CeedData ceed_data_f,
                                       Ceed ceed,
                                       UserMultProlongRestr prolong_restr_ctx);

// -----------------------------------------------------------------------------
// Jacobian setup
// -----------------------------------------------------------------------------
PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx);

// -----------------------------------------------------------------------------
// Solution output
// -----------------------------------------------------------------------------
PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U,
                            PetscInt increment,
                            PetscScalar load_increment);

PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dm_U,
                                        UserMult user, AppCtx app_ctx, Vec U,
                                        CeedElemRestriction elem_restr_diagnostic);

// -----------------------------------------------------------------------------
// libCEED Operators for MatShell
// -----------------------------------------------------------------------------
// This function uses libCEED to compute the local action of an operator
PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user);

// This function uses libCEED to compute the non-linear residual
PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to apply the Jacobian for assembly via a SNES
PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to compute the action of the Jacobian
PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the prolongation operator
PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the restriction operator
PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y);

// This function returns the computed diagonal of the operator
PetscErrorCode GetDiag_Ceed(Mat A, Vec D);

// This function calculates the strain energy in the final solution
PetscErrorCode ComputeStrainEnergy(DM dm_energy, UserMult user,
                                   CeedOperator op_energy, Vec X,
                                   PetscReal *energy);

// -----------------------------------------------------------------------------
// Boundary Functions
// -----------------------------------------------------------------------------
// Note: If additional boundary conditions are added, an update is needed in
//         elasticity.h for the boundaryOptions variable.

// BCMMS - boundary function
// Values on all points of the mesh is set based on given solution below
// for u[0], u[1], u[2]
PetscErrorCode BCMMS(PetscInt dim, PetscReal load_increment,
                     const PetscReal coords[], PetscInt num_comp_u,
                     PetscScalar *u, void *ctx);

// BCClamp - fix boundary values with affine transformation at fraction of load
//   increment
PetscErrorCode BCClamp(PetscInt dim, PetscReal load_increment,
                       const PetscReal coords[], PetscInt num_comp_u,
                       PetscScalar *u, void *ctx);

#endif //setup_h

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED