text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * 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.
 *
 */
/*
 *
 *
 *       Filename:  entropy.c
 *
 *    Description:  Calculate the entropy production for a given state, generator
 *    			and equilibrium state
 *
 *        Version:  1.0
 *        Created:  05/05/2014 20:53:01
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  entropy_production
 *  Description:
 *
 */
double entropy_production ( const gsl_vector* rho, const gsl_vector* rhoeq,
		const gsl_matrix* L  )
{
	/* l1, l2, l3 */
	double l[3] ; unsigned int i, j ;
	for ( i = 1 ; i < 3 ; i++ )
	{
		l[i] = 0 ;
		for ( j = 0 ; j < 3 ; j++ )
			l[i] += gsl_matrix_get(L,i,j)*VECTOR(rho,j) ;
	}

	/* L[rho] */
	double Lr = 0 ;
	for ( i = 1 ; i < 3 ; i++ )
		Lr += l[i]*VECTOR(rho,i) ;

	/* L[rhoeq] */
	double Leq = 0 ;
	for ( i = 1 ; i < 3 ; i++ )
		Leq += l[i]*VECTOR(rhoeq,i) ;

	/* r , req */
	double r, req ;
	r = req = 0 ;

	r = gsl_hypot3(VECTOR(rho,1),VECTOR(rho,2),VECTOR(rho,3)) ;
	req = gsl_hypot3(VECTOR(rhoeq,1),VECTOR(rhoeq,2),
		VECTOR(rhoeq,3)) ;

	/* internal entropy s */
	double s ;
	if ( r < 1 && req < 1 )
		s = -(gsl_sf_log((1+r)/(1-r))*Lr/r -gsl_sf_log((1+req)/(1-req))*Leq/req);
	else if ( r < 1 && req >= 0 )
		s = -gsl_sf_log((1+r)/(1-r))*Lr/r ;
	else
		s = 0 ;

	return s;
}		/* -----  end of function entropy_production  ----- */

/*
 *      FUNCTION
 *         Name:  entropy_of_state
 *  Description:  Calculate the Von Neumann entropy of state 'rho'
 *
 */
double entropy_of_state ( const gsl_vector* rho )
{
	double entr = 0 ;

	/* Finding the eigenvalues */
	gsl_eigen_herm_workspace* rho_ei = gsl_eigen_herm_alloc(2);
	gsl_matrix_complex* dens = gsl_matrix_complex_calloc (2,2);
	gsl_matrix_complex_set (dens, 0, 0, gsl_complex_rect(1+VECTOR(rho, 3),0));
	gsl_matrix_complex_set (dens,0,1,gsl_complex_rect(VECTOR(rho,1),-VECTOR(rho,2)));
	gsl_matrix_complex_set (dens,1,0,gsl_complex_rect(VECTOR(rho,1),VECTOR(rho,2)));
	gsl_matrix_complex_set (dens,1,1,gsl_complex_rect(1-VECTOR(rho,3),0));
	gsl_matrix_complex_scale (dens, gsl_complex_rect(0.5,0));
	gsl_vector* eigenvalues = gsl_vector_calloc(2) ;
	gsl_eigen_herm (dens, eigenvalues, rho_ei) ;

	/* Calculating entropy */
	double norm = gsl_hypot3( VECTOR(rho,1), VECTOR(rho,2),	VECTOR(rho,3) ) ;
	if ( gsl_fcmp(norm, 1, 1e-9) > 0 )
		entr = 0 ;
	else
		entr = - (VECTOR(eigenvalues,0)*gsl_sf_log(VECTOR(eigenvalues,0)) +
			VECTOR(eigenvalues,1)*gsl_sf_log(VECTOR(eigenvalues,1))) ;

	return (entr);
}		/* -----  end of function entropy_of_state  ----- */


text:
/****************************************************************
 *
 * vpreprocess: FreqFilter.c
 *
 * Copyright (C) Max Planck Institute
 * for Human Cognitive and Brain Sciences, Leipzig
 *
 * <lipsia@cbs.mpg.de>
 *
 * 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 * $Id: FreqFilter.c 3181 2008-04-01 15:19:44Z karstenm $
 *
 *****************************************************************/

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

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

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

extern gsl_matrix *DataMap(VImage *src,VImage mask,size_t nvox,int);

gsl_matrix *DataMap(VImage *src,VImage mask,size_t nvox,int tail)
{
  int j,b,r,c;
  int nrows=0,ncols=0,nt0=0;
  double tiny=1.0e-6;

  int nslices = VImageNBands(mask);
  VImageDimensions(src,nslices,&nt0,&nrows,&ncols);

  int nt = nt0 + 2*tail;
  gsl_matrix *A = gsl_matrix_calloc(nvox,nt);

  size_t n=0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	double u = VGetPixel(src[b],0,r,c);
	if (fabs(u) < tiny) continue;

	for (j=0; j<tail; j++) {
	  u = (double)VGetPixel(src[b],tail-j,r,c);
	  gsl_matrix_set(A,n,j,u);
	}
	for (j=tail; j<nt-tail; j++) {
	  u = (double)VGetPixel(src[b],j-tail,r,c);
	  gsl_matrix_set(A,n,j,u);
	}
	for (j=nt-tail; j<nt; j++) {
	  u = (double)VGetPixel(src[b],2*nt-3*tail-2-j,r,c);
	  gsl_matrix_set(A,n,j,u);
	}
	n++;
      }
    }
  }
  return A;
}

void VFreqFilter(VAttrList list,VFloat high,VFloat low,VBoolean stop,VFloat sharp)
{
  VString str1=NULL;
  double *in=NULL;
  fftw_complex *out;
  fftw_plan p1,p2;
  int i,j,nc,b,r,c;
  float freq=0,alpha=0;
  double *highp=NULL, *lowp=NULL;
  double x_high, x_low, x;
  double u=0,v=0,tiny=1.0e-6;

  /* dialog messages */
  str1 = (VString)VMalloc(sizeof(char)*4);
  if (stop) strcpy(str1, "stop");
  else strcpy(str1, "pass");
  if (low>0 && high>0)
    fprintf(stderr," band %s filter: 1/%.1f Hz and 1/%.1f Hz\n",str1,high,low);
  else {
    if (high>0) fprintf(stderr," high %s filter: 1/%.1f Hz\n",str1,high);
    if (low>0)  fprintf(stderr," low %s filter: 1/%.1f Hz\n",str1,low);
  }

  /* read functional data */
  int nrows=0,ncols=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&nt,&nrows,&ncols);

  /* read header info */
  double *D = (double *)VCalloc(8,sizeof(double));
  for (i=0; i<8; i++) D[i] = 1.0;
  VAttrList geolist = VGetGeoInfo(list);
  if (geolist != NULL) D = VGetGeoPixdim(geolist,D);
  double tr = D[4];

  VFloat xtr=-1;
  if (VGetAttr(VImageAttrList(src[0]),"repetition_time",NULL,VFloatRepn,&xtr) == VAttrFound) {
    tr = (double)xtr;
  }
  if (tr > 100) tr /= 1000.0;  /* convert to seconds */
  if (tr < 0.01) VError("FreqFilter: implausible TR (%f seconds)",tr);
  fprintf(stderr," image dimensions: %d x %d x %d,  nt= %d,  TR= %.3f secs\n",nslices,nrows,ncols,nt,tr);

  /* alloc memory */
  double nx = (double)nt;
  int tail = 50;
  if (nt<tail) tail=nt-2;
  nt  += 2 * tail;
  nc  = (nt / 2) + 1;
  in  = (double *)fftw_malloc(sizeof(double) * nt);
  out = fftw_malloc (sizeof (fftw_complex ) * nc);
  for (i=0; i<nt; i++) in[i] = 0;

  /* make plans */
  p1 = fftw_plan_dft_r2c_1d (nt,in,out,FFTW_ESTIMATE);
  p2 = fftw_plan_dft_c2r_1d (nt,out,in,FFTW_ESTIMATE);

  /* repetition time */
  alpha = (double)nt * tr;

  /* filter function */
  if (sharp > 0) {
    highp = (double *)malloc(sizeof(double) * nc);
    lowp  = (double *)malloc(sizeof(double) * nc);
    for (i=1; i <nc; i++) {
      highp[i] = 1.0 / (1.0 +  exp( (alpha/high -(double)i)*sharp )   );
      lowp[i]  = 1.0 / (1.0 +  exp( ((double)i - alpha/low)*sharp )   );
    }
  }

  /* create ROI mask */
  VImage mask = VCreateImage(nslices,nrows,ncols,VBitRepn);
  VFillImage(mask,VAllBands,0);
  size_t nvox=0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	double u = VGetPixel(src[b],0,r,c);
	if (fabs(u) < tiny) continue;
	VPixel(mask,b,r,c,VBit) = 1;
	nvox++;
      }
    }
  }

  /* map to double values */
  gsl_matrix *A = DataMap(src,mask,nvox,tail);
  double zmax = gsl_matrix_max(A);
  double zmin = gsl_matrix_min(A);

  /* main process */
  size_t n=0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	if (VPixel(mask,b,r,c,VBit) == 0) continue;
	for (i=0; i<nt; i++) {
	  in[i] = gsl_matrix_get(A,n,i);
	}

	/*  fft */
	fftw_execute(p1);

	/* remove specified frequencies */
	for (i=1; i <nc; i++) {

	  if (sharp > 0) {

	    /* highpass */
	    if (high > 0) x_high = highp[i];
	    else x_high = 1.0;

	    /* lowpass */
	    if (low > 0) x_low = lowp[i];
	    else x_low = 1.0;

	    x = x_high + x_low - 1.0;
	    if (stop) x = fabs(1.0-x);

	    out[i][0] *= x;
	    out[i][1] *= x;

	  }
	  else {

	    /* hard thresholding */
	    freq = 1.0/(double)i * alpha;     /* 1/frequency */

	    if ((!stop && (freq < low || (freq > high && high>0)))
		|| (stop && !(freq < low || (freq > high && high>0))))
	      out[i][0] = out[i][1] = 0;
	  }
	}

	/* inverse fft */
	fftw_execute(p2);
	for (i=0; i<nt; i++) {
	  gsl_matrix_set(A,n,i,in[i]/nx);
	}
	n++;
      }
    }
  }

  /* rescale */
  VRepnKind repn = VPixelRepn(src[0]);
  double smax = VPixelMaxValue(src[0]);
  double smin = VPixelMinValue(src[0]);

  double ymax = 32000;
  if (repn == VFloatRepn || repn == VDoubleRepn) ymax = 1000.0;

  for (b=0; b<nslices; b++) {
    VFillImage(src[b],VAllBands,0);
  }
  zmax = gsl_matrix_max(A);
  zmin = gsl_matrix_min(A);

  /* write output */
  n=0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (VPixel(mask,b,r,c,VBit) == 0) continue;

	/* map to output image */
	for (j=tail; j<nt-tail; j++) {
	  u = gsl_matrix_get(A,n,j);
	  if (fabs(u) > 0.0001) u = (u-zmin)/(zmax-zmin);
	  else u = 0;
	  if (u > 1.0) u = 1.0;
	  if (u < 0) u = 0;

	  v = u*ymax;
	  if (v > smax) v = smax;
	  if (v < smin) v = smin;
	  VSetPixel(src[b],j-tail,r,c,v);
	}
	n++;
      }
    }
  }
  VDestroyImage(mask);
  gsl_matrix_free(A);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2017 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <cblas.h>
#include <lapacke.h>
#include <rokko/cmatrix.h>

int imax(int x, int y) { return (x > y) ? x : y; }

int main(int argc, char** argv) {
  int n = 5;
  int i, j, info;
  double norm;
  double *w;
  double **a, **u, **t;

  if (argc > 1) n = atoi(argv[1]);

  /* generate matrix */
  a = alloc_dmatrix(n, n);
  for (j = 0; j < n; ++j) {
    for (i = 0; i < n; ++i) {
      mat_elem(a, i, j) = n - imax(i, j);
    }
  }
  printf("Matrix A: ");
  fprint_dmatrix(stdout, n, n, a);

  /* diagonalization */
  u = alloc_dmatrix(n, n);
  cblas_dcopy(n * n, mat_ptr(a), 1, mat_ptr(u), 1);
  w = alloc_dvector(n);
  info = LAPACKE_dsyev(LAPACK_COL_MAJOR, 'V', 'U', n, mat_ptr(u), n, vec_ptr(w));
  if (info != 0) {
    fprintf(stderr, "Error: dsyev fails with error code %d\n", info);
    exit(255);
  }
  printf("Eigenvalues: ");
  fprint_dvector(stdout, n, w);
  printf("Eigenvectors: ");
  fprint_dmatrix(stdout, n, n, u);

  /* orthogonality check */
  t = alloc_dmatrix(n, n);
  cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n, n, n, 1, mat_ptr(u), n,
              mat_ptr(u), n, 0, mat_ptr(t), n);
  for (i = 0; i < n; ++i) mat_elem(t, i, i) -= 1;
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', n, n, mat_ptr(t), n);
  printf("|| U^t U - I || = %e\n", norm);
  if (norm > 1e-10) {
    fprintf(stderr, "Error: orthogonality check\n");
    exit(255);
  }

  /* eigenvalue check */
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1, mat_ptr(a), n,
              mat_ptr(u), n, 0, mat_ptr(t), n);
  cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n, n, n, 1, mat_ptr(u), n,
              mat_ptr(t), n, 0, mat_ptr(a), n);
  for (i = 0; i < n; ++i) mat_elem(a, i, i) -= w[i];
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', n, n, mat_ptr(a), n);
  printf("|| U^t A U - diag(w) || = %e\n", norm);
  if (norm > 1e-10) {
    fprintf(stderr, "Error: eigenvalue check\n");
    exit(255);
  }

  free_dmatrix(a);
  free_dmatrix(u);
  free_dvector(w);
  free_dmatrix(t);
  return 0;
}

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf.h>

// Compute the psi polygamma function
int main (int argc, char **argv)
{
  if ( argc != 3 )
    {
      fprintf(stderr,"Incorrect arguments\n");
      return 1;
    }

  double x1 = atof(argv[1]);
  double x2 = atof(argv[2]);

  printf("%g\n",gsl_sf_psi_n(x1,x2));

  return 0;
}

max_stars_repo_name: sfarrens/cosmostat

max_issues_repo_name: sfarrens/cosmostat

max_forks_repo_name: sfarrens/cosmostat


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <stdbool.h>
#include "parmt_config.h"
#include "compearth.h"
#ifdef PARMT_USE_INTEL
#include <mkl.h>
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#include <ipps.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "iscl/array/array.h"
#include "iscl/linalg/linalg.h"
#include "iscl/memory/memory.h"
#include "iscl/statistics/statistics.h"
#include "iscl/sorting/sorting.h"

static double median_sorted_array(int n, const double *__restrict x);
static int setG64f(const int npts,
                   const double *__restrict__ Gxx,
                   const double *__restrict__ Gyy,
                   const double *__restrict__ Gzz,
                   const double *__restrict__ Gxy,
                   const double *__restrict__ Gxz,
                   const double *__restrict__ Gyz,
                   double *__restrict__ G);

#define LDG 8

int parmt_computeL1Magnitude64f(const int ldm,
                                const int nobs,
                                const int nmtSolve,
                                const int maxit, const double eps,
                                const double tol,
                                const int *__restrict__ signalPtr,
                                const int *__restrict__ lags,
                                const int *__restrict__ mtPtr,
                                const double *__restrict__ Gxx,
                                const double *__restrict__ Gyy,
                                const double *__restrict__ Gzz,
                                const double *__restrict__ Gxy,
                                const double *__restrict__ Gxz,
                                const double *__restrict__ Gyz,
                                const double *__restrict__ mts,
                                const double *__restrict__ d,
                                double *__restrict__ mags)
{
    const char *fcnm = "parmt_computeL1Magnitude64f\0";
    double *G, *est, *obs, *wts, xmag;
    double m6[8] __attribute__((aligned(64)));
    int i, i1, i2, ierr, imt, jmt, k, maxlag, nptsAll, nptsPad;
    bool luseLag;
    const double p = 1.0;
    const double one = 1.0;
    const double zero = 0.0;
    nptsAll = signalPtr[nobs];
    if (nobs < 1)
    {
        printf("%s: Error no observations\n", fcnm);
        return -1;
    }
    // Get the max number of lags
    maxlag = 0;
    luseLag = false;
    if (lags != NULL)
    {
        luseLag = true;
        maxlag = lags[array_absArgmax32i(nobs, lags, &ierr)];
    }
    if (maxlag == 0){luseLag = false;}
    nptsPad = nptsAll + nobs*maxlag;
    G = memory_calloc64f(nptsPad*LDG);
    obs = memory_calloc64f(nptsPad);
    est = memory_calloc64f(nptsPad);
    wts = memory_calloc64f(nptsPad);
    // Insert the Green's functions and data into the moment tensor matrix
    if (!luseLag)
    {
        for (k=0; k<nobs; k++)
        {
            i1 = signalPtr[k];
            i2 = signalPtr[k+1];
            for (i=i1; i<i2; i++)
            {
                G[i*LDG+0] = Gxx[i];
                G[i*LDG+1] = Gyy[i];
                G[i*LDG+2] = Gzz[i];
                G[i*LDG+3] = Gxy[i];
                G[i*LDG+4] = Gxz[i];
                G[i*LDG+5] = Gyz[i];
                obs[i] = d[i];
            }
        }
    }
    else
    {

    }
    for (imt=0; imt<nmtSolve; imt++)
    {
        // Extract the moment tensor
        jmt = mtPtr[imt];
        for (i=0; i<6; i++)
        {
            m6[0] = mts[ldm*jmt+0];
        }
        // Remove the scalar moment
        compearth_CMT2m0(1, 1, m6, &xmag);
        cblas_dscal(6, 1.0/xmag, m6, 1);
        // Create the linear model [G*m]*M0 = [est]*{M0} = {obs}
        cblas_dgemv(CblasRowMajor, CblasNoTrans,
                    nptsAll, 6, one, G, LDG, m6, 1, zero, est, 1);
        // Solve for magnitude in the L1 norm
        ierr = linalg_irls64f_work(nptsAll, 1,
                                   maxit, p,
                                   eps, tol,
                                   est, obs, &xmag, wts);
        mags[imt] = xmag;
 double xmagMw;
compearth_m02mw(1, CE_KANAMORI_1978, &xmag, &xmagMw);
printf("%f\n", xmagMw);
    }
    memory_free64f(&wts);
    memory_free64f(&G);
    memory_free64f(&obs);
    memory_free64f(&est);
    return 0;
}

int parmt_computeL1StationMagnitude64f(const int ldm, const int nmt,
                                       const int npts, const int nblock,
                                       const double *__restrict__ Gxx,
                                       const double *__restrict__ Gyy,
                                       const double *__restrict__ Gzz,
                                       const double *__restrict__ Gxy,
                                       const double *__restrict__ Gxz,
                                       const double *__restrict__ Gyz,
                                       const double *__restrict__ mts,
                                       const double *__restrict__ d,
                                       double *__restrict__ mags)
{
    const char *fcnm = "parmt_computeL1StationMagnitude64f\0";
    double *G, *M, *R, *Dmat, *res, *resSort, xmag, xopt;
    const double one = 1.0;
    const double zero = 0.0;
    const double negOne =-one;
    int *perm;
    bool lsorted;
    int i, ic, idx, ierr, imt, jdx, jmt, mblock, Mrows, Ncols;
    const int Kcols = 6; // Number of columns of G
    Mrows = npts;
    mblock = nblock;// + computePadding6f(npts);
    G = memory_calloc64f(LDG*npts);
    R = memory_calloc64f(mblock*npts);
    M = memory_calloc64f(mblock*8);
    Dmat = memory_calloc64f(mblock*npts);
    res = memory_calloc64f(npts);
    resSort = memory_calloc64f(npts);
    perm = memory_calloc32i(npts);
    for (i=0; i<npts; i++){perm[i] = i;}
        // Set the observations
        for (i=0; i<npts; i++)
        {
            for (ic=0; ic<nblock; ic++)
            {
                Dmat[mblock*i+ic] = d[i];
            }
        }
#ifdef __INTEL_COMPILER
    __assume_aligned(G, 64);
#endif
    ierr = setG64f(npts, Gxx, Gyy, Gzz, Gxy, Gxz, Gyz, G);
    if (ierr != 0)
    {
        printf("%s: Failed to set G\n", fcnm);
        return -1;
    }
int nsorted = 0;
    for (jmt=0; jmt<nmt; jmt=jmt+nblock)
    {
        Ncols = MIN(nblock, nmt - jmt); // Number of columns of M
        // Set the observations
        cblas_dcopy(npts*mblock, Dmat, 1, R, 1);
        for (i=0; i<6; i++)
        {
            for (ic=0; ic<Ncols; ic++)
            {
                imt = jmt + ic;
                idx = ldm*imt + i;
                jdx = mblock*i;
                M[jdx+ic] = mts[idx];///xmag;
            }
        }
        // Compute R = GM - D = GM - R
        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                    Mrows, Ncols, Kcols, one, G, LDG,
                    M, mblock, negOne, R, mblock);
        //lsorted = false;
        //for (i=0; i<npts; i++){perm[i] = i;}
        // Sort the columns of the residual matrix
        for (ic=0; ic<Ncols; ic++)
        {
            imt = jmt + ic;
            // Compute the absolute value of the residual
            for (i=0; i<npts; i++)
            {
                res[i] = fabs(R[ic+i*mblock]);
            }
            //cblas_dcopy(npts, &R[ic], mblock, res, 1);
            // Apply the old permutation
            sorting_applyPermutation64f_work(npts, perm, res, resSort);
            lsorted = sorting_issorted64f(npts, resSort, SORT_ASCENDING, &ierr);
            if (!lsorted)
            {
                sorting_argsort64f_work(npts, res, SORT_ASCENDING, perm);
                sorting_applyPermutation64f_work(npts, perm, res, resSort);
            }
            else
            {
                nsorted = nsorted + 1;
            }
            // Compute the median
            xopt = median_sorted_array(npts, resSort);
            //compearth_CMT2mw(1, 1, &mts[8*imt], &xmag);
            compearth_CMT2m0(1, 1, &mts[8*imt], &xmag);
            xopt = xopt*xmag;
            compearth_m02mw(1, CE_KANAMORI_1978, &xopt, &mags[imt]);
        }
    }
printf("%d %d\n", nsorted, nmt);
    memory_free32i(&perm);
    memory_free64f(&G);
    memory_free64f(&R);
    memory_free64f(&M);
    memory_free64f(&Dmat);
    memory_free64f(&res);
    memory_free64f(&resSort);
    return 0;
}

static double median_sorted_array(int n, const double *__restrict x)
{
    double xmed;
    int n2;
    const double half = 0.5;
    n2 = n/2;
    // Even -> average middle two elements
    if (fmod(n, 2) == 0)
    {
        xmed = half*(x[n2-1] + x[n2]);
    }
    // Median is middle element
    else
    {
        xmed = x[n2];
    }
    return xmed;
}

static int setG64f(const int npts,
                   const double *__restrict__ Gxx,
                   const double *__restrict__ Gyy,
                   const double *__restrict__ Gzz,
                   const double *__restrict__ Gxy,
                   const double *__restrict__ Gxz,
                   const double *__restrict__ Gyz,
                   double *__restrict__ G)
{
    const char *fcnm = "setG64f\0";
    int i;
    bool lalign;
    if (Gxx == NULL || Gyy == NULL || Gzz == NULL || Gxy == NULL ||
        Gxz == NULL || G == NULL || npts < 1)

    {
        if (Gxx == NULL){printf("%s: Error Gxx is NULL\n", fcnm);}
        if (Gyy == NULL){printf("%s: Error Gyy is NULL\n", fcnm);}
        if (Gzz == NULL){printf("%s: Error Gzz is NULL\n", fcnm);}
        if (Gxy == NULL){printf("%s: Error Gxy is NULL\n", fcnm);}
        if (Gxz == NULL){printf("%s: Error Gxz is NULL\n", fcnm);}
        if (Gyz == NULL){printf("%s: ERror Gyz is NULL\n", fcnm);}
        if (npts < 1)
        {
            printf("%s: No points in Green's functions\n", fcnm);
        }
        return -1;
    }
    lalign = true;
    if (memory_isAligned(Gxx, 64) != 1 || memory_isAligned(Gyy, 64) != 1 ||
        memory_isAligned(Gzz, 64) != 1 || memory_isAligned(Gxy, 64) != 1 ||
        memory_isAligned(Gxz, 64) != 1 || memory_isAligned(Gyz, 64) != 1 ||
        memory_isAligned(G, 64) != 1)
    {
        lalign = false;
    }
    if (lalign)
    {
        #pragma omp simd aligned(G, Gxx, Gyy, Gzz, Gxy, Gxz, Gyz: 64)
        for (i=0; i<npts; i++)
        {
            G[LDG*i+0] = Gxx[i];
            G[LDG*i+1] = Gyy[i];
            G[LDG*i+2] = Gzz[i];
            G[LDG*i+3] = Gxy[i];
            G[LDG*i+4] = Gxz[i];
            G[LDG*i+5] = Gyz[i];
        }
    }
    else
    {
        #pragma omp simd
        for (i=0; i<npts; i++)
        {
            G[LDG*i+0] = Gxx[i];
            G[LDG*i+1] = Gyy[i];
            G[LDG*i+2] = Gzz[i];
            G[LDG*i+3] = Gxy[i];
            G[LDG*i+4] = Gxz[i];
            G[LDG*i+5] = Gyz[i];
        }
    }
    return 0;
}


text:
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include "time_utility_c.h"

/** INS Variables **/
/* Matrix */
gsl_matrix *WEII;
gsl_matrix *TBIC;
gsl_matrix *TDCI;
gsl_matrix *TEIC;
gsl_matrix *TBDC;
gsl_matrix *TBICI;
gsl_matrix *TLI;

/* Vector */
gsl_vector *EVBI;
gsl_vector *EVBID;
gsl_vector *ESBI;
gsl_vector *ESBID;
gsl_vector *RICI;
gsl_vector *RICID;
gsl_vector *TBIC_Q;
gsl_vector *TBIDC_Q;
gsl_vector *SBIIC;
gsl_vector *VBIIC;
gsl_vector *SBEEC;
gsl_vector *VBEEC;
gsl_vector *WBICI;
gsl_vector *EGRAVI;
gsl_vector *VBECD;
gsl_vector *INS_I_ATT_ERR;
gsl_vector *TESTV;
gsl_vector *TMP_old;
gsl_vector *VBIIC_old;
gsl_vector *POS_ERR;
gsl_vector *GRAVGI;
gsl_vector *TBDCQ;
gsl_vector *VBIIC_old_old;
gsl_vector *PHI_C;
gsl_vector *DELTA_VEL_C;
gsl_vector *PHI_LOW_C;
gsl_vector *PHI_HIGH_C;

/* Double */
double dbic;
double dvbec;
double alphacx;
double betacx;
double thtvdcx;
double psivdcx;
double alppcx;
double phipcx;
double loncx;
double latcx;
double altc;
double phibdcx;
double thtbdcx;
double psibdcx;
double ins_pos_err;
double ins_vel_err;
double ins_tilt_err;
double ins_pose_err;
double ins_vele_err;
double ins_phi_err;
double ins_tht_err;
double ins_psi_err;

/* Unsigned int */
unsigned int gpsupdate;
unsigned int liftoff;
unsigned int ideal;

GPS_TIME gpstime;
UTC_TIME utctime;
/*******************/

/* Control variables */
double fmasse;
double mdot;
double fmass0;
double xcg_0;
double xcg_1;
double isp;
double xcg;
double thrust;
double mass_ratio;
double reference_point;
double d;
double theta_a_cmd;
double theta_b_cmd;
double theta_c_cmd;
double theta_d_cmd;
double lx;
// pitch
double thterror;
double perrori;
double perrorp;
double perror_old;
double pitchiout_old;
double pN;
double pdout_old;
double kpp;
double kpi;
double kpd;
double kppp;
double pitchcmd;
double pitchcmd_new;
double pitchcmd_old;
double pitchcmd_out_old;
// roll
double rollerror;
double rerrori;
double rerrorp;
double rerror_old;
double rolliout_old;
double rN;
double rdout_old;
double krp;
double kri;
double krd;
double krpp;
double rollcmd;
// yaw
double yawerror;
double yerrori;
double yerrorp;
double yerror_old;
double yawiout_old;
double yN;
double ydout_old;
double kyp;
double kyi;
double kyd;
double kypp;
double yawcmd;
// aoa
double aoaerror;
double aoaerrori;
double aoaerrorp;
double aoaerror_old;
double aoaiout_old;
double aoaN;
double aoadout_old;
double kaoap;
double kaoai;
double kaoad;
double kaoapp;
double aoaerror_smooth_old;
double aoaerror_new;
double aoaerror_out_old;
double aoacmd;

gsl_vector *IBBB0;
gsl_vector *IBBB1;
gsl_vector *IBBB2;
gsl_vector *CONTROLCMD;
gsl_vector *CMDQ;
gsl_vector *TCMDQ;
gsl_vector *delta_euler;
gsl_vector *euler;
gsl_vector *WBICBT;
/*********************/


text:
#pragma once

#include "text_parser.h"
#include "text_renderer.h"
#include "terminal_io.h"
#include "geometry.h"

#include <string>
#include <vector>

#include <gsl/span>

namespace terminal_editor {

enum class Color {
    Black	        = 30,
    Red	            = 31,
    Green	        = 32,
    Yellow	        = 33,
    Blue	        = 34,
    Magenta	        = 35,
    Cyan	        = 36,
    White	        = 37,
    Bright_Black	= 90,
    Bright_Red	    = 91,
    Bright_Green	= 92,
    Bright_Yellow	= 93,
    Bright_Blue	    = 94,
    Bright_Magenta	= 95,
    Bright_Cyan	    = 96,
    Bright_White	= 97,
};

enum class Style {
    Bold        = 1,
    Normal      = 22,
};

struct Attributes {
    Color fgColor;
    Color bgColor;
    Style style;
};

class ScreenBuffer;

class ScreenCanvas {
    ScreenBuffer& m_screenBuffer;
    Point m_origin;  ///< Point in screen coordinates that all drawing functions are relative to.
    Rect m_clipRect; ///< In screen coordinates.
public:
    /// clipRect will be clipped to screenBuffer bounds.
    ScreenCanvas(ScreenBuffer& screenBuffer, Point origin, Rect clipRect);

    /// Sub canvas will be clipped to this canvas size.
    /// @param rect     Sub rectangle of this canvas. Rect is relative to this canvas' origin.
    ScreenCanvas getSubCanvas(Rect rect);

    /// Clears canvas to given color.
    void clear(Color bgColor) {
        auto localRect = m_clipRect;
        localRect.move(-m_origin.asSize()); // In local coordinates.
        fill(localRect, bgColor);
    }

    /// Draw filled rectangle.
    void fill(Rect rect, Color bgColor);

    /// Draw rectangle.
    void fillRect(Rect rect, bool doubleEdge, bool fill, Attributes attributes);

    /// Draws given text on the canvas.
    /// Text is clipped to boundaries of the canvas.
    /// So only graphemes fully inside are drawn.
    /// @param text     Input string, doesn't have to be valid or printable UTF-8.
    void print(Point pt, const std::string& text, Attributes normal, Attributes invalid, Attributes replacement);

    /// Draws given text on the canvas.
    /// Text is clipped to boundaries of the canvas.
    /// So only graphemes fully inside are drawn.
    /// @param graphemes    Graphemes to draw.
    void print(Point pt, gsl::span<const Grapheme> graphemes, Attributes normal, Attributes invalid, Attributes replacement);
};

class ScreenBuffer {
public:
    struct Character {
        std::string text; ///< UTF-8 text to draw. If empty, then this place will be drawn by preceeding character with width > 1.
        int width;        ///< Width of text once it will be rendered.
        Attributes attributes;

        bool operator==(const Character& other) const {
            if (text != other.text) return false;
            if (width != other.width) return false;
            if (attributes.fgColor != other.attributes.fgColor) return false;
            if (attributes.bgColor != other.attributes.bgColor) return false;
            if (attributes.style != other.attributes.style) return false;
            return true;
        }
    };

private:
    Size size;
    std::vector<Character> characters;
    std::vector<Character> previousCharacters;
    bool fullRepaintNeeded;     ///< If true while screen will be repainted. Otherwise only changes from previousCharacters will be repainted.

public:
    ScreenBuffer()
        : fullRepaintNeeded(true) {
    }

    ScreenCanvas getCanvas() {
        return ScreenCanvas(*this, Point{0, 0}, getSize());
    }

    /// Sets fullRepaintNeeded flag. Used when contents of the screen were changed without ScreenBuffer.
    void setFullRepaintNeeded() {
        fullRepaintNeeded = true;
    }

    Size getSize() const {
        return size;
    }

    int getWidth() const {
        return size.width;
    }

    int getHeight() const {
        return size.height;
    }

    /// Resizes this screen buffer.
    void resize(int w, int h);

    /// Clears screen to given color.
    void clear(Color bgColor);

    /// Draws a filled rectangle with given color.
    /// Rectangle is first clipped to fit the scren buffer.
    void fillRect(Rect rect, Color bgColor);

    /// Draws given text on the screen.
    /// Throws if text is not entirely on the screen.
    /// @param text     Input string, doesn't have to be valid or printable UTF-8.
    void print(int x, int y, const std::string& text, Attributes attributes);

    /// Draws given graphemes on the screen.
    /// Throws if text is not entirely on the screen.
    /// @param graphemes    Graphemes to draw.
    void print(int x, int y, gsl::span<const Grapheme> graphemes, Attributes attributes);

    /// Draws this screen buffer to the console.
    void present();
};

/// Returns Grapheme that corresponds to given simpleCharacter.
/// @param simpleCharacter      Must be a valid UTF-8 string, that converts to a printable character of width 1.
template<int N>
Grapheme simpleGrapheme(const char (&simpleChar)[N]) {
    return { GraphemeKind::NORMAL, simpleChar, "", 1, {simpleChar, N - 1} };
}

/// Draws a rectangle with borders.
/// Rectangle is clipped by the clipRect. clipRect must be wholy inside screen buffer.
void draw_rect(ScreenBuffer& screenBuffer, Rect clipRect, Rect rect, bool doubleEdge, bool fill, Attributes attributes);

/// Measures given text on the terminal
/// @param eventQueue   Event queue to use for listening for the result.
/// @param codePoints   List of code points to measure.
/// @return Length of the code points. Can be zero.
int measureText(EventQueue& eventQueue, gsl::span<uint32_t> codePoints);

} // namespace terminal_editor

ext: h

max_stars_repo_name: Zbyl/terminal-editor

max_issues_repo_name: Zbyl/terminal-editor

max_forks_repo_name: Zbyl/terminal-editor


text:
#ifndef CTETRA_INPUT_H
#define CTETRA_INPUT_H

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

typedef void (*InputFn)(double k[3], gsl_vector *values);

typedef void (*UEInputFn)(double k[3], gsl_vector *evalues, gsl_matrix_complex *evecs);

#endif // CTETRA_INPUT_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include "atomic.h"
#include "monotonic.h"

#include <memcached/types.h>
#include <gsl/gsl>

#include <functional>
#include <unordered_map>

struct DocKey;

namespace Collections {

// The reserved name of the system owned, default collection.
const char* const DefaultCollectionName = "_default";
static std::string_view DefaultCollectionIdentifier(DefaultCollectionName);

const char* const DefaultScopeName = "_default";
static std::string_view DefaultScopeIdentifier(DefaultScopeName);

// The SystemEvent keys are given some human readable tags to make disk or
// memory dumps etc... more helpful.
const char* const CollectionEventDebugTag = "_collection";
const char* const ScopeEventDebugTag = "_scope";

// Couchstore private file name for manifest data
const char CouchstoreManifest[] = "_local/collections_manifest";

// Name of file where the manifest will be kept on persistent buckets
const char* const ManifestFile = "collections.manifest";
static std::string_view ManifestFileName(ManifestFile);

// Length of the string excluding the zero terminator (i.e. strlen)
const size_t CouchstoreManifestLen = sizeof(CouchstoreManifest) - 1;

using ManifestUid = WeaklyMonotonic<uint64_t>;

// Struct/Map used in summary stat collecting (where we do vb accumulation)
struct AccumulatedStats {
    AccumulatedStats& operator+=(const AccumulatedStats& other);
    uint64_t itemCount{0};
    uint64_t diskSize{0};
    uint64_t opsStore{0};
    uint64_t opsDelete{0};
    uint64_t opsGet{0};
};
using Summary = std::unordered_map<CollectionID, AccumulatedStats>;

struct ManifestUidNetworkOrder {
    explicit ManifestUidNetworkOrder(ManifestUid uid) : uid(htonll(uid)) {
    }
    ManifestUid to_host() const {
        return ManifestUid(ntohll(uid));
    }
    ManifestUid::value_type uid;
};
static_assert(sizeof(ManifestUidNetworkOrder) == 8,
              "ManifestUidNetworkOrder must have fixed size of 8 bytes as "
              "written to disk.");

/**
 * Return a ManifestUid from a C-string.
 * A valid ManifestUid is a C-string where each character satisfies
 * std::isxdigit and can be converted to a ManifestUid by std::strtoull.
 *
 * @param uid C-string uid
 * @param len a length for validation purposes
 * @throws std::invalid_argument if uid is invalid
 */
ManifestUid makeUid(const char* uid, size_t len = 16);

/**
 * Return a ManifestUid from a std::string
 * A valid ManifestUid is a std::string where each character satisfies
 * std::isxdigit and can be converted to a ManifestUid by std::strtoull.
 *
 * @param uid std::string
 * @throws std::invalid_argument if uid is invalid
 */
static inline ManifestUid makeUid(const std::string& uid) {
    return makeUid(uid.c_str());
}

/**
 * Return a CollectionID from a C-string.
 * A valid CollectionID is a C-string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul.
 *
 * @param uid C-string uid
 * @throws std::invalid_argument if uid is invalid
 */
static inline CollectionID makeCollectionID(const char* uid) {
    // CollectionID is 8 characters max and smaller than a ManifestUid
    return gsl::narrow_cast<CollectionID>(makeUid(uid, 8));
}

/**
 * Return a CollectionID from a std::string
 * A valid CollectionID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 *
 * @param uid std::string
 * @throws std::invalid_argument if uid is invalid
 */
static inline CollectionID makeCollectionID(const std::string& uid) {
    return makeCollectionID(uid.c_str());
}

/**
 * Return a ScopeID from a C-string.
 * A valid CollectionID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 * @param uid C-string uid
 * @return std::invalid_argument if uid is invalid
 */
static inline ScopeID makeScopeID(const char* uid) {
    // ScopeId is 8 characters max and smaller than a ManifestUid
    return gsl::narrow_cast<ScopeID>(makeUid(uid, 8));
}

/**
 * Return a ScopeID from a std::string
 * A valid ScopeID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 * @param uid std::string
 * @return std::invalid_argument if uid is invalid
 */
static inline ScopeID makeScopeID(const std::string& uid) {
    return makeScopeID(uid.c_str());
}

/**
 * The metadata of a single collection
 *
 * Default construction yields the default collection
 */
struct CollectionMetaData {
    ScopeID sid{ScopeID::Default}; // The scope that the collection belongs to
    CollectionID cid{CollectionID::Default}; // The collection's ID
    std::string name{DefaultCollectionName}; // The collection's name
    cb::ExpiryLimit maxTtl{}; // The collection's maxTTL

    bool operator==(const CollectionMetaData& other) const {
        return sid == other.sid && cid == other.cid && name == other.name &&
               maxTtl == other.maxTtl;
    }
};

/**
 * The metadata of a single scope
 *
 * Default construction yields the default scope
 */
struct ScopeMetaData {
    ScopeID sid{ScopeID::Default}; // The scope's ID
    std::string name{DefaultScopeName}; // The scope's name

    bool operator==(const ScopeMetaData& other) const {
        return sid == other.sid && name == other.name;
    }
};

/**
 * For creation of collection SystemEvents - The SystemEventFactory
 * glues the CollectionID into the event key (so create of x doesn't
 * collide with create of y). This method yields the 'keyExtra' parameter
 *
 * @param collection The value to turn into a string
 * @return the keyExtra parameter to be passed to SystemEventFactory
 */
std::string makeCollectionIdIntoString(CollectionID collection);

/**
 * For creation of collection SystemEvents - The SystemEventFactory
 * glues the CollectionID into the event key (so create of x doesn't
 * collide with create of y). This method basically reverses
 * makeCollectionIdIntoString so we can get a CollectionID from a
 * SystemEvent key
 *
 * @param key DocKey from an Item in the System namespace and is a collection
 *        event
 * @return the ID which was in the event
 */
CollectionID getCollectionIDFromKey(const DocKey& key);

/// Same as getCollectionIDFromKey but for events changing scopes
ScopeID getScopeIDFromKey(const DocKey& key);

/**
 * Callback function for processing against dropped collections in an ephemeral
 * vb, returns true if the key at seqno should be dropped
 *
 * @param DocKey the key of the item we should process
 * @param int64_t the seqno of the item
 */
using IsDroppedEphemeralCb = std::function<bool(const DocKey&, int64_t)>;

/**
 * A function for determining if a collection is visible
 */
using IsVisibleFunction =
        std::function<bool(ScopeID, std::optional<CollectionID>)>;

namespace VB {
enum class ManifestUpdateStatus {
    Success,
    // The new Manifest has a 'UID' that is < current.
    Behind,
    // The new Manifest has a 'UID' that is == current, but adds/drops
    // collections/scopes.
    EqualUidWithDifferences,
    // The new Manifest changes scopes or collections immutable properties, e.g.
    // current has {id:8, name:"c1"} and new has {id:8,name:"c2"}.
    ImmutablePropertyModified
};
std::string to_string(ManifestUpdateStatus);

/// values required by the flusher to calculate new collection statistics
struct StatsForFlush {
    uint64_t itemCount;
    size_t diskSize;
    uint64_t highSeqno;
};

// Following classes define the metadata that will be held by the Manager but
// referenced from VB::Manifest
class CollectionSharedMetaData;
class CollectionSharedMetaDataView {
public:
    CollectionSharedMetaDataView(std::string_view name,
                                 ScopeID scope,
                                 cb::ExpiryLimit maxTtl);
    CollectionSharedMetaDataView(const CollectionSharedMetaData&);
    std::string to_string() const;
    std::string_view name;
    const ScopeID scope;
    const cb::ExpiryLimit maxTtl;
};

// The type stored by the Manager SharedMetaDataTable
class CollectionSharedMetaData : public RCValue {
public:
    CollectionSharedMetaData(std::string_view name,
                             ScopeID scope,
                             cb::ExpiryLimit maxTtl);
    CollectionSharedMetaData(const CollectionSharedMetaDataView& view);
    bool operator==(const CollectionSharedMetaDataView& view) const;
    bool operator!=(const CollectionSharedMetaDataView& view) const {
        return !(*this == view);
    }
    bool operator==(const CollectionSharedMetaData& meta) const;
    bool operator!=(const CollectionSharedMetaData& meta) const {
        return !(*this == meta);
    }

    const std::string name;
    const ScopeID scope;
    const cb::ExpiryLimit maxTtl;
};
std::ostream& operator<<(std::ostream& os,
                         const CollectionSharedMetaData& meta);

class ScopeSharedMetaData;
class ScopeSharedMetaDataView {
public:
    ScopeSharedMetaDataView(const ScopeSharedMetaData&);
    ScopeSharedMetaDataView(std::string_view name) : name(name) {
    }
    std::string to_string() const;
    std::string_view name;
};

// The type stored by the Manager SharedMetaDataTable
class ScopeSharedMetaData : public RCValue {
public:
    ScopeSharedMetaData(const ScopeSharedMetaDataView& view);
    bool operator==(const ScopeSharedMetaDataView& view) const;
    bool operator!=(const ScopeSharedMetaDataView& view) const {
        return !(*this == view);
    }
    bool operator==(const ScopeSharedMetaData& meta) const;
    bool operator!=(const ScopeSharedMetaData& meta) const {
        return !(*this == meta);
    }

    const std::string name;
};
std::ostream& operator<<(std::ostream& os, const ScopeSharedMetaData& meta);

} // namespace VB

} // end namespace Collections


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

int main(int argc, char * argv[])
{
  if (argc != 3 || strcmp(argv[1], "-h") == 0)
  {
    printf("Usage: ./gslDeterm N_MATRIX(<=1000) ITERATIONS(~1000)\n");
    return 1;
  }

  srand(time(NULL));
  int MAT_DIM = 6;
  int N_MATRIX = atoi(argv[1]);
  int ITERATIONS = atoi(argv[2]);
  int i,j,k, count;
  double start, end, t_time = 0.0, dummy;

  printf("N_MATRIX: %d, ITERATIONS: %d\n", N_MATRIX, ITERATIONS);

  gsl_permutation *p;
  gsl_matrix *m_list[N_MATRIX];
  gsl_vector *x = gsl_vector_alloc (MAT_DIM);
  int signum;

  double b_data[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
  gsl_vector_view b = gsl_vector_view_array (b_data, 6);

  // Do this ITERATIONS times, timing only the loop with the determinant
  // finding
  for (count=0; count<ITERATIONS; count++)
  {
    // Initialise a list of gsl arrays
    for (i=0; i<N_MATRIX; i++)
     m_list[i] = gsl_matrix_alloc(MAT_DIM,MAT_DIM);

    // Initialise values to random integers 0 -> 99
    for (i=0; i<N_MATRIX; i++)
      for (j=0; j<MAT_DIM; j++)
        for (k=0; k<MAT_DIM; k++)
          gsl_matrix_set (m_list[i], j, k, rand()%100);

    p = gsl_permutation_alloc(m_list[0]->size1);

    //start time: find the determinant of each matrix in
    //m_list
    start = (double)clock() /(double) CLOCKS_PER_SEC;

    for (i=0; i<N_MATRIX; i++)
    {
      gsl_linalg_LU_decomp(m_list[i], p, &signum);
      gsl_linalg_LU_solve(m_list[i], p, &b.vector, x);
    }

    end = (double)clock() / (double) CLOCKS_PER_SEC;
    t_time += end - start;
  }

  printf("GSL: %d matrices, %d iterations, cpu time: %fs\n",
                    N_MATRIX, ITERATIONS, t_time);

  gsl_permutation_free(p);
  for (i=0; i<N_MATRIX; i++)
    gsl_matrix_free(m_list[i]);
  gsl_vector_free (x);

  return 0;
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar


text:
#pragma once

#include <gsl.h>

#include "TMP_Helper.h"
#include "HE_Assert.h"
#include "HE_Math.h"

#include <type_traits>

namespace HE
{
	struct MemoryBlock
	{
		void* ptr;
		size_t length;

		void* begin() const { return ptr; }
		void* end() const { return static_cast<char*>(ptr) + length; }
	};
	using Blk = MemoryBlock;

	// Max alignment
	constexpr size_t PlatformMaxAlignment = Math::Max(alignof(void*), alignof(size_t));
	static_assert(Math::IsPow2(PlatformMaxAlignment), "PlatformMaxAlignment is not a power of 2, as should be");

	namespace Private
	{
		template<class T>
		using try_allocate = std::enable_if_t<std::is_same<Blk, decltype(std::declval<T>().allocate(std::declval<size_t>()))>::value>;

		// Custom aligned allocation's alignment parameter generally has to be a power of 2 and bigger than the allocator's default alignment,
		// in order to respect the allocator's semantics
		template<class T>
		using try_aligned_allocate = std::enable_if_t<std::is_same<Blk, decltype(std::declval<T>().allocate(std::declval<size_t>(), std::declval<size_t>()))>::value>;

		template<class T>
		using try_deallocate = decltype(std::declval<T>().deallocate(std::declval<Blk>()));

		template<class T>
		using try_deallocateAll = decltype(std::declval<T>().deallocateAll());

		template<class T>
		using try_owns = std::enable_if_t<std::is_same<bool, decltype(std::declval<T>().owns(std::declval<Blk>()))>::value>;

		template<class T>
		using try_it = std::enable_if_t<std::is_same<T, decltype(T::it)>::value>;
	}

	// Allocator
	// Must have (for allocator of type T):
	// - MemoryBlock T::allocate(size_t)
	// - void T::deallocate(MemoryBlock) noexcept
	template<class T>
	using is_allocator = and_ < has_op<T, Private::try_allocate>, has_op<T, Private::try_deallocate> > ;
	template<class... T>
	constexpr bool IsAllocator()
	{
		return and_<is_allocator<T>...>::value;
	}

	// StatelessAllocator
	// - Allocator
	// - StateSize<T>::value == 0
	// - is_same<decltype(T::it), T>
	template<class T>
	using is_stateless_allocator = and_<is_allocator<T>, equal_<StateSize<T>, std::integral_constant<size_t, 0>>>;
	template<class... T>
	constexpr bool IsStatelessAllocator()
	{
		return and_<is_stateless_allocator<T>...>::value;
	}

	// OwningAllocator
	// Must have (for allocator of type T)
	// - HE::is_allocator<T>()
	// - bool T::owns(MemoryBlock)
	template<class T>
	using is_owning_allocator = and_<is_allocator<T>, has_op<T, Private::try_owns>>;
	template<class... T>
	constexpr bool IsOwningAllocator()
	{
		return and_<is_owning_allocator<T>...>::value;
	}

	// AlignedAllocator
	// Must have (for allocator of type T)
	// - HE::is_allocator<T>()
	// - bool T::allocate(size_t, size_t)
	template<class T>
	using is_aligned_allocator = and_<is_allocator<T>, has_op<T, Private::try_aligned_allocate>>;
	template<class... T>
	constexpr bool IsAlignedAllocator()
	{
		return and_<is_aligned_allocator<T>...>::value;
	}

	// Generic allocator functions
	template< class Type, class Allocator, class Enable = std::enable_if_t<is_allocator<Allocator>::value>>
	Blk allocate(Allocator&& a)
	{
		return a.allocate(sizeof(Type));
	}

	template<class Type, class Allocator, class Enable = std::enable_if_t<is_allocator<Allocator>::value>>
	Blk allocate(Allocator&& a, size_t count)
	{
		return a.allocate(count*sizeof(Type), alignof(Type));
	}

	class NullAllocator
	{
	public:
		static constexpr size_t alignment = 64 * 1024;
		static NullAllocator it;

		Blk allocate(size_t);
		Blk allocate(size_t, size_t);
		void deallocate(Blk) noexcept;
		void deallocateAll() noexcept;
		bool owns(Blk);
	};

	// Returns the beginning of the buffer is the buffer is big enough, even if it
	// was already allocated. This is because the allocator does no tracking. The client
	// is responsible for making sure memory doesn't get corrupted
	// Probably shouldn't be used as the main allocator of the Fallback allocator, but pairs
	// well with Segregator
	template<size_t N>
	class alignas(PlatformMaxAlignment) LightInlineAllocator
	{
	public:
		static constexpr size_t alignment = PlatformMaxAlignment;

		Blk allocate(size_t n)
		{
			if (n <= N)
			{
				return{ m_buffer, n };
			}
			else
			{
				return{ nullptr, 0 };
			}
		}

		Blk allocate(size_t n, size_t a)
		{
			EXPECTS(Math::IsPow2(alignment) && a >= alignment);
			auto const p = reinterpret_cast<char*>(Math::RoundUpToMultipleOf(reinterpret_cast<size_t>(&m_buffer[0]), a));
			Blk const b{ p, n };
			if (b.end() <= m_buffer + N)
			{
				return b;
			}
			else
			{
				return{ nullptr, 0 };
			}
		}

		bool owns(Blk b)
		{
			return b.begin() >= m_buffer && b.end() <= m_buffer + N;
		}

		void deallocate(Blk) noexcept {}

	private:
		char m_buffer[N];
	};

	class MallocAllocator
	{
	public:
		static constexpr size_t alignment = PlatformMaxAlignment;
		static MallocAllocator it;

		Blk allocate(size_t n);
		void deallocate(Blk) noexcept;
	};

	class AlignedMallocAllocator
	{
	public:
		static constexpr size_t alignment = PlatformMaxAlignment;
		static AlignedMallocAllocator it;

		Blk allocate(size_t n);
		Blk allocate(size_t n, size_t alignment);
		void deallocate(Blk) noexcept;
	};

	template< class Primary, class Fallback >
	class FallbackAllocator;

	namespace Private
	{
		template < class Primary, class Fallback, class Enable = void>
		struct FallbackAllocatorImpl {	};

		template< class Primary, class Fallback >
		using FallbackStatelessCond = std::enable_if_t<
			and_<
			equal_<StateSize<Primary>, std::integral_constant<size_t, 0>>,
			equal_<StateSize<Fallback>, std::integral_constant<size_t, 0>>
			>::value
		>;

		template < class Primary, class Fallback >
		struct FallbackAllocatorImpl<Primary, Fallback,
			FallbackStatelessCond<Primary, Fallback> >
		{
			static FallbackAllocator<Primary, Fallback> it;
		};
	}

	template< class Primary, class Fallback >
	class FallbackAllocator :
		private Primary,
		private Fallback,
		private Private::FallbackAllocatorImpl<Primary, Fallback>
	{
	protected:
		using P = Primary;
		using F = Fallback;

		static_assert(IsOwningAllocator<P>(), "Primary allocator does not meet the HE::OwningAllocator concept");
		static_assert(IsAllocator<F>(), "Fallback allocator does not meet the HE::Allocator concept");

	public:
		static constexpr size_t alignment = Math::Min(Primary::alignment, Fallback::alignment);

		Blk allocate(size_t n)
		{
			auto const blk = P::allocate(n);
			if (!blk.ptr) return F::allocate(n);
			return blk;
		}

		template<class Enable = std::enable_if_t<and_<is_aligned_allocator<Primary>, is_aligned_allocator<Fallback>>::value>>
		Blk allocate(size_t n, size_t alignment)
		{
			auto const blk = P::allocate(n, alignment);
			if (!blk.ptr) return F::allocate(n, alignment);
			return blk;
		}

		void deallocate(Blk b) noexcept
		{
			if (P::owns(b))
				P::deallocate(b);
			else
				F::deallocate(b);
		}

		template<class Enable = std::enable_if_t<and_<has_op<Primary, Private::try_deallocateAll>, has_op<Fallback, Private::try_deallocateAll>>::value>>
		void deallocateAll() noexcept
		{
			Primary::deallocateAll();
			Fallback::deallocateAll();
		}

		template<class Enable = std::enable_if_t<is_owning_allocator<Fallback>::value>>
		bool owns(Blk b)
		{
			return Primary::owns(b) || Fallback::owns(b);
		}
	};

	template< class Primary, class Fallback >
	FallbackAllocator<Primary, Fallback> Private::FallbackAllocatorImpl<Primary, Fallback, Private::FallbackStatelessCond<Primary, Fallback>>::it;

	namespace Allocator
	{
		constexpr size_t unbounded = static_cast<size_t>(-1);
	}

	// An allocator that allocates nodes on a Parent allocator and internally keeps blocks
	// of memory on deallocation instead of actually deallocating them, but only
	// if they have a certain size
	// Important: Actual deallocations to the parent allocator are not guaranteed to be ordered as
	// deallocated in the Freelist allocator, and therefore should not be used with allocators
	// that require ordered deallocations, such as StackAllocator

	// minSize: Minimum size of the allocation to be considered in range
	// maxSize: Maxsimum size of the allocation to be considered in range
	// batchCount: Number of allocations on a fresh "in range" allocation. Basically fills up the Freelist with batchCount nodes on allocation
	// maxNodes: Max number of nodes the freelist can keep
	template< class Parent,
		size_t MinSize,
		size_t MaxSize = MinSize,
		size_t MaxNodes = Allocator::unbounded,
		class Enable = std::enable_if_t<is_allocator<Parent>::value>
		>
	class FreelistAllocator
		: private Parent
	{
		static_assert(MaxSize >= MinSize, "FreelistAllocator's MaxSize should be higher or equal to MinSize");
		static_assert(MaxSize >= sizeof(void*), "FreelistAllocator's MaxSize and MinSize should be higher or equal than sizeof(void*)");

	public:
		static constexpr size_t alignment = Parent::alignment;

		Blk allocate(size_t n)
		{
			return allocateImpl(n);
		}

		template<class Enable = std::enable_if_t<is_aligned_allocator<Parent>::value>>
		Blk allocate(size_t n, size_t alignment)
		{
			return allocateImpl(n, alignment);
		}

		void deallocate(Blk b) noexcept
		{
			if ((MaxNodes == Allocator::unbounded || m_nNodesCount != MaxNodes) && inRange(b.length))
			{
				auto const next = m_pFreelistRoot;
				m_pFreelistRoot = static_cast<Node*>(b.ptr);
				m_pFreelistRoot->next = next;
				++m_nNodesCount;
			}
			else
			{
				Parent::deallocate(b);
			}
		}

		// Only O(1) if the Parent allocator supports deallocateAll
		// If the Freelist is unbounded and the Parent doesn't support deallocateAll,
		// the Freelist will do a best effort of deallocating all the nodes in its freelist
		// in O(n)
		void deallocateAll() noexcept
		{
			deallocateAllImpl();
		}

		static constexpr bool has_fast_deallocateAll() { return has_op<Parent, Private::try_deallocateAll>::value; }

		bool owns(Blk b)
		{
			return Parent::owns(b);
		}

	private:
		struct Node
		{
			Node* next;
		};
		Node* m_pFreelistRoot{ nullptr };
		size_t m_nNodesCount{ 0 };

		bool inRange(size_t n) const
		{
			if (MinSize == MaxSize) return n == MaxSize;

			return (MinSize == 0 || n >= MinSize) && n <= MaxSize;
		}

		template<class... Args>
		Blk allocateImpl(size_t n, Args... args)
		{
			if (!inRange(n)) return Parent::allocate(n, args...);

			if (!m_pFreelistRoot)
			{
				auto const b = Parent::allocate(MaxSize, args...);
				return{ b.ptr, n };
			}
			else
			{
				Blk const result{ static_cast<void*>(m_pFreelistRoot), n };
				m_pFreelistRoot = m_pFreelistRoot->next;
				--m_nNodesCount;
				return result;
			}
		}

		void deallocateAllImpl(std::enable_if_t<has_op<Parent, Private::try_deallocateAll>::value>* = nullptr) noexcept
		{
			Parent::deallocateAll();
			m_pFreelistRoot = nullptr;
		}

		// In this case, we can only guarantee a complete deallocation if the freelist is unbounded
		void deallocateAllImpl(std::enable_if_t<
			and_<
			not_<has_op<Parent, Private::try_deallocateAll>>,
			has_op<Parent, Private::try_deallocate>,
			equal_<std::integral_constant<size_t, MaxNodes>, std::integral_constant<size_t, Allocator::unbounded>>
			>::value>* = nullptr) noexcept
		{
			auto next = m_pFreelistRoot;
			while (next)
			{
				Blk const b{ static_cast<void*>(next), MaxSize };
				next = next->next;
				Parent::deallocate(b);
			}
			m_pFreelistRoot = nullptr;
		}
	};

	template<class Parent, class PrefixType, class SuffixType>
	class AffixAllocator;

	namespace Private
	{
		template < class Parent, class PrefixType, class SuffixType, class Enable = void>
		struct AffixParentImpl : protected Parent {	};

		template < class Parent, class PrefixType, class SuffixType>
		struct AffixParentImpl<Parent, PrefixType, SuffixType, std::enable_if_t<StateSize<Parent>::value == 0>>
			: protected Parent
		{
			static AffixAllocator<Parent, PrefixType, SuffixType> it;
		};

		template<class PrefixType>
		struct PrefixImpl
		{
			template<class Enable = std::enable_if_t<StateSize<PrefixType>::value != 0>>
			static PrefixType& Prefix(Blk& b)
			{
				return reinterpret_cast<PrefixType*>(b.ptr)[-1];
			}
		};

		template<>
		struct PrefixImpl<void>{};

		template<class PrefixType, class SuffixType>
		struct SuffixImpl
		{
			template<class Enable = std::enable_if_t<StateSize<SuffixType>::value != 0>>
			static SuffixType& Suffix(Blk& b)
			{

				auto const p = static_cast<char*>(b.ptr) + b.length;
				return *reinterpret_cast<SuffixType*>(p);
			}

		protected:
			static size_t totalAllocationSize(size_t s)
			{
				if (StateSize<SuffixType>::value == 0)
				{
					return s + StateSize<PrefixType>::value;
				}
				else
				{
					return Math::RoundUpToMultipleOf(s + StateSize<PrefixType>::value, alignof(SuffixType)) + StateSize<SuffixType>::value;
				}
			}
		};

		template<class PrefixType>
		struct SuffixImpl<PrefixType, void>
		{
		protected:
			static size_t totalAllocationSize(size_t s)
			{
				return s + StateSize<PrefixType>::value;
			}
		};
	}

	template<class Parent, class PrefixType, class SuffixType = void>
	class AffixAllocator : public Private::AffixParentImpl<Parent, PrefixType, SuffixType>,
		public Private::PrefixImpl<PrefixType>,
		public Private::SuffixImpl<PrefixType, SuffixType>
	{
	public:
		static constexpr size_t alignment = StateSize<PrefixType>::value ? alignof(PrefixType) : Parent::alignment;

		Blk allocate(size_t n)
		{
			auto const result = Parent::allocate(totalAllocationSize(n));
			if (!result.ptr) return result;

			return{ static_cast<char*>(result.ptr) + StateSize<PrefixType>::value, n };
		}

		bool owns(Blk b)
		{
			return Parent::owns(actualAllocation(b));
		}

		void deallocate(Blk b)
		{
			Parent::deallocate(actualAllocation(b));
		}

	private:
		// Takes a requested allocation, and returns the actual allocation that was request to the parent
		static Blk actualAllocation(Blk b)
		{
			if (!b.ptr) return{ nullptr, 0 };
			return{ static_cast<char*>(b.ptr) - StateSize<PrefixType>::value, totalAllocationSize(b.length) };
		}
	};

	template<class Parent, class PrefixType, class SuffixType>
	AffixAllocator<Parent, PrefixType, SuffixType> Private::AffixParentImpl<Parent, PrefixType, SuffixType, std::enable_if_t<StateSize<Parent>::value == 0>>::it;

	template < size_t Threshold, class SmallAllocator, class LargeAllocator >
	class SegregateAllocator;

	namespace Private
	{
		template < size_t Threshold, class SmallAllocator, class LargeAllocator, class Enable = void>
		struct SegregateAllocatorImpl {	};

		template< class SmallAllocator, class LargeAllocator >
		using SegregateAllocatorStatelessCond = std::enable_if_t<
			and_<
			equal_<StateSize<SmallAllocator>, std::integral_constant<size_t, 0>>,
			equal_<StateSize<LargeAllocator>, std::integral_constant<size_t, 0>>
			>::value
		>;

		template < size_t Threshold, class SmallAllocator, class LargeAllocator >
		struct SegregateAllocatorImpl<Threshold, SmallAllocator, LargeAllocator,
			SegregateAllocatorStatelessCond<SmallAllocator, LargeAllocator> >
		{
			static SegregateAllocator<Threshold, SmallAllocator, LargeAllocator> it;
		};
	}

	// Seggregator
	template<size_t Threshold,
		class SmallAllocator,
		class LargeAllocator >
	class SegregateAllocator
		: private SmallAllocator,
		private LargeAllocator
	{
	public:
		constexpr static size_t alignment = Math::Min(SmallAllocator::alignment, LargeAllocator::alignment);

		Blk allocate(size_t n)
		{
			return allocate_Impl(n);
		}

		Blk allocate(size_t n, size_t alignment_requirement)
		{
			EXPECTS(alignment_requirement >= alignment && Math::IsPow2(alignment_requirement));
			return allocate_Impl(n, alignment_requirement);
		}

		bool owns(Blk b)
		{
			return n <= Threshold ? SmallAllocator::owns(b) : LargeAllocator::owns(b);
		}

		void deallocate(Blk b)
		{
			return n <= Threshold ? SmallAllocator::deallocate(b) : LargeAllocator::deallocate(b);
		}

		void deallocateAll()
		{
			SmallAllocator::deallocateAll();
			LargeAllocator::deallocateAll();
		}

	private:
		template<class... Args>
		Blk allocate_Impl(size_t n, Args... args)
		{
			return n <= Threshold ? SmallAllocator::allocate(n, args...) : LargeAllocator::allocate(n, args...);
		}
	};

	template < size_t Threshold, class SmallAllocator, class LargeAllocator >
	SegregateAllocator<Threshold, SmallAllocator, LargeAllocator> Private::SegregateAllocatorImpl<Threshold, SmallAllocator, LargeAllocator, Private::SegregateAllocatorStatelessCond<SmallAllocator, LargeAllocator>>::it;
}


text:
//
// Created by Huy Vo on 2019-06-29.
//

#ifndef PACMENSL_SRC_PETSCWRAP_PETSCWRAP_H_
#define PACMENSL_SRC_PETSCWRAP_PETSCWRAP_H_
#include <petsc.h>
#include <iostream>
#include "Sys.h"
namespace pacmensl {

template<typename PetscT>
class Petsc {
 protected:
  PetscT dat = nullptr;

  int Destroy(PetscT *obj);

 public:

  Petsc() {
    static_assert(std::is_convertible<PetscT, Vec>::value || std::is_convertible<PetscT, Mat>::value
                      || std::is_convertible<PetscT, IS>::value || std::is_convertible<PetscT, VecScatter>::value
                      || std::is_convertible<PetscT, KSP>::value || std::is_convertible<PetscT, TS>::value,
                  "pacmensl::Petsc can only wrap PETSc objects.");
  }

  PetscT *mem() { return &dat; }

  const PetscT *mem() const { return &dat; }

  bool IsEmpty() { return (dat == nullptr); }

  operator PetscT() { return dat; }

  ~Petsc() {
    Destroy(&dat);
  }
};

template<>
int Petsc<Vec>::Destroy(Vec *obj);

template<>
int Petsc<Mat>::Destroy(Mat *obj);

template<>
int Petsc<IS>::Destroy(IS *obj);

template<>
int Petsc<VecScatter>::Destroy(VecScatter *obj);

template<>
int Petsc<KSP>::Destroy(KSP *obj);

template<>
int Petsc<TS>::Destroy(TS *obj);

PacmenslErrorCode ExpandVec(Petsc<Vec> &p, const std::vector<PetscInt> &new_indices, const PetscInt new_local_size);
PacmenslErrorCode ExpandVec(Vec &p, const std::vector<PetscInt> &new_indices, const PetscInt new_local_size);
}
#endif //PACMENSL_SRC_PETSCWRAP_PETSCWRAP_H_

max_stars_repo_name: voduchuy/pacmensl

max_issues_repo_name: voduchuy/pacmensl

max_forks_repo_name: voduchuy/pacmensl


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_results_Filter_h_
#define SQ_INCLUDE_GUARD_results_Filter_h_

#include "core/Field.h"
#include "core/typeutil.h"
#include "parser/FilterSpec.h"

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

namespace sq::results {

struct Filter;
using FilterPtr = std::unique_ptr<Filter>;

struct Filter {
  /**
   * Create a Filter for the given spec.
   */
  SQ_ND static FilterPtr create(const parser::FilterSpec &spec);

  /**
   * Apply this filter to a Result.
   */
  SQ_ND virtual Result operator()(Result &&result) const = 0;

  virtual ~Filter() = default;
  Filter() = default;
  Filter(const Filter &) = delete;
  Filter(Filter &&) = delete;
  Filter &operator=(const Filter &) = delete;
  Filter &operator=(Filter &&) = delete;
};

} // namespace sq::results

#endif // SQ_INCLUDE_GUARD_results_Filter_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
/** \file oposim_fns.c
    \brief Funciones auxiliares
*/

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <time.h>
#include <stdlib.h>

#include "oposim.h"

static gsl_rng *r = NULL; /**< \a r es el generador de números aleatorios */

/** \fn void random_init(char opc)
\brief Inicializador del generador de números aleatorios.
\param opc a character.
 */
void random_init(char opc)
{

  if (!r) {
    switch(opc) {
    case 'm':
      r = gsl_rng_alloc(gsl_rng_mt19937);
      break;
    case 't':
      r = gsl_rng_alloc(gsl_rng_taus);
      break;
    case 'g':
      r = gsl_rng_alloc(gsl_rng_gfsr4);
      break;
    case 'x':
      r = gsl_rng_alloc(gsl_rng_ranlxs0);
      break;
    default:
      r = gsl_rng_alloc(gsl_rng_mrg);
    }
    printf("Usando RNG %s\n", gsl_rng_name(r));
    gsl_rng_set(r, time(NULL));
  }

}

/** \fn void random_free()
\brief Destructor del generador de números aleatorios
*/
void random_free()
{

  if (!r) {
    gsl_rng_free(r);
  }

}

/** \fn int s_experimento(int temas, int bolas, int estudiados)
\brief Realiza una simulación.
\param temas Número de temas en la oposición.
\param bolas Número de bolas que se extraen.
\param estudiados Número de temas estudiados.
\return Número temas estudiados que estarían entre los extraídos.
*/
int s_experimento(int temas, int bolas, int estudiados)
{

  //printf("bolas extraidas: ");
  int b_extraidas[bolas];
  int bola;

  for (int i = 0; i < bolas; i++) {
    do {
      bola = (int) gsl_rng_uniform_int(r, temas) + 1;
      for (int j = 0; j < i; j++) {
	if (b_extraidas[j] == bola) {
	  bola = 0;
	  break;
	}
      }
    } while (!bola);
    //printf("%d ", bola);
    b_extraidas[i] = bola;
  }

  int res = 0;
  for (int i = 0; i < bolas; i++) {
    if(b_extraidas[i] <= estudiados) res++;
  }

  //printf("\n");
  return res;

}

/** \fn int *experimento(int temas, int bolas, int estudiados, int reps)
\brief Realiza varias simulaciones y devuelve los resultados.
\param temas Número de temas en la oposición.
\param bolas Número de bolas que se extraen.
\param estudiados Número de temas estudiados.
\param reps Número de veces que se realizará el experimento.
\return Histograma con las veces que el número de temas estudiados estaba entre las bolas extraídas.
 */
int *experimento(int temas, int bolas, int estudiados, int reps)
{

  int *res = (int *)calloc(bolas + 1, sizeof(int));

  for (int i = 0; i < reps; i++)
    res[s_experimento(temas, bolas, estudiados)]++;

  return res;

}

/** \fn double *frecs(int *histograma, int l)
\brief Calcula las frecuencias relativas.
\param histograma array de frecuencias absolutas
\param l longitud del histograma
\return array con las frecuencias relativas
 */
double *frecs(int *histograma, int l)
{

  double suma = 0;
  for (int i = 0; i < l; i++)
    suma += histograma[i];

  double *res = (double *)calloc(l, sizeof(double));

  for (int i = 0; i < l; i++)
    res[i] = histograma[i] / suma;

  return res;

}

max_stars_repo_name: jonatanlv/oposim

max_issues_repo_name: jonatanlv/oposim

max_forks_repo_name: jonatanlv/oposim


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

//Spline creator
//n     -> number of points
//x     -> x-axis
//y     -> f(x)-axis
//y0,yf -> values of f(x) to use beyond the interpolation range
ccl_f1d_t *ccl_f1d_t_new(int n,double *x,double *y,double y0,double yf,
			 ccl_f1d_extrap_t extrap_lo_type,
			 ccl_f1d_extrap_t extrap_hi_type, int *status)
{
  ccl_f1d_t *spl=malloc(sizeof(ccl_f1d_t));
  if(spl==NULL) {
    *status = CCL_ERROR_MEMORY;
    return NULL;
  }

  spl->spline=gsl_spline_alloc(gsl_interp_akima,n);
  if (spl->spline == NULL) {
    *status = CCL_ERROR_MEMORY;
    free(spl);
    return NULL;
  }
  int parstatus=gsl_spline_init(spl->spline,x,y,n);
  if(parstatus) {
    *status = CCL_ERROR_SPLINE;
    gsl_spline_free(spl->spline);
    free(spl);
    return NULL;
  }

  if(*status==0) {
    spl->y0=y0;
    spl->yf=yf;
    spl->x_ini=x[0];
    spl->x_end=x[n-1];
    spl->y_ini=y[0];
    spl->y_end=y[n-1];
    spl->extrap_lo_type=extrap_lo_type;
    spl->extrap_hi_type=extrap_hi_type;
  }

  // Compute derivatives
  // Low-end
  if(spl->extrap_lo_type == ccl_f1d_extrap_const) {
    spl->der_lo = 0;
  }
  else if(spl->extrap_lo_type == ccl_f1d_extrap_linx_liny) {
    spl->der_lo = (y[1]-y[0])/(x[1]-x[0]);
  }
  else if(spl->extrap_lo_type == ccl_f1d_extrap_linx_logy) {
    if(y[1]*y[0]<=0)
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_lo = log(y[1]/y[0])/(x[1]-x[0]);
  }
  else if(spl->extrap_lo_type == ccl_f1d_extrap_logx_liny) {
    if(x[1]*x[0]<=0)
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_lo = (y[1]-y[0])/log(x[1]/x[0]);
  }
  else if(spl->extrap_lo_type == ccl_f1d_extrap_logx_logy) {
    if((y[1]*y[0]<=0) || (x[1]*x[0]<=0))
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_lo = log(y[1]/y[0])/log(x[1]/x[0]);
  }
  else // No extrapolation
    spl->der_lo = 0;

  // High-end
  if(spl->extrap_hi_type == ccl_f1d_extrap_const) {
    spl->der_hi = 0;
  }
  else if(spl->extrap_hi_type == ccl_f1d_extrap_linx_liny) {
    spl->der_hi = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
  }
  else if(spl->extrap_hi_type == ccl_f1d_extrap_linx_logy) {
    if(y[n-1]*y[n-2]<=0)
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_hi = log(y[n-1]/y[n-2])/(x[n-1]-x[n-2]);
  }
  else if(spl->extrap_hi_type == ccl_f1d_extrap_logx_liny) {
    if(x[n-1]*x[n-2]<=0)
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_hi = (y[n-1]-y[n-2])/log(x[n-1]/x[n-2]);
  }
  else if(spl->extrap_hi_type == ccl_f1d_extrap_logx_logy) {
    if((y[n-1]*y[n-2]<=0) || (x[n-1]*x[n-2]<=0))
      *status = CCL_ERROR_SPLINE;
    else
      spl->der_hi = log(y[n-1]/y[n-2])/log(x[n-1]/x[n-2]);
  }
  else // No extrapolation
    spl->der_hi = 0;

  // If one of the derivatives could not be calculated
  // then return NULL.
  if (*status){
    ccl_f1d_t_free(spl);
    spl = NULL;
  }

  return spl;
}

//Evaluates spline at x checking for bound errors
double ccl_f1d_t_eval(ccl_f1d_t *spl,double x)
{
  // Extrapolation
  // Lin_x, lin_y
  // f(x) = f_n + (x-x_n) * der
  //   der = (f_n - f_{n-1}) / (x_n - x_{n-1})
  //

  // Log_x, lin_y
  // f(x) = f_n + log(x/x_n) * der
  // der = (f_n - f_{n-1}) / log(x_n/x_{n-1})

  // Lin_x, log_y
  // f(x) = f_n * exp[ (x - x_n) * der ]
  // der = log(f_n/f_{n-1})/(x_n-x_{n-1})

  // Log_x, log_y
  // f(x) = f_n * exp[ log(x/x_n) * der ]
  // der = log(f_n/f_{n-1})/log(x_n/x_{n-1})
  if(x<=spl->x_ini) {
    if(spl->extrap_lo_type==ccl_f1d_extrap_const)
      return spl->y0;
    else if(spl->extrap_lo_type==ccl_f1d_extrap_linx_liny) {
      return spl->y_ini + spl->der_lo * (x - spl->x_ini);
    }
    else if(spl->extrap_lo_type==ccl_f1d_extrap_logx_liny) {
      return spl->y_ini + spl->der_lo * log(x/spl->x_ini);
    }
    else if(spl->extrap_lo_type==ccl_f1d_extrap_linx_logy) {
      return spl->y_ini * exp(spl->der_lo * (x - spl->x_ini));
    }
    else if(spl->extrap_lo_type==ccl_f1d_extrap_logx_logy) {
      return spl->y_ini * pow(x/spl->x_ini, spl->der_lo);
    }
    else {
      ccl_raise_gsl_warning(CCL_ERROR_SPLINE_EV,
                            "ccl_f1d.c: ccl_f1d_t_eval(): "
                            "x-value below range.");
      return NAN;
    }
  }
  else if(x>=spl->x_end) {
    if(spl->extrap_hi_type==ccl_f1d_extrap_const)
      return spl->yf;
    else if(spl->extrap_hi_type==ccl_f1d_extrap_linx_liny) {
      return spl->y_end + spl->der_hi * (x - spl->x_end);
    }
    else if(spl->extrap_hi_type==ccl_f1d_extrap_logx_liny) {
      return spl->y_end + spl->der_hi * log(x/spl->x_end);
    }
    else if(spl->extrap_hi_type==ccl_f1d_extrap_linx_logy) {
      return spl->y_end * exp(spl->der_hi * (x - spl->x_end));
    }
    else if(spl->extrap_hi_type==ccl_f1d_extrap_logx_logy) {
      return spl->y_end * pow(x/spl->x_end, spl->der_hi);
    }
    else {
      ccl_raise_gsl_warning(CCL_ERROR_SPLINE_EV,
                            "ccl_f1d.c: ccl_f1d_t_eval(): "
                            "x-value above range.");
      return NAN;
    }
  }
  else {
    double y;
    int stat=gsl_spline_eval_e(spl->spline,x,NULL,&y);
    if (stat!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(stat, "ccl_f1d.c: ccl_f1d_t_eval(): "
                            "x-value outside range.");
      return NAN;
    }
    return y;
  }
}

//Spline destructor
void ccl_f1d_t_free(ccl_f1d_t *spl)
{
  if (spl != NULL) {
    gsl_spline_free(spl->spline);
  }
  free(spl);
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
/* gsl_histogram2d_calloc_range.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library 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 library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
/***************************************************************
 *
 * File gsl_histogram2d_calloc_range.c:
 * Routine to create a variable binning 2D histogram providing
 * the input range vectors. Need GSL library and header.
 * Do range check and allocate the histogram data.
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram2d.h>
/*
 * Routine that create a 2D histogram using the given
 * values for X and Y ranges
 */
gsl_histogram2d *
gsl_histogram2d_calloc_range (size_t nx, size_t ny,
                              double *xrange,
                              double *yrange)
{
  size_t i, j;
  gsl_histogram2d *h;

  /* check arguments */

  if (nx == 0)
    {
      GSL_ERROR_VAL ("histogram length nx must be positive integer",
                        GSL_EDOM, 0);
    }

  if (ny == 0)
    {
      GSL_ERROR_VAL ("histogram length ny must be positive integer",
                        GSL_EDOM, 0);
    }

  /* init ranges */

  for (i = 0; i < nx; i++)
    {
      if (xrange[i] >= xrange[i + 1])
        {
          GSL_ERROR_VAL ("histogram xrange not in increasing order",
                            GSL_EDOM, 0);
        }
    }

  for (j = 0; j < ny; j++)
    {
      if (yrange[j] >= yrange[j + 1])
        {
          GSL_ERROR_VAL ("histogram yrange not in increasing order"
                            ,GSL_EDOM, 0);
        }
    }

  /* Allocate histogram  */

  h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d));

  if (h == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram struct",
                        GSL_ENOMEM, 0);
    }

  h->xrange = (double *) malloc ((nx + 1) * sizeof (double));

  if (h->xrange == 0)
    {
      /* exception in constructor, avoid memory leak */
      free (h);

      GSL_ERROR_VAL ("failed to allocate space for histogram xrange",
                        GSL_ENOMEM, 0);
    }

  h->yrange = (double *) malloc ((ny + 1) * sizeof (double));

  if (h->yrange == 0)
    {
      /* exception in constructor, avoid memory leak */
      free (h);

      GSL_ERROR_VAL ("failed to allocate space for histogram yrange",
                        GSL_ENOMEM, 0);
    }

  h->bin = (double *) malloc (nx * ny * sizeof (double));

  if (h->bin == 0)
    {
      /* exception in constructor, avoid memory leak */
      free (h->xrange);
      free (h->yrange);
      free (h);

      GSL_ERROR_VAL ("failed to allocate space for histogram bins",
                        GSL_ENOMEM, 0);
    }

  /* init histogram */

  /* init ranges */

  for (i = 0; i <= nx; i++)
    {
      h->xrange[i] = xrange[i];
    }

  for (j = 0; j <= ny; j++)
    {
      h->yrange[j] = yrange[j];
    }

  /* clear contents */

  for (i = 0; i < nx; i++)
    {
      for (j = 0; j < ny; j++)
        {
          h->bin[i * ny + j] = 0;
        }
    }

  h->nx = nx;
  h->ny = ny;

  return h;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

int
func (double t, const double y[], double f[],
      void *params)
{
  double mu = *(double *)params;
  f[0] = y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

int
main (void)
{
  double mu = 10;
  gsl_odeiv2_system sys = {func, jac, 2, &mu};

  gsl_odeiv2_driver * d =
    gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
				  1e-6, 1e-6, 0.0);
  int i;
  double t = 0.0, t1 = 100.0;
  double y[2] = { 1.0, 0.0 };

  for (i = 1; i <= 100; i++)
    {
      double ti = i * t1 / 100.0;
      int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

      if (status != GSL_SUCCESS)
	{
	  printf ("error, return value=%d\n", status);
	  break;
	}

      printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

  gsl_odeiv2_driver_free (d);
  return 0;
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "ellipsoid/ellipsoid.h"
#include "tanhsinh.h"

#undef __FUNCT__
#define __FUNCT__ "NormConstantIntFixed"
PetscErrorCode NormConstantIntFixed(EllipsoidalSystem *e, PetscInt n, PetscInt p, PetscInt prec, PetscInt nPts, PetscReal *intVals, PetscReal *normConst)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  FuncInfo2 ctx1 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = 1 };
  FuncInfo2 ctx2 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = 1 };
  FuncInfo2 ctx3 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = -1 };
  FuncInfo2 ctx4 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = -1 };

  PetscReal integrals[4];

  mpfr_t mpfrzero;
  mpfr_t mpfrone;
  mpfr_inits(mpfrzero, mpfrone, NULL);
  mpfr_set_d(mpfrzero, 0.0, MPFR_RNDN);
  mpfr_set_d(mpfrone, 1.0, MPFR_RNDN);

  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+0, &ctx1);
  intVals[0] = integrals[0];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+1, &ctx2);
  intVals[1] = integrals[1];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+2, &ctx3);
  intVals[2] = integrals[2];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+3, &ctx4);
  intVals[3] = integrals[3];
  *normConst = 8.0*(integrals[2]*integrals[1] - integrals[0]*integrals[3]); flopCount += 4;

  mpfr_clears(mpfrzero, mpfrone, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "NormConstantIntFixedSE"
PetscErrorCode NormConstantIntFixedSE(EllipsoidalSystem *e, PetscInt n, PetscInt p, PetscInt prec, PetscInt nPts, PetscReal *normConst)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  FuncInfo2 ctx1 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = 1 };
  FuncInfo2 ctx2 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = 1 };
  FuncInfo2 ctx3 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = -1 };
  FuncInfo2 ctx4 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = -1 };

  PetscReal integrals[4];

  mpfr_t mpfrzero;
  mpfr_t mpfrone;
  mpfr_inits(mpfrzero, mpfrone, NULL);
  mpfr_set_d(mpfrzero, 0.0, MPFR_RNDN);
  mpfr_set_d(mpfrone, 1.0, MPFR_RNDN);

  ierr = SEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+0, &ctx1);

  ierr = SEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+1, &ctx2);

  ierr = SEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+2, &ctx3);

  ierr = SEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+3, &ctx4);

  *normConst = 8.0*(integrals[2]*integrals[1] - integrals[0]*integrals[3]); flopCount += 4;

  mpfr_clears(mpfrzero, mpfrone, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "NormConstantIntFixedERF"
PetscErrorCode NormConstantIntFixedERF(EllipsoidalSystem *e, PetscInt n, PetscInt p, PetscInt prec, PetscInt nPts, PetscReal *normConst)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  FuncInfo2 ctx1 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = 1 };
  FuncInfo2 ctx2 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = 1 };
  FuncInfo2 ctx3 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = -1 };
  FuncInfo2 ctx4 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = -1 };

  PetscReal integrals[4];

  mpfr_t mpfrzero;
  mpfr_t mpfrone;
  mpfr_inits(mpfrzero, mpfrone, NULL);
  mpfr_set_d(mpfrzero, 0.0, MPFR_RNDN);
  mpfr_set_d(mpfrone, 1.0, MPFR_RNDN);

  ierr = ERFQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+0, &ctx1);
  //ierr = ERFQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) inte, mpfrzero, mpfrone, prec, nPts, integrals+0, &ctx1);
  //printf("inte: %15.15f\n", integrals[0]);
  ierr = ERFQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+1, &ctx2);

  ierr = ERFQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+2, &ctx3);

  ierr = ERFQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+3, &ctx4);

  *normConst = 8.0*(integrals[2]*integrals[1] - integrals[0]*integrals[3]); flopCount += 4;

  mpfr_clears(mpfrzero, mpfrone, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "NormPlot"
PetscErrorCode NormPlot()
{
  const PetscInt MPFR_PREC = 128;
  const PetscInt NUM_SOLUTIONS = 100;
  const PetscInt POINTS_MIN = 4;
  const PetscInt POINTS_STEP = 2;
  const PetscInt prec = 16;
  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;
  PetscReal solExact;
  PetscReal intExact[4];
  PetscReal solutions[NUM_SOLUTIONS];
  PetscReal errors   [NUM_SOLUTIONS];
  PetscReal intSols[4*NUM_SOLUTIONS];
  PetscReal intErrors[4*NUM_SOLUTIONS];
  PetscReal sol2, sol3, sol4;
  PetscLogEvent flopCounts[NUM_SOLUTIONS];
  EllipsoidalSystem e;
  PetscErrorCode ierr;
  PetscInt i;
  PetscEventPerfInfo info;
  PetscInt n = 21;
  PetscInt p = 22;
  PetscFunctionBegin;

  mpfr_set_default_prec(4*MPFR_PREC);
  initEllipsoidalSystem(&e, a, b, c, MPFR_PREC);

  /* calculate "exact" solution */
  mpfr_t exactsol;
  mpfr_init(exactsol);
  ierr = calcNormalizationMPFR(&e, n, p, &exactsol);
  ierr = calcNormalization(&e, n, p, &solExact);CHKERRQ(ierr);
  ierr = calcNormalization2(&e, n, p, intExact, &solExact);CHKERRQ(ierr);
  solExact = mpfr_get_d(exactsol, MPFR_RNDN);
  printf("old norm constant: %15.15f\n", solExact);

  /* calculate approximate solutions and record flops */
  char text[40] = "%d points";
  char sText[40];
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    sprintf(sText, text, POINTS_MIN + POINTS_STEP*i);
    ierr = PetscLogEventRegister(sText, 0, flopCounts+i);CHKERRQ(ierr);
    ierr = PetscLogEventBegin(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = NormConstantIntFixed(&e, n, p, prec, POINTS_MIN + POINTS_STEP*i, intSols+4*i, solutions+i);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    printf("new norm constant: %15.15f\n", solutions[i]);
  }

  /* calculate errors */
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    errors[i] = PetscAbsReal((solExact - solutions[i])/solExact);
    intErrors[4*i+0] = PetscAbsReal((intExact[0] - intSols[4*i+0])/intExact[0]);
    intErrors[4*i+1] = PetscAbsReal((intExact[1] - intSols[4*i+1])/intExact[1]);
    intErrors[4*i+2] = PetscAbsReal((intExact[2] - intSols[4*i+2])/intExact[2]);
    intErrors[4*i+3] = PetscAbsReal((intExact[3] - intSols[4*i+3])/intExact[3]);
    printf("errors[%d] = %15.15f\n", i, errors[i]);
  }

  FILE *fp1 = fopen("out/normInt1Prec.txt", "w");
  FILE *fp2 = fopen("out/normInt2Prec.txt", "w");
  FILE *fp3 = fopen("out/normInt3Prec.txt", "w");
  FILE *fp4 = fopen("out/normInt4Prec.txt", "w");
  FILE *fp = fopen("out/normWorkPrec.txt", "w");
  fprintf(fp, "points flops error\n");
  fprintf(fp1, "points error\n");
  fprintf(fp2, "points error\n");
  fprintf(fp3, "points error\n");
  fprintf(fp4, "points error\n");
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, flopCounts[i], &info);CHKERRQ(ierr);
    fprintf(fp, "%d %4.4e %4.4e\n", POINTS_MIN + POINTS_STEP*i, info.flops, errors[i]);
    fprintf(fp1, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+0]);
    fprintf(fp2, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+1]);
    fprintf(fp3, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+2]);
    fprintf(fp4, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+3]);
  }
  fclose(fp);
  fclose(fp1);
  fclose(fp2);
  fclose(fp3);
  mpfr_clear(exactsol);

  PetscFunctionReturn(0);
}

/*
NormPlot2 - plots convergence of a normalization
            constant for three different ellipsoidal
	    systems to compare convergence
 */

#undef __FUNCT__
#define __FUNCT__ "NormPlot2"
PetscErrorCode NormPlot2()
{
  const PetscInt NUM_SOLUTIONS = 43;
  const PetscInt POINTS_MIN = 4;
  const PetscInt POINTS_STEP = 2;
  const PetscInt prec = 16;
  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;
  PetscReal solExact;
  PetscReal intExact[4];
  PetscReal solutions[NUM_SOLUTIONS];
  PetscReal errors   [NUM_SOLUTIONS];
  PetscReal intSols[4*NUM_SOLUTIONS];
  PetscReal intErrors[4*NUM_SOLUTIONS];
  PetscReal sol2, sol3, sol4;
  PetscLogEvent flopCounts[NUM_SOLUTIONS];
  EllipsoidalSystem e1, e2, e3;
  EllipsoidalSystem *e;
  PetscErrorCode ierr;
  PetscInt i;
  PetscEventPerfInfo info;
  PetscInt n = 3;
  PetscInt p = 5;
  PetscFunctionBegin;

  ierr = initEllipsoidalSystem(&e1, 2.5, 2.0, 1.0, 64);CHKERRQ(ierr);
  ierr = initEllipsoidalSystem(&e2, 2.5, 2.0, 1.5, 64);CHKERRQ(ierr);
  ierr = initEllipsoidalSystem(&e3, 2.5, 2.0, 1.89, 64);CHKERRQ(ierr);

  for(PetscInt num=0; num<3; ++num) {
    if(num==0)
      e = &e1;
    else if(num==1)
      e = &e2;
    else
      e = &e3;
    /* calculate "exact" solution */
    ierr = calcNormalization(e, n, p, &solExact);CHKERRQ(ierr);
    ierr = calcNormalization2(e, n, p, intExact, &solExact);CHKERRQ(ierr);
    printf("old norm constant: %15.15f\n", solExact);

    /* calculate approximate solutions and record flops */
    char text[40] = "%d points";
    char sText[40];
    for(i=0; i<NUM_SOLUTIONS; ++i) {
      sprintf(sText, text, POINTS_MIN + POINTS_STEP*i);
      ierr = PetscLogEventRegister(sText, 0, flopCounts+i);CHKERRQ(ierr);
      ierr = PetscLogEventBegin(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
      ierr = NormConstantIntFixed(e, n, p, prec, POINTS_MIN + POINTS_STEP*i, intSols+4*i, solutions+i);CHKERRQ(ierr);
      ierr = PetscLogEventEnd(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
      printf("new norm constant: %15.15f\n", solutions[i]);
    }

    /* calculate errors */
    for(i=0; i<NUM_SOLUTIONS; ++i) {
      errors[i] = PetscAbsReal((solExact - solutions[i])/solExact);
      intErrors[4*i+0] = PetscAbsReal((intExact[0] - intSols[4*i+0])/intExact[0]);
      intErrors[4*i+1] = PetscAbsReal((intExact[1] - intSols[4*i+1])/intExact[1]);
      intErrors[4*i+2] = PetscAbsReal((intExact[2] - intSols[4*i+2])/intExact[2]);
      intErrors[4*i+3] = PetscAbsReal((intExact[3] - intSols[4*i+3])/intExact[3]);
      printf("errors[%d] = %15.15f\n", i, errors[i]);
    }
    FILE *fp1, *fp2, *fp3, *fp4, *fp;
    if(num==0) {
      fp1 = fopen("out/e1normInt1Prec.txt", "w");
      fp2 = fopen("out/e1normInt2Prec.txt", "w");
      fp3 = fopen("out/e1normInt3Prec.txt", "w");
      fp4 = fopen("out/e1normInt4Prec.txt", "w");
      fp = fopen("out/e1normWorkPrec.txt", "w");
    }
    else if(num==1) {
      fp1 = fopen("out/e2normInt1Prec.txt", "w");
      fp2 = fopen("out/e2normInt2Prec.txt", "w");
      fp3 = fopen("out/e2normInt3Prec.txt", "w");
      fp4 = fopen("out/e2normInt4Prec.txt", "w");
      fp = fopen("out/e2normWorkPrec.txt", "w");
    }
    else {
      fp1 = fopen("out/e3normInt1Prec.txt", "w");
      fp2 = fopen("out/e3normInt2Prec.txt", "w");
      fp3 = fopen("out/e3normInt3Prec.txt", "w");
      fp4 = fopen("out/e3normInt4Prec.txt", "w");
      fp = fopen("out/e3normWorkPrec.txt", "w");
    }
    fprintf(fp, "points flops error\n");
    fprintf(fp1, "points error\n");
    fprintf(fp2, "points error\n");
    fprintf(fp3, "points error\n");
    fprintf(fp4, "points error\n");
    for(i=0; i<NUM_SOLUTIONS; ++i) {
      ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, flopCounts[i], &info);CHKERRQ(ierr);
      fprintf(fp, "%d %4.4e %4.4e\n", POINTS_MIN + POINTS_STEP*i, info.flops, errors[i]);
      fprintf(fp1, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+0]);
      fprintf(fp2, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+1]);
      fprintf(fp3, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+2]);
      fprintf(fp4, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, intErrors[4*i+3]);
    }
    fclose(fp);
    fclose(fp1);
    fclose(fp2);
    fclose(fp3);

  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "NormPlotSE"
PetscErrorCode NormPlotSE()
{
  const PetscInt NUM_SOLUTIONS = 60;
  const PetscInt POINTS_MIN = 4;
  const PetscInt POINTS_STEP = 2;
  const PetscInt prec = 16;
  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;
  PetscReal solOld;
  PetscReal solutions[NUM_SOLUTIONS];
  PetscReal errors   [NUM_SOLUTIONS];
  PetscReal sol2, sol3, sol4;
  PetscLogEvent flopCounts[NUM_SOLUTIONS];
  EllipsoidalSystem e;
  PetscErrorCode ierr;
  PetscInt i;
  PetscEventPerfInfo info;
  PetscInt n = 3;
  PetscInt p = 5;
  PetscFunctionBegin;

  initEllipsoidalSystem(&e, a, b, c, 64);

  /* calculate approximate solutions and record flops */
  char text[40] = "%d points";
  char sText[40];
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    sprintf(sText, text, POINTS_MIN + POINTS_STEP*i);
    ierr = PetscLogEventRegister(sText, 0, flopCounts+i);CHKERRQ(ierr);
    ierr = PetscLogEventBegin(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = NormConstantIntFixedSE(&e, n, p, prec, POINTS_MIN + POINTS_STEP*i, solutions+i);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    printf("new norm constant (N=%d): %15.15f\n", POINTS_MIN+POINTS_STEP*i, solutions[i]);
  }

  /* calculate "exact" solution */
  ierr = calcNormalization(&e, n, p, &solOld);CHKERRQ(ierr);
  printf("old norm constant: %15.15f\n", solOld);

  /* calculate errors */
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    errors[i] = PetscAbsReal((solOld - solutions[i])/solOld);
    printf("errors[%d] = %15.15f\n", i, errors[i]);
  }

  FILE *fp = fopen("out/normWorkPrecSE.txt", "w");
  fprintf(fp, "points flops error\n");
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, flopCounts[i], &info);CHKERRQ(ierr);
    fprintf(fp, "%d %4.4e %4.4e\n", POINTS_MIN + POINTS_STEP*i, info.flops, errors[i]);
  }
  fclose(fp);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "NormPlotERF"
PetscErrorCode NormPlotERF()
{
  const PetscInt NUM_SOLUTIONS = 50;
  const PetscInt POINTS_MIN = 4;
  const PetscInt POINTS_STEP = 2;
  const PetscInt prec = 16;
  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;
  PetscReal solOld;
  PetscReal solutions[NUM_SOLUTIONS];
  PetscReal errors   [NUM_SOLUTIONS];
  PetscReal sol2, sol3, sol4;
  PetscLogEvent flopCounts[NUM_SOLUTIONS];
  EllipsoidalSystem e;
  PetscErrorCode ierr;
  PetscInt i;
  PetscEventPerfInfo info;
  PetscInt n = 5;
  PetscInt p = 3;
  PetscFunctionBegin;

  initEllipsoidalSystem(&e, a, b, c, 64);

  /* calculate "exact" solution */
  ierr = calcNormalization(&e, n, p, &solOld);CHKERRQ(ierr);
  printf("old norm constant: %15.15f\n", solOld);

  /* calculate approximate solutions and record flops */
  char text[40] = "%d points";
  char sText[40];
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    sprintf(sText, text, POINTS_MIN + POINTS_STEP*i);
    ierr = PetscLogEventRegister(sText, 0, flopCounts+i);CHKERRQ(ierr);
    ierr = PetscLogEventBegin(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = NormConstantIntFixedERF(&e, n, p, prec, POINTS_MIN + POINTS_STEP*i, solutions+i);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    printf("new norm constant (N=%d): %15.15f\n", POINTS_MIN+POINTS_STEP*i, solutions[i]);
  }

  /* calculate errors */
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    errors[i] = PetscAbsReal((solOld - solutions[i])/solOld);
    printf("errors[%d] = %15.15f\n", i, errors[i]);
  }

  FILE *fp = fopen("out/normWorkPrecERF.txt", "w");
  fprintf(fp, "points flops error\n");
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, flopCounts[i], &info);CHKERRQ(ierr);
    fprintf(fp, "%d %4.4e %4.4e\n", POINTS_MIN + POINTS_STEP*i, info.flops, errors[i]);
  }
  fclose(fp);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "main"
PetscErrorCode main(int argc, char **argv)
{

  const PetscInt NUM_DIGITS = 128;
  const PetscInt PRECISION  = 32;
  mpfr_t a, b;
  mpfr_t val;
  PetscErrorCode ierr;
  PetscFunctionBeginUser;
  const PetscInt nPts = 10;

  /* initializations */
  ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);
  ierr = PetscLogDefaultBegin();CHKERRQ(ierr);

  mpfr_set_default_prec(4*NUM_DIGITS);
  mpfr_inits(a, b, val, NULL);
  mpfr_set_d(a, 2.0, MPFR_RNDN);
  mpfr_set_d(b, 3.0, MPFR_RNDN);

  /*
  int fac = 1;
  for(int i=0; i<10; ++i) {
    printf("\n########################################################\n");
    ierr = DEQuad(inte, a, b, PRECISION, fac*nPts, &val, NULL);CHKERRQ(ierr);
    fac *= 2;
  }
  */

  ierr = NormPlot();CHKERRQ(ierr);
  //ierr = NormPlot2();CHKERRQ(ierr);
  //ierr = NormPlotSE();CHKERRQ(ierr);
  //ierr = NormPlotERF();CHKERRQ(ierr);

  ierr = PetscFinalize();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation