text:
// MIT License
//
// Copyright (c) 2020 SunnyCase
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
#include "object.h"
#include "utility.h"
#include <chino/io.h>
#include <chino/threading.h>
#include <gsl/gsl-lite.hpp>

namespace chino::io
{
#if defined(_MSC_VER)
#pragma section(".CHDRV$A", long, read) // Begin drivers
#pragma section(".CHDRV$C", long, read) // Drivers
#pragma section(".CHDRV$Z", long, read) // End drivers
#define EXPORT_DRIVER(x) __declspec(allocate(".CHDRV$C")) const ::chino::io::driver *CHINO_CONCAT(_drv_, __COUNTER__) = &x

#pragma section(".CHHDR$A", long, read) // Begin hardware device registrations
#pragma section(".CHHDR$C", long, read) // Hardware device registrations
#pragma section(".CHHDR$Z", long, read) // End hardware device registrations
#define EXPORT_HARDWARE(x) __declspec(allocate(".CHHDR$C")) const ::chino::io::hardware_device_registration *CHINO_CONCAT(_hdr_, __COUNTER__) = &x
#elif defined(__GNUC__)
#define EXPORT_DRIVER(x) __attribute__((unused, section(".chdrv"))) const ::chino::io::driver *CHINO_CONCAT(_drv_, __COUNTER__) = &x
#define EXPORT_HARDWARE(x) __attribute__((unused, section(".chhdr"))) const ::chino::io::hardware_device_registration *CHINO_CONCAT(_hdr_, __COUNTER__) = &x
#else
#error "Unsupported compiler"
#endif

struct driver;
struct device;
struct file;

#define DEFINE_DEV_TYPE(x, v) x = v,
enum class device_type
{
#include "device_types.def"
    COUNT
};
#undef DEFINE_DEV_TYPE

enum class driver_type
{
    hardware,
    console
};

struct hardware_device_registration
{
    const driver &drv;
};

typedef result<void, error_code> (*driver_add_device_t)(const driver &drv, const hardware_device_registration &hdr);
typedef result<void, error_code> (*driver_attach_device_t)(const driver &drv, device &bottom_dev, std::string_view args);
typedef result<file *, error_code> (*driver_open_device_t)(device &dev, std::string_view filename, create_disposition create_disp);
typedef result<void, error_code> (*driver_close_device_t)(file &file);
typedef result<size_t, error_code> (*driver_read_device_t)(file &file, gsl::span<gsl::byte> buffer);
typedef result<void, error_code> (*driver_write_device_t)(file &file, gsl::span<const gsl::byte> buffer);

struct driver_operations
{
    driver_add_device_t add_device;
    driver_attach_device_t attach_device;
    driver_open_device_t open_device;
    driver_close_device_t close_device;
    driver_read_device_t read_device;
    driver_write_device_t write_device;
};

struct driver
{
    driver_type type = driver_type::hardware;
    std::string_view name;
    driver_operations ops;
};

struct device : ob::object
{
    const driver &drv;
    device_type type;
    threading::sched_spinlock syncroot;

    template <class T = uint8_t>
    T &extension() noexcept { return *reinterpret_cast<T *>(reinterpret_cast<uint8_t *>(this) + sizeof(device)); }
};

struct device_extension
{
    device &dev() noexcept { return *reinterpret_cast<device *>(reinterpret_cast<uint8_t *>(this) - sizeof(device)); }
};

struct file : ob::object
{
    device &dev;
    size_t offset;
    threading::sched_spinlock syncroot;

    template <class T = uint8_t>
    T &extension() noexcept { return *reinterpret_cast<T *>(reinterpret_cast<uint8_t *>(this) + sizeof(file)); }
};

struct file_extension
{
    io::file &file() noexcept { return *reinterpret_cast<io::file *>(reinterpret_cast<uint8_t *>(this) - sizeof(device)); }
};

result<device *, error_code> create_device(const driver &drv, device_type type, size_t extension_size) noexcept;
result<device *, error_code> create_device(std::string_view name, const driver &drv, device_type type, size_t extension_size) noexcept;
result<file *, error_code> create_file(device &dev, size_t extension_size) noexcept;

result<file *, error_code> open_file(device &dev, access_mask access, std::string_view filename, create_disposition create_disp = create_disposition::open_existing) noexcept;
result<size_t, error_code> read_file(file &file, gsl::span<gsl::byte> buffer) noexcept;
result<void, error_code> write_file(file &file, gsl::span<const gsl::byte> buffer) noexcept;
result<void, error_code> close_file(file &file) noexcept;
}

max_stars_repo_name: chino-os/chino-os

max_issues_repo_name: dotnetGame/chino-os

max_forks_repo_name: dotnetGame/chino-os


text:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//
// 			     //////////////////////////////////////////////////////////
//  			     //							     //
// 			     //   	        hybridMANTIS v1.0		     //
// 			     //                fastDETECT2 - C code                  //
//			     //		   (optical photons transport)		     //
//  			     //							     //
//                           //               used for Load Balancing                //
//			     //							     //
//			     //////////////////////////////////////////////////////////
//
//
//
//
// ****Disclaimer****
//  This software and documentation (the "Software") were developed at the Food and Drug Administration (FDA) by employees of the Federal Government in
//  the course of their official duties. Pursuant to Title 17, Section 105 of the United States Code, this work is not subject to copyright protection
//  and is in the public domain. Permission is hereby granted, free of charge, to any person obtaining a copy of the Software, to deal in the Software
//  without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, or sell copies of the
//  Software or derivatives, and to permit persons to whom the Software is furnished to do so. FDA assumes no responsibility whatsoever for use by other
//  parties of the Software, its source code, documentation or compiled executables, and makes no guarantees, expressed or implied, about its quality,
//  reliability, or any other characteristic. Further, use of this code in no way implies endorsement by the FDA or confers any advantage in regulatory
//  decisions. Although this software can be redistributed and/or modified freely, we ask that any derivative works bear some notice that they are
//  derived from it, and any modified versions bear some notice that they have been modified.
//
//	Detailed comments are available in "hybridMANTIS_cuda_ver1_0.cu" and "hybridMANTIS_c_ver1_0.c" files.
//
//	Associated publication: Sharma Diksha, Badal Andreu and Badano Aldo, "hybridMANTIS: a CPU-GPU Monte Carlo method for modeling indirect x-ray detectors with
//				columnar scintillators". Physics in Medicine and Biology, 57(8), pp. 2357–2372 (2012)
//
//
//	File:   	hybridMANTIS_c_ver1_0_LB.c
//	Author: 	Diksha Sharma (US Food and Drug Administration)
//	Email: 		diksha.sharma@fda.hhs.gov
//	Last updated:  	Apr 13, 2012
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////
//
//      Header libraries
//
/////////////////////////////////////////

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

/////////////////////////////////////////
//
//      Global variables
//
/////////////////////////////////////////

#define max_photon_per_EDE 900000	// maximum number of optical photons that can be generated per energy deposition event (EDE)

#ifndef USING_CUDA
	#define mybufsizeT 2304000	// CPU buffer size: # of events sent to the CPU
#endif

/////////////////////////////////////////
//
//      Include kernel program
//
/////////////////////////////////////////
#include "kernel_cuda_c_ver1_0_LB.cu"

////////////////////////////////////////////////////////////////////////////
//				MAIN PROGRAM			          //
////////////////////////////////////////////////////////////////////////////

#ifndef USING_CUDA

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// cpuoptlb():   Performs optical transport for finding optimal load using load balancing in the CPU
//	  	 Input arguments: cptime, cpubufsize
//
//		 cptime: 	time taken by GPU to call this routine
//		 cpubufsize:  	CPU buffer size defined by user in PENELOPE tally code. Only to be used for load balancing. This is different from 'mybufsize'.
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	void cpuoptlb_(double *cptime, int *cpubufsize)
	{

		float dcos[3]={0}; 		// directional cosines
		float normal[3]={0}; 		// normal to surface in case of TIR
		float pos[3] = {0}; 		// new position
		float old_pos[3] = {0};  	// source coordinates

		// command line arguments
		float xdetector, ydetector, radius, height, n_C, n_IC, top_absfrac, bulk_abscoeff, beta, d_min, lbound_x, lbound_y, ubound_x, ubound_y, d_max, yield, sensorRefl;
		int pixelsize, num_primary, min_optphotons, max_optphotons, num_bins;

		int nbytes = (*cpubufsize)*sizeof(struct start_info);
		struct start_info *structa;
		structa = (struct start_info*) malloc(nbytes);
		if( structa == NULL )
			printf("\n Struct start_info array CANNOT BE ALLOCATED - %d !!", (*cpubufsize));

		// get cpu time
		clock_t start, end;
		double num_sec;

		// get current time stamp to initialize seed input for RNG
		time_t seconds;
		seconds = time (NULL);
		struct timeval tv;
		gettimeofday(&tv,NULL);

		float rr=0.0f, theta=0.0f;
		float r=0.0f;		// random number
		float norm=0.0f;
		int jj=0;

		// initialize random number generator (RANECU)
		int seed_input = 271828182 ; // ranecu seed input
		int seed[2];

		// gsl variables
		const gsl_rng_type * Tgsl;
		gsl_rng * rgsl;
		double mu_gsl;
		int my_index=0;
		int result_algo = 0;
		unsigned long long int *num_rebound;

		// output image variables
		int xdim = 0;
		int ydim = 0;
		int indexi=0, indexj=0;

	      	// create a generator chosen by the environment variable GSL_RNG_TYPE
	       	gsl_rng_env_setup();
	       	Tgsl = gsl_rng_default;
	       	rgsl = gsl_rng_alloc (Tgsl);

 		// copy to local variables from PENELOPE buffers
		xdetector = inputargs_.detx;		// x dimension of detector (in um). x in (0,xdetector)
		ydetector = inputargs_.dety;		// y dimension of detector (in um). y in (0,ydetector)
		height = inputargs_.detheight;		// height of column and thickness of detector (in um). z in range (-H/2, H/2)
		radius = inputargs_.detradius;		// radius of column (in um).
		n_C = inputargs_.detnC;			// refractive index of columns
		n_IC = inputargs_.detnIC;		// refractive index of intercolumnar material
		top_absfrac = inputargs_.dettop;	// column's top surface absorption fraction (0.0, 0.5, 0.98)
		bulk_abscoeff = inputargs_.detbulk;	// column's bulk absorption coefficient (in um^-1) (0.001, 0.1 cm^-1)
		beta = inputargs_.detbeta;		// roughness coefficient of column walls
		d_min = inputargs_.detdmin;		// minimum distance a photon can travel when transmitted from a column
		d_max = inputargs_.detdmax;
		lbound_x = inputargs_.detlboundx;	// x lower bound of region of interest of output image (in um)
		lbound_y = inputargs_.detlboundy;	// y lower bound (in um)
		ubound_x = inputargs_.detuboundx;	// x upper bound (in um)
		ubound_y = inputargs_.detuboundy;	// y upper bound (in um)
		yield = inputargs_.detyield;		// yield (/eV)
		pixelsize = inputargs_.detpixel;	// 1 pixel = pixelsize microns (in um)
		sensorRefl = inputargs_.detsensorRefl;	// Non-Ideal sensor reflectivity (%)
		num_primary = inputargs_.mynumhist;	// total number of primaries to be simulated
		min_optphotons = inputargs_.minphotons;	// minimum number of optical photons detected to be included in PHS
		max_optphotons = inputargs_.maxphotons;	// maximum number of optical photons detected to be included in PHS
		num_bins = inputargs_.mynumbins;	// number of bins for genrating PHS

		// dimensions of PRF image
		xdim = ceil((ubound_x - lbound_x)/pixelsize);
		ydim = ceil((ubound_y - lbound_y)/pixelsize);
		unsigned long long int myimage[xdim][ydim];

		// memory for storing histogram of # photons detected/primary
		int *h_num_detected_prim = 0;
		h_num_detected_prim = (int*)malloc(sizeof(int)*num_primary);

		for(indexj=0; indexj < num_primary; indexj++)
		  h_num_detected_prim[indexj] = 0;

		// start the clock
		start = clock();

	for(my_index = 0; my_index < (*cpubufsize); my_index++)		// iterate over x-rays
	{

		// reset the global counters
		num_generatedT = 0;
		num_detectT=0;
		num_abs_topT=0;
		num_abs_bulkT=0;
		num_lostT=0;
		num_outofcolT=0;
		num_theta1T=0;
		photon_distanceT=0.0f;

		// copying fortran buffer into *structa

		// units in the penelope output file are in cm. Convert to microns.
		structa[my_index].str_x = optical_.xbufopt[my_index] * 10000.0f;	// x-coordinate of interaction event.
		structa[my_index].str_y = optical_.ybufopt[my_index] * 10000.0f;	// y-coordinate
		structa[my_index].str_z = optical_.zbufopt[my_index] * 10000.0f;	// z-coordinate
		structa[my_index].str_E = optical_.debufopt[my_index];			// energy deposited
		structa[my_index].str_histnum = optical_.nbufopt[my_index];		// x-ray history number

		// sample # optical photons based on light yield and energy deposited for this interaction event (using Poisson distribution)
		mu_gsl = (double)structa[my_index].str_E * yield;
		structa[my_index].str_N = gsl_ran_poisson(rgsl,mu_gsl);

		if(structa[my_index].str_N > max_photon_per_EDE)
		{
			printf("\n\n CPU str_n exceeds max photons. program is exiting - %d !! \n\n",structa[my_index].str_N);
			exit(0);
		}

		num_rebound = (unsigned long long int*) malloc(structa[my_index].str_N*sizeof(unsigned long long int));
		if(num_rebound == NULL)
			printf("\n Error allocating num_rebound memory !\n");

		// initialize the RANECU generator in a position far away from the previous history:
		seed_input = (int)(seconds/3600+tv.tv_usec);			// seed input=seconds passed since 1970+current time in micro secs
		init_PRNG(my_index, 50000, seed_input, seed);      		// intialize RNG

		for(jj=0; jj<structa[my_index].str_N; jj++)
			num_rebound[jj] = 0;

		// reset the vectors
		dcos[0]=0.0f; dcos[1]=0.0f; dcos[2]=0.0f;
		normal[0]=0.0f; normal[1]=0.0f; normal[2]=0.0f;

		// set starting location of photon
		pos[0] = structa[my_index].str_x; pos[1] = structa[my_index].str_y; pos[2] = structa[my_index].str_z;
		old_pos[0] = structa[my_index].str_x; old_pos[1] = structa[my_index].str_y; old_pos[2] = structa[my_index].str_z;

		// initializing the direction cosines for the first particle in each core
		r = (ranecu(seed) * 2.0f) - 1.0f; // random number between (-1,1)

		while(fabs(r) <= 0.01f)
		 {
		   	r = (ranecu(seed) * 2.0f) - 1.0f;
		 }

		dcos[2] = r;		// random number between (-1,1)
		rr = sqrt(1.0f-r*r);
		theta=ranecu(seed)*twopipen;
		dcos[0]=rr*cos(theta);
		dcos[1]=rr*sin(theta);

		norm = sqrt(dcos[0]*dcos[0] + dcos[1]*dcos[1] + dcos[2]*dcos[2]);

		if ((norm < (1.0f - epsilon)) || (norm > (1.0f + epsilon)))	// normalize
		 {
			dcos[0] = dcos[0]/norm;
			dcos[1] = dcos[1]/norm;
			dcos[2] = dcos[2]/norm;
		 }

		local_counterT=0;
		while(local_counterT < structa[my_index].str_N)
		 {

			absorbedT = 0;
			detectT = 0;
			bulk_absT = 0;

			// set starting location of photon
			pos[0] = structa[my_index].str_x; pos[1] = structa[my_index].str_y; pos[2] = structa[my_index].str_z;
			old_pos[0] = structa[my_index].str_x; old_pos[1] = structa[my_index].str_y; old_pos[2] = structa[my_index].str_z;
			num_generatedT++;
			result_algo = 0;

			while(result_algo == 0)
			 {
			  	result_algo = algoT(normal, old_pos, pos, dcos, num_rebound, seed, structa[my_index], &myimage[0][0], xdetector, ydetector, radius, height, n_C, n_IC, top_absfrac, bulk_abscoeff, beta, d_min, pixelsize, lbound_x, lbound_y, ubound_x, ubound_y, sensorRefl, d_max, ydim, h_num_detected_prim);
			 }

		 }

		// release resources
		free(num_rebound);

	}	// my_index loop ends

	// end the clock
	end = clock();

	num_sec = (double)((end - start)/CLOCKS_PER_SEC);
        *cptime = num_sec;

	// release resources
	free(structa);
	free(h_num_detected_prim);

		return;
	}	// C main() ends

#endif

ext: c

max_stars_repo_name: diamfda/hybridmantis

max_issues_repo_name: diamfda/hybridmantis

max_forks_repo_name: diamfda/hybridmantis


text:
#include "coneOS.h"
#include <cblas.h>
#include <math.h>

// x = b*a
inline void setAsScaledArray(double *x, const double * a,const double b,int len) {
  int i;
  for( i=0;i<len;++i ) x[i] = b*a[i];
}

// a*= b
inline void scaleArray(double * a,const double b,int len){
  //b_dscal(len,b,a,1);
  cblas_dscal(len,b,a,1);
}

// x'*y
inline double innerProd(const double * x, const double * y, int len){
  //return b_ddot(len, x,1,y,1);
  return cblas_ddot(len, x,1,y,1);
}

// ||v||_2
inline double calcNorm(const double * v,int len){
  //return b_dnrm2(len, v, 1);
  return cblas_dnrm2(len, v, 1);
}

inline double calcNormInf(const double *a, int l){
	double tmp, max = 0.0;
	int i;
	for ( i=0; i<l; ++i){
		tmp = fabs(a[i]);
		if(tmp > max) max = tmp;
	}
	return max;
}

// ||v||_2^2
inline double calcNormSq(const double * v,int len){
  double nrm = calcNorm(v,len);
  return nrm*nrm;
}

// daxpy a += sc*b
inline void addScaledArray(double * a, const double * b, int n, const double sc){
  //b_daxpy(n, sc, b, 1, a, 1);
  cblas_daxpy(n, sc, b, 1, a, 1);
}

inline double calcNormDiff(const double *a,const double *b, int l) {
  double nmDiff = 0.0, tmp;
  int i;
  for ( i=0; i<l; ++i ){
    tmp = (a[i] - b[i]);
    nmDiff += tmp * tmp;
  }
  return sqrt(nmDiff);
}

inline double calcNormInfDiff(const double *a, const double *b, int l) {
	double tmp, max = 0.0;
	int i;
	for ( i=0; i<l; ++i){
		tmp = fabs(a[i] - b[i]);
		if(tmp > max) max = tmp;
	}
	return max;
}

inline void accumByAtrans(Data * d, const double *x, double *y)
{
	cblas_dgemv(CblasColMajor, CblasTrans, d->m, d->n, 1.0, d->Ax, d->m, x, 1, 1.0, y, 1);
}

inline void accumByA(Data * d, const double *x, double *y)
{
 	cblas_dgemv(CblasColMajor, CblasNoTrans, d->m, d->n, 1.0, d->Ax, d->m, x, 1, 1.0, y, 1);
}

/*

inline void b_daxpy(const int n, const double alpha, const double *x, const int incx, double *y, const int incy)
{
  daxpy(&n,&alpha,x,&incx,y,&incy);
  //cblas_daxpy(n,alpha,x,incx,y,incy);
}

inline void b_dsyr(char Uplo, const int N, const double alpha, const double *X, const int incX, double *A, const int lda)
{
  //dsyr(&Uplo,&N,&alpha,X,&incX,A,&lda);
  cblas_dsyr(CblasColMajor, Uplo, N, alpha, X, incX, A, lda);
}

inline void b_dscal(const int N, const double alpha, double *X, const int incX)
{
  //dscal(&N,&alpha,X,&incX);
  cblas_dscal(N,alpha,X,incX);
}

inline double b_ddot(const int n, const double *x, const int incx, const double *y, const int incy)
{
  //return ddot(&n,x,&incx,y,&incy);
  return cblas_ddot(n,x,incx,y,incy);
}

inline double b_dnrm2(const int N, const double *X, const int incX)
{
  return cblas_dnrm2(N,X,incX);
  //return dnrm2(&N,X,&incX);
}

inline void b_dgemm(
              char TransA,
              char TransB,
              const int M,
              const int N,
              const int K,
              const double alpha,
              const double *A,
              const int lda,
              const double *B,
              const int ldb,
              const double beta,
              double *C,
              const int ldc){
  dgemm(&TransA, &TransB, &M, &N, &K,&alpha,A,&lda,B,&ldb,&beta,C,&ldc);
  //cblas_dgemm();
}

inline void b_dgemv(char Trans,
                    const int m,
                    const int n,
                    const double alpha,
                    const double  *a,
                    const int lda,
                    const double  *x,
                    const int incx,
                    const double beta,
                    double  *y,
                    const int incy){
  dgemv(&Trans,&m,&n,&alpha,a,&lda,x,&incx,&beta,y,&incy);
}

inline void b_dtrsv(char Uplo, char TransA, char Diag, const int N, const double *A, const int lda, double *X, const int incX)
{
  dtrsv(&Uplo, &TransA, &Diag, &N, A, &lda, X, &incX);
}

inline void b_dsymv(char Uplo, const int N, const double alpha, const double *A,
                    const int lda, const double *X, const int incX, const double beta,
                    double *Y, const int incY){
  dsymv(&Uplo, &N, &alpha, A, &lda, X, &incX, beta, Y, &incY);
}

*/

max_stars_repo_name: cvxgrp/coneos

max_issues_repo_name: cvxgrp/coneos

max_forks_repo_name: cvxgrp/coneos


text:
/*
FRACTAL - A program growing fractals to benchmark parallelization and drawing
libraries.

Copyright 2009-2021, Javier Burguete Tolosa.

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 Javier Burguete Tolosa ``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 Javier Burguete Tolosa 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.
*/

/**
 * \file graphic.c
 * \brief Source file to define the graphic drawing data and functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2009-2021, Javier Burguete Tolosa.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#include <png.h>
#include <ft2build.h>
#include FT_FREETYPE_H
#include <gtk/gtk.h>
#include <GL/glew.h>

#include "config.h"
#include "fractal.h"
#include "image.h"
#include "text.h"
#include "graphic.h"

unsigned int window_width = 480;        ///< Graphic window width.
unsigned int window_height = 480;       ///< Graphic window height.

/**
 * Function to init the graphic data.

 * \return 1 on success, 0 on error.
 */
int
graphic_init (Graphic * graphic,        ///< Draw struct.
              char *logo_name)  ///< Logo PNG file name.
{
  const GLfloat projection_matrix[16] = {
    1., 0., 0., 0.,
    0, 1., 0., 0.,
    0., 0, 0., 0.,
    0., 0., 0., 1.
  };
  const char *vs_3D_source =
    "attribute highp vec3 position;"
    "attribute lowp vec3 color;"
    "varying lowp vec3 fcolor;"
    "uniform highp mat4 matrix;"
    "void main ()"
    "{gl_Position = matrix * vec4 (position, 1.f); fcolor = color;}";
  const char *fs_source =
    "varying lowp vec3 fcolor;"
    "void main () {gl_FragColor = vec4 (fcolor, 1.f);}";
  const char *vertex_name = "position";
  const char *color_name = "color";
  const char *matrix_name = "matrix";
  // GLSL version
  const char *version = "#version 120\n";       // OpenGL 2.1
  const char *vs_3D_sources[3] = { version, INIT_GL_GLES, vs_3D_source };
  const char *fs_sources[3] = { version, INIT_GL_GLES, fs_source };
  const char *error_message;
  GLint k;
  GLuint vs, fs;
  GLenum glew_status;

#if DEBUG
  printf ("graphic_init: start\n");
  fflush (stdout);
#endif

  // Initing GLEW
#if DEBUG
  printf ("Initing GLEW\n");
  fflush (stdout);
#endif
  glew_status = glewInit ();
  if (glew_status != GLEW_OK)
    {
      printf ("ERROR! glewInit: %s\n", glewGetErrorString (glew_status));
      return 0;
    }

#if DEBUG
  printf ("graphic_init: compiling fragment shader\n");
  fflush (stdout);
#endif
  fs = glCreateShader (GL_FRAGMENT_SHADER);
  glShaderSource (fs, 3, fs_sources, NULL);
  glCompileShader (fs);
  glGetShaderiv (fs, GL_COMPILE_STATUS, &k);
  if (!k)
    {
      error_message = "unable to compile the fragment shader";
      goto exit_on_error;
    }

#if DEBUG
  printf ("graphic_init: compiling 3D vertex shader\n");
  fflush (stdout);
#endif
  vs = glCreateShader (GL_VERTEX_SHADER);
  glShaderSource (vs, 3, vs_3D_sources, NULL);
  glCompileShader (vs);
  glGetShaderiv (vs, GL_COMPILE_STATUS, &k);
  if (!k)
    {
      glDeleteShader (fs);
      error_message = "unable to compile the 3D vertex shader";
      goto exit_on_error;
    }

  graphic->program_3D = glCreateProgram ();
  glAttachShader (graphic->program_3D, vs);
  glAttachShader (graphic->program_3D, fs);
  glLinkProgram (graphic->program_3D);
  glDeleteShader (vs);
  glDeleteShader (fs);
  glGetProgramiv (graphic->program_3D, GL_LINK_STATUS, &k);
  if (!k)
    {
      error_message = "unable to link the program 3D";
      goto exit_on_error;
    }

  graphic->attribute_3D_position
    = glGetAttribLocation (graphic->program_3D, vertex_name);
  if (graphic->attribute_3D_position == -1)
    {
      error_message = "could not bind attribute";
      goto exit_on_error;
    }

  graphic->attribute_3D_color
    = glGetAttribLocation (graphic->program_3D, color_name);
  if (graphic->attribute_3D_color == -1)
    {
      error_message = "could not bind attribute";
      goto exit_on_error;
    }

  graphic->uniform_3D_matrix
    = glGetUniformLocation (graphic->program_3D, matrix_name);
  if (graphic->uniform_3D_matrix == -1)
    {
      error_message = "could not bind uniform";
      goto exit_on_error;
    }

  memcpy (graphic->projection_matrix, projection_matrix, 16 * sizeof (GLfloat));

#if DEBUG
  printf ("graphic_init: initing logo\n");
  fflush (stdout);
#endif
  if (logo_name)
    {
      graphic->logo = image_new (logo_name);
      if (!graphic->logo)
        {
          error_message = "unable to open the logo";
          goto exit_on_error;
        }
      if (!image_init (graphic->logo))
        {
          error_message = "unable to init the logo";
          goto exit_on_error;
        }
    }
  else
    graphic->logo = NULL;

#if DEBUG
  printf ("graphic_init: initing text\n");
  fflush (stdout);
#endif
  if (!text_init (graphic->text))
    {
      error_message = "unable to init the text";
      goto exit_on_error;
    }

#if DEBUG
  printf ("graphic_init: end\n");
  fflush (stdout);
#endif
  return 1;

exit_on_error:
  printf ("ERROR! %s\n", error_message);
#if DEBUG
  printf ("graphic_init: end\n");
  fflush (stdout);
#endif
  return 0;
}

/**
 * Function to free the memory used by graphic drawing functions.
 */
void
graphic_destroy (Graphic * graphic)     ///< Graphic struct.
{
#if DEBUG
  printf ("graphic_destroy: start\n");
  fflush (stdout);
#endif
  text_destroy (graphic->text);
  image_destroy (graphic->logo);
#if DEBUG
  printf ("graphic_destroy: end\n");
  fflush (stdout);
#endif
}

/**
 * Function to draw the fractal.
 */
void
graphic_render (Graphic * graphic)      ///< Graphic struct.
{
  // Rectangle matrix
  const Point3D square_vertices[4] = {
    {{0., 0., 0.}, {0., 0., 0.}},
    {{length, 0., 0.}, {0., 0., 0.}},
    {{length, width, 0.}, {0., 0., 0.}},
    {{0., width, 0.}, {0., 0., 0.}}
  };
  const GLushort square_indices[4] = { 0, 1, 2, 3 };
  const GLfloat black[4] = { 0., 0., 0., 1. };
  const char *str_version = "Fractal 3.4.15";
  float cp, sp, ct, st, w, h, sx, sy;
  GLuint vbo_square, ibo_square, vbo_points;

#if DEBUG
  printf ("graphic_render: start\n");
  fflush (stdout);
#endif

  // Drawing a white background
  glClearColor (1., 1., 1., 0.);
  glClear (GL_COLOR_BUFFER_BIT);

#if DEBUG
  printf ("graphic_render: ending if no fractal\n");
  fflush (stdout);
#endif

  // Ending if no fractal
  if (!medium)
    goto end_draw;

#if DEBUG
  printf ("graphic_render: checking the fractal type\n");
  fflush (stdout);
#endif

  // Checking if 3D or 2D fractal
  glUseProgram (graphic->program_3D);
  if (fractal_3D)
    {
      // Projection matrix
      sincosf (graphic->phi, &sp, &cp);
      sincosf (graphic->theta, &st, &ct);
      w = graphic->xmax - graphic->xmin;
      h = graphic->ymax - graphic->ymin;
      graphic->projection_matrix[0] = 2. * cp / w;
      graphic->projection_matrix[1] = 2. * ct * sp / h;
      graphic->projection_matrix[4] = -2. * sp / w;
      graphic->projection_matrix[5] = 2. * ct * cp / h;
      graphic->projection_matrix[9] = 2. * st / h;
      graphic->projection_matrix[12] = -1.;
      graphic->projection_matrix[13] = -1. - 2. * graphic->ymin / h;
      glUniformMatrix4fv (graphic->uniform_3D_matrix, 1, GL_FALSE,
                          graphic->projection_matrix);

      // Drawing a black rectangle
      glGenBuffers (1, &vbo_square);
      glBindBuffer (GL_ARRAY_BUFFER, vbo_square);
      glBufferData (GL_ARRAY_BUFFER, sizeof (square_vertices), square_vertices,
                    GL_DYNAMIC_DRAW);
      glGenBuffers (1, &ibo_square);
      glBindBuffer (GL_ELEMENT_ARRAY_BUFFER, ibo_square);
      glBufferData (GL_ELEMENT_ARRAY_BUFFER, sizeof (square_indices),
                    square_indices, GL_DYNAMIC_DRAW);
      glBindBuffer (GL_ARRAY_BUFFER, vbo_square);
      glEnableVertexAttribArray (graphic->attribute_3D_position);
      glVertexAttribPointer (graphic->attribute_3D_position,
                             3, GL_FLOAT, GL_FALSE, sizeof (Point3D), NULL);
      glEnableVertexAttribArray (graphic->attribute_3D_color);
      glVertexAttribPointer (graphic->attribute_3D_color,
                             3,
                             GL_FLOAT,
                             GL_FALSE,
                             sizeof (Point3D),
                             (GLvoid *) offsetof (Point3D, c));
      glBindBuffer (GL_ELEMENT_ARRAY_BUFFER, ibo_square);
      glDrawElements (GL_LINE_LOOP, 4, GL_UNSIGNED_SHORT, 0);
      glDeleteBuffers (1, &ibo_square);
      glDeleteBuffers (1, &vbo_square);
      glDisableVertexAttribArray (graphic->attribute_3D_color);
      glDisableVertexAttribArray (graphic->attribute_3D_position);
    }
  else
    {
      // Projection matrix
      graphic->projection_matrix[0] = 2. / width;
      graphic->projection_matrix[5] = 2. / height;
      graphic->projection_matrix[12] = graphic->projection_matrix[13] = -1.;
      graphic->projection_matrix[1] = graphic->projection_matrix[4]
        = graphic->projection_matrix[9] = 0.;
      glUniformMatrix4fv (graphic->uniform_3D_matrix, 1, GL_FALSE,
                          graphic->projection_matrix);
    }

#if DEBUG
  printf ("graphic_render: drawing the fractal points\n");
  fflush (stdout);
#endif

  // Drawing the fractal points
  glGenBuffers (1, &vbo_points);
  glBindBuffer (GL_ARRAY_BUFFER, vbo_points);
  glBufferData (GL_ARRAY_BUFFER, npoints * sizeof (Point3D), point,
                GL_DYNAMIC_DRAW);
  glEnableVertexAttribArray (graphic->attribute_3D_position);
  glVertexAttribPointer (graphic->attribute_3D_position,
                         3, GL_FLOAT, GL_FALSE, sizeof (Point3D), NULL);
  glEnableVertexAttribArray (graphic->attribute_3D_color);
  glVertexAttribPointer (graphic->attribute_3D_color,
                         3,
                         GL_FLOAT,
                         GL_FALSE,
                         sizeof (Point3D), (GLvoid *) offsetof (Point3D, c));
  glDrawArrays (GL_POINTS, 0, npoints);
  glDeleteBuffers (1, &vbo_points);
  glDisableVertexAttribArray (graphic->attribute_3D_color);
  glDisableVertexAttribArray (graphic->attribute_3D_position);

end_draw:

  // OpenGL properties
  glEnable (GL_BLEND);
  glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

#if DEBUG
  printf ("graphic_render: drawing the logo\n");
  fflush (stdout);
#endif

  // Drawing the logo
  image_draw (graphic->logo, 0, 0, window_width, window_height);

#if DEBUG
  printf ("graphic_render: displaying the program version\n");
  fflush (stdout);
#endif

  // Displaying the program version
  sx = 0.15 * 12. / window_width;
  sy = 0.15 * 12. / window_height;
  text_draw (graphic->text, (char *) str_version,
             0.99 - 7. * strlen (str_version) * sx, -0.99, sx, sy, black);

  // Disabling OpenGL properties
  glDisable (GL_BLEND);

#if DEBUG
  printf ("graphic_render: displaying the draw\n");
  fflush (stdout);
#endif

#if DEBUG
  printf ("graphic_render: end\n");
  fflush (stdout);
#endif
}

/**
 * Function to save the draw in a PNG file.
 */
void
graphic_save (char *file_name)  ///< File name.
{
  png_struct *png;
  png_info *info;
  png_byte **row_pointers;
  GLubyte *pixels;
  FILE *file;
  unsigned int i, row_bytes, pointers_bytes, pixels_bytes;

#if DEBUG
  printf ("graphic_save: start\n");
  fflush (stdout);
#endif

  // Creating the PNG header
  png = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
  if (!png)
    return;
  info = png_create_info_struct (png);
  if (!info)
    return;
  file = fopen (file_name, "wb");
  if (!file)
    return;
  if (setjmp (png_jmpbuf (png)))
    {
      printf ("Error png_init_io\n");
      exit (0);
    }
  png_init_io (png, file);
  if (setjmp (png_jmpbuf (png)))
    {
      printf ("Error png_set_IHDR\n");
      exit (0);
    }
  png_set_IHDR (png,
                info,
                window_width,
                window_height,
                8,
                PNG_COLOR_TYPE_RGBA,
                PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
  if (setjmp (png_jmpbuf (png)))
    {
      printf ("Error png_write_info\n");
      exit (0);
    }
  png_write_info (png, info);

  // Getting the OpenGL pixels
  glViewport (0, 0, window_width, window_height);
  row_bytes = 4 * window_width;
  pixels_bytes = row_bytes * window_height;
  pixels = (GLubyte *) g_slice_alloc (pixels_bytes);
  glReadPixels (0, 0, window_width, window_height, GL_RGBA, GL_UNSIGNED_BYTE,
                pixels);

  // Saving the pixels in the PNG order
  pointers_bytes = window_height * sizeof (png_byte *);
  row_pointers = (png_byte **) g_slice_alloc (pointers_bytes);
  for (i = 0; i < window_height; ++i)
    {
      row_pointers[i] = (png_byte *) g_slice_alloc (row_bytes);
      memcpy (row_pointers[i], pixels + (window_height - 1 - i) * row_bytes,
              row_bytes);
    }
  if (setjmp (png_jmpbuf (png)))
    {
      printf ("Error png_write_image\n");
      exit (0);
    }
  png_write_image (png, row_pointers);

  // Freeing memory
  for (i = 0; i < window_height; ++i)
    g_slice_free1 (row_bytes, row_pointers[i]);
  g_slice_free1 (pointers_bytes, row_pointers);
  g_slice_free1 (pixels_bytes, pixels);

  // Saving the file
  if (setjmp (png_jmpbuf (png)))
    {
      printf ("Error png_write_end\n");
      exit (0);
    }
  png_write_end (png, NULL);
  fclose (file);

  // Freeing memory
  png_destroy_write_struct (&png, &info);

#if DEBUG
  printf ("graphic_save: end\n");
  fflush (stdout);
#endif
}

max_stars_repo_name: jburguete/fractal

max_issues_repo_name: jburguete/fractal

max_forks_repo_name: jburguete/fractal


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 <nlohmann/json_fwd.hpp>
#include <gsl/gsl>

#include <string>

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

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

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

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

} // namespace cb::breakpad

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

ext: h


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <memory>
#include "halley/data_structures/vector.h"
#include <gsl/gsl>

#include "halley/core/api/audio_api.h"

struct OggVorbis_File;

#if defined(_MSC_VER) || defined(__clang__)
using OggOffsetType = int64_t;
#else
using OggOffsetType = long int;
#endif

namespace Halley {
	class ResourceData;
	class ResourceDataReader;

	class VorbisData {
	public:
		VorbisData(std::shared_ptr<ResourceData> resource, bool open);
		~VorbisData();

		size_t read(gsl::span<Vector<float>> dst);
		size_t read(AudioMultiChannelSamples dst, size_t nChannels);

		size_t getNumSamples() const; // Per channel
		int getSampleRate() const;
		int getNumChannels() const;

		void close();
		void reset();
		void seek(double t);
		void seek(size_t sample);
		size_t tell() const;

		size_t getSizeBytes() const;

	private:
		void open();
		static size_t vorbisRead(void* ptr, size_t size, size_t nmemb, void* datasource);
		static int vorbisSeek(void *datasource, OggOffsetType offset, int whence);
		static int vorbisClose(void *datasource);
		static long vorbisTell(void *datasource);

		std::shared_ptr<ResourceData> resource;
		std::shared_ptr<ResourceDataReader> stream;
		OggVorbis_File* file;

		bool streaming;
		long long pos;
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


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

#ifndef __GSL_MATRIX_SHORT_H__
#define __GSL_MATRIX_SHORT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_short.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  short * data;
  gsl_block_short * block;
  int owner;
} gsl_matrix_short;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_view;

typedef _gsl_matrix_short_view gsl_matrix_short_view;

typedef struct
{
  gsl_matrix_short matrix;
} _gsl_matrix_short_const_view;

typedef const _gsl_matrix_short_const_view gsl_matrix_short_const_view;

/* Allocation */

GSL_FUN gsl_matrix_short *
gsl_matrix_short_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_short *
gsl_matrix_short_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_short *
gsl_matrix_short_alloc_from_block (gsl_block_short * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_short *
gsl_matrix_short_alloc_from_matrix (gsl_matrix_short * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_short *
gsl_vector_short_alloc_row_from_matrix (gsl_matrix_short * m,
                                        const size_t i);

GSL_FUN gsl_vector_short *
gsl_vector_short_alloc_col_from_matrix (gsl_matrix_short * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_short_free (gsl_matrix_short * m);

/* Views */

GSL_FUN _gsl_matrix_short_view
gsl_matrix_short_submatrix (gsl_matrix_short * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_row (gsl_matrix_short * m, const size_t i);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_column (gsl_matrix_short * m, const size_t j);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_diagonal (gsl_matrix_short * m);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_subdiagonal (gsl_matrix_short * m, const size_t k);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_superdiagonal (gsl_matrix_short * m, const size_t k);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_subrow (gsl_matrix_short * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_short_view
gsl_matrix_short_subcolumn (gsl_matrix_short * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_short_view
gsl_matrix_short_view_array (short * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_short_view
gsl_matrix_short_view_array_with_tda (short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_short_view
gsl_matrix_short_view_vector (gsl_vector_short * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_short_view
gsl_matrix_short_view_vector_with_tda (gsl_vector_short * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_short_const_view
gsl_matrix_short_const_submatrix (const gsl_matrix_short * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_row (const gsl_matrix_short * m,
                            const size_t i);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_column (const gsl_matrix_short * m,
                               const size_t j);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_diagonal (const gsl_matrix_short * m);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_subdiagonal (const gsl_matrix_short * m,
                                    const size_t k);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_superdiagonal (const gsl_matrix_short * m,
                                      const size_t k);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_subrow (const gsl_matrix_short * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_short_const_view
gsl_matrix_short_const_subcolumn (const gsl_matrix_short * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_short_const_view
gsl_matrix_short_const_view_array (const short * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_short_const_view
gsl_matrix_short_const_view_array_with_tda (const short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector (const gsl_vector_short * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_short_const_view
gsl_matrix_short_const_view_vector_with_tda (const gsl_vector_short * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_short_set_zero (gsl_matrix_short * m);
GSL_FUN void gsl_matrix_short_set_identity (gsl_matrix_short * m);
GSL_FUN void gsl_matrix_short_set_all (gsl_matrix_short * m, short x);

GSL_FUN int gsl_matrix_short_fread (FILE * stream, gsl_matrix_short * m) ;
GSL_FUN int gsl_matrix_short_fwrite (FILE * stream, const gsl_matrix_short * m) ;
GSL_FUN int gsl_matrix_short_fscanf (FILE * stream, gsl_matrix_short * m);
GSL_FUN int gsl_matrix_short_fprintf (FILE * stream, const gsl_matrix_short * m, const char * format);

GSL_FUN int gsl_matrix_short_memcpy(gsl_matrix_short * dest, const gsl_matrix_short * src);
GSL_FUN int gsl_matrix_short_swap(gsl_matrix_short * m1, gsl_matrix_short * m2);

GSL_FUN int gsl_matrix_short_swap_rows(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_short_swap_columns(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_short_swap_rowcol(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_short_transpose (gsl_matrix_short * m);
GSL_FUN int gsl_matrix_short_transpose_memcpy (gsl_matrix_short * dest, const gsl_matrix_short * src);

GSL_FUN short gsl_matrix_short_max (const gsl_matrix_short * m);
GSL_FUN short gsl_matrix_short_min (const gsl_matrix_short * m);
GSL_FUN void gsl_matrix_short_minmax (const gsl_matrix_short * m, short * min_out, short * max_out);

GSL_FUN void gsl_matrix_short_max_index (const gsl_matrix_short * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_short_min_index (const gsl_matrix_short * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_short_minmax_index (const gsl_matrix_short * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_short_isnull (const gsl_matrix_short * m);
GSL_FUN int gsl_matrix_short_ispos (const gsl_matrix_short * m);
GSL_FUN int gsl_matrix_short_isneg (const gsl_matrix_short * m);
GSL_FUN int gsl_matrix_short_isnonneg (const gsl_matrix_short * m);

GSL_FUN int gsl_matrix_short_add (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_FUN int gsl_matrix_short_sub (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_FUN int gsl_matrix_short_mul_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_FUN int gsl_matrix_short_div_elements (gsl_matrix_short * a, const gsl_matrix_short * b);
GSL_FUN int gsl_matrix_short_scale (gsl_matrix_short * a, const double x);
GSL_FUN int gsl_matrix_short_add_constant (gsl_matrix_short * a, const double x);
GSL_FUN int gsl_matrix_short_add_diagonal (gsl_matrix_short * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_short_get_row(gsl_vector_short * v, const gsl_matrix_short * m, const size_t i);
GSL_FUN int gsl_matrix_short_get_col(gsl_vector_short * v, const gsl_matrix_short * m, const size_t j);
GSL_FUN int gsl_matrix_short_set_row(gsl_matrix_short * m, const size_t i, const gsl_vector_short * v);
GSL_FUN int gsl_matrix_short_set_col(gsl_matrix_short * m, const size_t j, const gsl_vector_short * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL short   gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x);
GSL_FUN INLINE_DECL short * gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const short * gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
short
gsl_matrix_short_get(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_short_set(gsl_matrix_short * m, const size_t i, const size_t j, const short x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
short *
gsl_matrix_short_ptr(gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (short *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const short *
gsl_matrix_short_const_ptr(const gsl_matrix_short * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_SHORT_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/**
 *
 * @file example_dpotrf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_dmain.h"

static int check_factorization(int, int, double*, double*, int);

static void GENMAT_SYM_FULL(int m, int n, double *A)
{
	srand48(time(NULL));

	int j;
	for (j = 0; j < m; ++j ) {
		int i;
		for( i = j; i < n; ++i ) {
			double dran = drand48();
			A[j*m+i] = A[i*m+j] = dran;
		}
  	}
	for(j = 0; j < m; ++j)
		A[j*m+j] += 10 * m;
}

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

	int M	  = 1000;
    int N     = 1000;
    int LDA   = 1000;
    int info_factorization;

    double *A1   = (double *)malloc(LDA*N*sizeof(double));
    double *A2   = (double *)malloc(LDA*N*sizeof(double));
	#pragma omp register ([LDA*N]A2)

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Initialize A1 and A2 for Symmetric Positive Matrix */
	GENMAT_SYM_FULL(LDA, N, A1);

	int i;
	for(i = 0; i < N*LDA; ++i){
		A2[i] = A1[i];
	}
    /* Plasma routines */
	PLASMA_desc *T;
	PLASMA_Alloc_Workspace_dgels(M, N, &T);
    PLASMA_dgeqrf(M, N, A2, LDA, T);

    /* Check the factorization */
    info_factorization = check_factorization( M, N, A1, A2, LDA);

    if ( info_factorization != 0 )
       printf("-- Error in DPOTRF example ! \n");
    else
       printf("-- Run of DPOTRF example successful ! \n");

    free(A1); free(A2);

    return 0;
}

static int check_factorization(int M, int N, double *A1, double *A2, int LDA)
{
    double Anorm, Rnorm;
	double Anorm1, Rnorm1;
    double alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

//    double *Residual = (double *)malloc(M*N*sizeof(double));
//    double *L1       = (double *)malloc(M*N*sizeof(double));
//    double *L2       = (double *)malloc(M*N*sizeof(double));
    double *work              = (double *)malloc(N*sizeof(double));

	LAPACKE_dgeqrf(LAPACK_COL_MAJOR, M, N, A1, M, work);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A1, M, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A2, LDA, work);

	Rnorm1 = PLASMA_dlange(PlasmaInfNorm, M, N, A1, LDA);
	Anorm1 = PLASMA_dlange(PlasmaInfNorm, M, N, A2, LDA);

	printf("|Rnorm-Rnorm1|: %e, |Anorm-Anorm1|: %e\n", fabs(Rnorm-Rnorm1), fabs(Anorm-Anorm1));
    printf("============\n");
    printf("Checking the QR Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",fabs(Rnorm-Anorm)/(Anorm));

    if ( isnan(fabs(Rnorm-Anorm)/(Anorm)) || (fabs(Rnorm-Anorm)/(Anorm) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

//    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}


text:
/*
 *  Player - One Hell of a Robot Server
 *  Copyright (C) 2000  Brian Gerkey   &  Kasper Stoy
 *                      gerkey@usc.edu    kaspers@robotics.usc.edu
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library 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
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */
/**************************************************************************
 * Desc: Useful pdf functions
 * Author: Andrew Howard
 * Date: 10 Dec 2002
 * CVS: $Id: pf_pdf.c 6348 2008-04-17 02:53:17Z gerkey $
 *************************************************************************/

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
//#include <gsl/gsl_rng.h>
//#include <gsl/gsl_randist.h>

#include "amcl/pf/pf_pdf.h"

// Random number generator seed value
static unsigned int pf_pdf_seed;

/**************************************************************************
 * Gaussian
 *************************************************************************/

// Create a gaussian pdf
pf_pdf_gaussian_t *pf_pdf_gaussian_alloc(pf_vector_t x, pf_matrix_t cx)
{
  pf_matrix_t cd;
  pf_pdf_gaussian_t *pdf;

  pdf = calloc(1, sizeof(pf_pdf_gaussian_t));

  pdf->x = x;
  pdf->cx = cx;
  //pdf->cxi = pf_matrix_inverse(cx, &pdf->cxdet);

  // Decompose the convariance matrix into a rotation
  // matrix and a diagonal matrix.
  pf_matrix_unitary(&pdf->cr, &cd, pdf->cx);
  pdf->cd.v[0] = sqrt(cd.m[0][0]);
  pdf->cd.v[1] = sqrt(cd.m[1][1]);
  pdf->cd.v[2] = sqrt(cd.m[2][2]);

  // Initialize the random number generator
  //pdf->rng = gsl_rng_alloc(gsl_rng_taus);
  //gsl_rng_set(pdf->rng, ++pf_pdf_seed);
  srand48(++pf_pdf_seed);

  return pdf;
}

// Destroy the pdf
void pf_pdf_gaussian_free(pf_pdf_gaussian_t *pdf)
{
  //gsl_rng_free(pdf->rng);
  free(pdf);
  return;
}

/*
// Compute the value of the pdf at some point [x].
double pf_pdf_gaussian_value(pf_pdf_gaussian_t *pdf, pf_vector_t x)
{
  int i, j;
  pf_vector_t z;
  double zz, p;

  z = pf_vector_sub(x, pdf->x);

  zz = 0;
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      zz += z.v[i] * pdf->cxi.m[i][j] * z.v[j];

  p =  1 / (2 * M_PI * pdf->cxdet) * exp(-zz / 2);

  return p;
}
*/

// Generate a sample from the pdf.
pf_vector_t pf_pdf_gaussian_sample(pf_pdf_gaussian_t *pdf)
{
  int i, j;
  pf_vector_t r;
  pf_vector_t x;

  // Generate a random vector
  for (i = 0; i < 3; i++)
  {
    //r.v[i] = gsl_ran_gaussian(pdf->rng, pdf->cd.v[i]);
    r.v[i] = pf_ran_gaussian(pdf->cd.v[i]);
  }

  for (i = 0; i < 3; i++)
  {
    x.v[i] = pdf->x.v[i];
    for (j = 0; j < 3; j++)
      x.v[i] += pdf->cr.m[i][j] * r.v[j];
  }

  return x;
}

// Draw randomly from a zero-mean Gaussian distribution, with standard
// deviation sigma.
// We use the polar form of the Box-Muller transformation, explained here:
//   http://www.taygeta.com/random/gaussian.html
double pf_ran_gaussian(double sigma)
{
  double x1, x2, w, r;

  do
  {
    do { r = drand48(); } while (r==0.0);
    x1 = 2.0 * r - 1.0;
    do { r = drand48(); } while (r==0.0);
    x2 = 2.0 * r - 1.0;
    w = x1*x1 + x2*x2;
  } while(w > 1.0 || w==0.0);

  return(sigma * x2 * sqrt(-2.0*log(w)/w));
}


text:
#ifndef AMICI_MODEL_DIMENSIONS_H
#define AMICI_MODEL_DIMENSIONS_H

#include <gsl/gsl-lite.hpp>
#include <vector>

namespace amici {

/**
 * @brief Container for model dimensions.
 *
 * Holds number of states, observables, etc.
 */
struct ModelDimensions {
    /** Default ctor */
    ModelDimensions() = default;

    /**
     * @brief Constructor with model dimensions
     * @param nx_rdata Number of state variables
     * @param nxtrue_rdata Number of state variables of the non-augmented model
     * @param nx_solver Number of state variables with conservation laws applied
     * @param nxtrue_solver Number of state variables of the non-augmented model
     * with conservation laws applied
     * @param nx_solver_reinit Number of state variables with conservation laws
     * subject to reinitialization
     * @param np Number of parameters
     * @param nk Number of constants
     * @param ny Number of observables
     * @param nytrue Number of observables of the non-augmented model
     * @param nz Number of event observables
     * @param nztrue Number of event observables of the non-augmented model
     * @param ne Number of events
     * @param nJ Number of objective functions
     * @param nw Number of repeating elements
     * @param ndwdx Number of nonzero elements in the `x` derivative of the
     * repeating elements
     * @param ndwdp Number of nonzero elements in the `p` derivative of the
     * repeating elements
     * @param ndwdw Number of nonzero elements in the `w` derivative of the
     * repeating elements
     * @param ndxdotdw Number of nonzero elements in the \f$ w\f$ derivative of \f$ xdot\f$
     * @param ndJydy Number of nonzero elements in the \f$ y\f$ derivative of \f$ dJy\f$ (shape `nytrue`)
     * @param nnz Number of nonzero elements in Jacobian
     * @param ubw Upper matrix bandwidth in the Jacobian
     * @param lbw Lower matrix bandwidth in the Jacobian
     */
    ModelDimensions(
            const int nx_rdata, const int nxtrue_rdata, const int nx_solver,
            const int nxtrue_solver, const int nx_solver_reinit, const int np,
            const int nk, const int ny,
            const int nytrue, const int nz, const int nztrue, const int ne,
            const int nJ, const int nw, const int ndwdx, const int ndwdp,
            const int ndwdw, const int ndxdotdw, std::vector<int> ndJydy,
            const int nnz, const int ubw, const int lbw)
        : nx_rdata(nx_rdata), nxtrue_rdata(nxtrue_rdata), nx_solver(nx_solver),
          nxtrue_solver(nxtrue_solver), nx_solver_reinit(nx_solver_reinit),
          np(np), nk(nk),
          ny(ny), nytrue(nytrue), nz(nz), nztrue(nztrue),
          ne(ne), nw(nw), ndwdx(ndwdx), ndwdp(ndwdp), ndwdw(ndwdw),
          ndxdotdw(ndxdotdw), ndJydy(std::move(ndJydy)),
          nnz(nnz), nJ(nJ), ubw(ubw), lbw(lbw) {
        Expects(nxtrue_rdata >= 0);
        Expects(nxtrue_rdata <= nx_rdata);
        Expects(nxtrue_solver >= 0);
        Expects(nx_solver <= nx_rdata);
        Expects(nxtrue_solver <= nx_solver);
        Expects(nx_solver_reinit >= 0);
        Expects(nx_solver_reinit <= nx_solver);
        Expects(np >= 0);
        Expects(nk >= 0);
        Expects(nytrue <= ny);
        Expects(nytrue >= 0);
        Expects(nztrue >= 0);
        Expects(nztrue <= nz);
        Expects(ne >= 0);
        Expects(nw >= 0);
        Expects(ndwdx >= 0);
        Expects(ndwdx <= nw * nx_solver);
        Expects(ndwdp >= 0);
        Expects(ndwdp <= nw * np);
        Expects(ndwdw >= 0);
        Expects(ndwdw <= nw * nw);
        Expects(ndxdotdw >= 0);
        Expects(nnz >= 0);
        Expects(nJ >= 0);
        Expects(ubw >= 0);
        Expects(lbw >= 0);
    }

    /** Number of states */
    int nx_rdata{0};

    /** Number of states in the unaugmented system */
    int nxtrue_rdata{0};

    /** Number of states with conservation laws applied */
    int nx_solver{0};

    /**
     * Number of states in the unaugmented system with conservation laws
     * applied
     */
    int nxtrue_solver{0};

    /** Number of solver states subject to reinitialization */
    int nx_solver_reinit{0};

    /** Number of parameters */
    int np{0};

    /** Number of constants */
    int nk{0};

    /** Number of observables */
    int ny{0};

    /** Number of observables in the unaugmented system */
    int nytrue{0};

    /** Number of event outputs */
    int nz{0};

    /** Number of event outputs in the unaugmented system */
    int nztrue{0};

    /** Number of events */
    int ne{0};

    /** Number of common expressions */
    int nw{0};

    /**
     * Number of nonzero elements in the `x` derivative of the
     * repeating elements
     */
    int ndwdx {0};

    /**
     * Number of nonzero elements in the `p` derivative of the
     * repeating elements
     */
    int ndwdp {0};

    /**
     * Number of nonzero elements in the `w` derivative of the
     * repeating elements
     */
    int ndwdw {0};

    /** Number of nonzero elements in the \f$ w \f$ derivative of \f$ xdot \f$ */
    int ndxdotdw {0};

    /**
     * Number of nonzero elements in the \f$ y \f$ derivative of
     * \f$ dJy \f$ (dimension `nytrue`)
     */
    std::vector<int> ndJydy;

    /** Number of nonzero entries in Jacobian */
    int nnz{0};

    /** Dimension of the augmented objective function for 2nd order ASA */
    int nJ{0};

    /** Upper bandwidth of the Jacobian */
    int ubw{0};

    /** Lower bandwidth of the Jacobian */
    int lbw{0};
};

} // namespace amici

#endif // AMICI_MODEL_DIMENSIONS_H

max_stars_repo_name: AMICI-dev/AMICI

max_issues_repo_name: AMICI-dev/AMICI

max_forks_repo_name: AMICI-dev/AMICI


text:
#ifndef setupdm_h
#define setupdm_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscfe.h>
#include "../include/structs.h"

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

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

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

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

#endif // setupdm_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#pragma once

#include "gl_helpers.h"
#include "GL_Loader.h"
#include "MappedGL.h"
#include "Range.h"

#include <handy/Guard.h>

#include <gsl/span>

#include <type_traits>
#include <vector>

namespace ad
{

struct [[nodiscard]] VertexArrayObject : public ResourceGuard<GLuint>
{
    /// \note Deleting a bound VAO reverts the binding to zero
    VertexArrayObject() :
        ResourceGuard<GLuint>{reserve(glGenVertexArrays),
                              [](GLuint aIndex){glDeleteVertexArrays(1, &aIndex);}}
    {}
};

/// \TODO understand when glDisableVertexAttribArray should actually be called
///       (likely not specially before destruction, but more when rendering other objects
///        since it is a current(i.e. global) VAO state)
/// \Note Well note even that: Activated vertex attribute array are per VAO, so changing VAO
//        Already correctly handles that.
struct [[nodiscard]] VertexBufferObject : public ResourceGuard<GLuint>
{
    VertexBufferObject() :
        ResourceGuard<GLuint>{reserve(glGenBuffers),
                              [](GLuint aIndex){glDeleteBuffers(1, &aIndex);}}
    {}
};

/// \brief A VertexArray with a vector of VertexBuffers
struct [[nodiscard]] VertexSpecification
{
    VertexSpecification(VertexArrayObject aVertexArray={},
                        std::vector<VertexBufferObject> aVertexBuffers={}) :
        mVertexArray{std::move(aVertexArray)},
        mVertexBuffers{std::move(aVertexBuffers)}
    {}

    VertexArrayObject mVertexArray;
    std::vector<VertexBufferObject> mVertexBuffers;
};

/// \brief Describes the attribute access from the shader (layout id, value type, normalization)
struct Attribute
{
    enum class Access
    {
        Float,
        Integer,
    };

    constexpr Attribute(GLuint aValue) :
        mIndex(aValue)
    {}

    constexpr Attribute(GLuint aValue, Access aAccess, bool aNormalize=false) :
        mIndex(aValue),
        mTypeInShader(aAccess),
        mNormalize(aNormalize)
    {}

    GLuint mIndex;
    Access mTypeInShader{Access::Float};
    bool mNormalize{false};
};

/// \brief The complete description of an attribute expected by OpenGL
struct AttributeDescription : public Attribute
{
    GLuint mDimension;
    size_t mOffset;
    GLenum mDataType;
};

std::ostream & operator<<(std::ostream &aOut, const AttributeDescription & aDescription);

typedef std::initializer_list<AttributeDescription> AttributeDescriptionList;

template <class T_vertex>
VertexBufferObject loadVertexBuffer(const VertexArrayObject & aVertexArray,
                                    const AttributeDescriptionList & aAttributes,
                                    const gsl::span<T_vertex> aVertices,
                                    GLuint aAttributeDivisor = 0)
{
    glBindVertexArray(aVertexArray);

    if (aAttributeDivisor)
    {
        for(const auto & attribute : aAttributes)
        {
            glVertexAttribDivisor(attribute.mIndex, aAttributeDivisor);
        }
    }

    return makeLoadedVertexBuffer(aAttributes, aVertices);
}

template <class T_vertex>
void appendToVertexSpecification(VertexSpecification & aSpecification,
                                 const AttributeDescriptionList & aAttributes,
                                 const gsl::span<T_vertex> aVertices,
                                 GLuint aAttributeDivisor = 0)
{
    aSpecification.mVertexBuffers.push_back(
            loadVertexBuffer(aSpecification.mVertexArray,
                             aAttributes,
                             std::move(aVertices),
                             aAttributeDivisor));
}

/// \brief Create a VertexBufferObject with provided attributes, load it with data.
///
/// This is the lowest level overload, with explicit attribute description and raw data pointer.
/// Other overloads end-up calling it.
VertexBufferObject makeLoadedVertexBuffer(AttributeDescriptionList aAttributes,
                                          GLsizei aStride,
                                          size_t aSize,
                                          const GLvoid * aData);

template <class T_vertex>
VertexBufferObject makeLoadedVertexBuffer(AttributeDescriptionList aAttributes,
                                          const gsl::span<T_vertex> aVertices)
{
    return makeLoadedVertexBuffer(aAttributes,
                                  sizeof(T_vertex),
                                  aVertices.size_bytes(),
                                  aVertices.data());
}

/***
 * Buffer re-specification
 *
 * see: https://www.khronos.org/opengl/wiki/Buffer_Object_Streaming#Buffer_re-specification
 ***/

inline void respecifyBuffer(const VertexBufferObject & aVBO, const GLvoid * aData, GLsizei aSize)
{
    glBindBuffer(GL_ARRAY_BUFFER, aVBO);

    // Orphan the previous buffer
    glBufferData(GL_ARRAY_BUFFER, aSize, NULL, GL_STATIC_DRAW);

    // Copy value to new buffer
    glBufferSubData(GL_ARRAY_BUFFER, 0, aSize, aData);
}

/// \brief Respecify the buffer with the same size (allowing potential optimizations)
inline void respecifyBuffer(const VertexBufferObject & aVBO, const GLvoid * aData)
{
    GLint size;
    glGetBufferParameteriv(GL_ARRAY_BUFFER, GL_BUFFER_SIZE, &size);

    respecifyBuffer(aVBO, aData, size);
}

template <class T_vertex>
void respecifyBuffer(const VertexBufferObject & aVBO, const gsl::span<T_vertex> aVertices)
{
    respecifyBuffer(aVBO,
                    aVertices.data(),
                    static_cast<GLsizei>(aVertices.size_bytes()));
}

} // namespace ad

max_stars_repo_name: FranzPoize/graphics

max_issues_repo_name: FranzPoize/graphics

max_forks_repo_name: FranzPoize/graphics


text:
#ifndef CTETRA_SUM_H
#define CTETRA_SUM_H

#include <stdlib.h>
#include <stdbool.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#include "input.h"
#include "ecache.h"
#include "evcache.h"
#include "submesh.h"
#include "fermi.h"
#include "weights.h"

double SumEnergy(double *E_Fermi, InputFn Efn, int na, int nb, int nc, int num_bands, double num_electrons, gsl_matrix *R, bool use_cache);

double SumEnergyFixedFermi(double E_Fermi, InputFn Efn, int na, int nb, int nc, int num_bands, gsl_matrix *R, bool use_cache);

double with_cache_SumEnergyFixedFermi(double E_Fermi, EnergyCache *Ecache);

double** partial_num_states(UEInputFn UEfn, int na, int nb, int nc, int num_bands, double num_total_electrons, gsl_matrix *R, double **Es, int num_E, double *E_Fermi, double **num_states_Fermi);

double one_band_n(double E_Fermi, int band_index, EvecCache *evCache);

#endif // CTETRA_SUM_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
/* -*- c++ -*- */
/*
 * Copyright 2008,2012 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio 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, or (at your option)
 * any later version.
 *
 * GNU Radio 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 GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef	INCLUDED_WAVELET_SQUASH_FF_IMPL_H
#define	INCLUDED_WAVELET_SQUASH_FF_IMPL_H

#include <gnuradio/wavelet/api.h>
#include <gnuradio/wavelet/squash_ff.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

namespace gr {
  namespace wavelet {

    class WAVELET_API squash_ff_impl : public squash_ff
    {
      size_t  d_inum;
      size_t  d_onum;
      double *d_igrid;
      double *d_iwork;
      double *d_ogrid;

      gsl_interp_accel *d_accel;
      gsl_spline       *d_spline;

    public:
      squash_ff_impl(const std::vector<float> &igrid,
		     const std::vector<float> &ogrid);

      ~squash_ff_impl();

      int work(int noutput_items,
	       gr_vector_const_void_star &input_items,
	       gr_vector_void_star &output_items);
    };

  } /* namespace wavelet */
} /* namespace gr */

#endif /* INCLUDED_WAVELET_WAVELET_FF_IMPL_H */


text:
/* Author: G. Jungman + modifications from O. Teytaud
 */
#ifndef __GSL_QRNG_H__
#define __GSL_QRNG_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>

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

__BEGIN_DECLS

/* Once again, more inane C-style OOP... kill me now. */

/* Structure describing a type of generator.
 */
typedef struct
{
  const char * name;
  unsigned int max_dimension;
  size_t (*state_size) (unsigned int dimension);
  int (*init_state) (void * state, unsigned int dimension);
  int (*get) (void * state, unsigned int dimension, double x[]);
}
gsl_qrng_type;

/* Structure describing a generator instance of a
 * specified type, with generator-specific state info
 * and dimension-specific info.
 */
typedef struct
{
  const gsl_qrng_type * type;
  unsigned int dimension;
  size_t state_size;
  void * state;
}
gsl_qrng;

/* Supported generator types.
 */
GSL_VAR const gsl_qrng_type * gsl_qrng_niederreiter_2;
GSL_VAR const gsl_qrng_type * gsl_qrng_sobol;
GSL_VAR const gsl_qrng_type * gsl_qrng_halton;
GSL_VAR const gsl_qrng_type * gsl_qrng_reversehalton;

/* Allocate and initialize a generator
 * of the specified type, in the given
 * space dimension.
 */
gsl_qrng * gsl_qrng_alloc (const gsl_qrng_type * T, unsigned int dimension);

/* Copy a generator. */
int gsl_qrng_memcpy (gsl_qrng * dest, const gsl_qrng * src);

/* Clone a generator. */
gsl_qrng * gsl_qrng_clone (const gsl_qrng * q);

/* Free a generator. */
void gsl_qrng_free (gsl_qrng * q);

/* Intialize a generator. */
void gsl_qrng_init (gsl_qrng * q);

/* Get the standardized name of the generator. */
const char * gsl_qrng_name (const gsl_qrng * q);

/* ISN'T THIS CONFUSING FOR PEOPLE?
  WHAT IF SOMEBODY TRIES TO COPY WITH THIS ???
  */
size_t gsl_qrng_size (const gsl_qrng * q);

void * gsl_qrng_state (const gsl_qrng * q);

/* Retrieve next vector in sequence. */
INLINE_DECL int gsl_qrng_get (const gsl_qrng * q, double x[]);

#ifdef HAVE_INLINE
INLINE_FUN int gsl_qrng_get (const gsl_qrng * q, double x[])
{
  return (q->type->get) (q->state, q->dimension, x);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* !__GSL_QRNG_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


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

#ifndef __GSL_INTERP2D_H__
#define __GSL_INTERP2D_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_interp.h>

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

__BEGIN_DECLS

typedef struct {
    const char* name;
    unsigned int min_size;
    void * (*alloc)(size_t xsize, size_t ysize);
    int    (*init)(void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize);
    int    (*eval)(const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z);
    int    (*eval_deriv_x) (const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z_p);
    int    (*eval_deriv_y) (const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z_p);
    int    (*eval_deriv_xx) (const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z_pp);
    int    (*eval_deriv_xy) (const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z_pp);
    int    (*eval_deriv_yy) (const void *, const double xa[], const double ya[], const double za[], size_t xsize, size_t ysize, double x, double y, gsl_interp_accel*, gsl_interp_accel*, double* z_pp);
    void   (*free)(void *);
} gsl_interp2d_type;

typedef struct {
    const gsl_interp2d_type * type; /* interpolation type */
    double xmin;                    /* minimum value of x for which data have been provided */
    double xmax;                    /* maximum value of x for which data have been provided */
    double ymin;                    /* minimum value of y for which data have been provided */
    double ymax;                    /* maximum value of y for which data have been provided */
    size_t xsize;                   /* number of x values provided */
    size_t ysize;                   /* number of y values provided */
    void * state;                   /* internal state object specific to the interpolation type */
} gsl_interp2d;

/* available types */
GSL_VAR const gsl_interp2d_type * gsl_interp2d_bilinear;
GSL_VAR const gsl_interp2d_type * gsl_interp2d_bicubic;

GSL_FUN gsl_interp2d * gsl_interp2d_alloc(const gsl_interp2d_type * T, const size_t xsize,
                                  const size_t ysize);

GSL_FUN const char * gsl_interp2d_name(const gsl_interp2d * interp);
GSL_FUN size_t gsl_interp2d_min_size(const gsl_interp2d * interp);
GSL_FUN size_t gsl_interp2d_type_min_size(const gsl_interp2d_type * T);
GSL_FUN int gsl_interp2d_set(const gsl_interp2d * interp, double zarr[],
                     const size_t i, const size_t j, const double z);
GSL_FUN double gsl_interp2d_get(const gsl_interp2d * interp, const double zarr[],
                        const size_t i, const size_t j);
GSL_FUN size_t gsl_interp2d_idx(const gsl_interp2d * interp,
                        const size_t i, const size_t j);
GSL_FUN int gsl_interp2d_init(gsl_interp2d * interp, const double xa[], const double ya[],
                      const double za[], const size_t xsize, const size_t ysize);
GSL_FUN void gsl_interp2d_free(gsl_interp2d * interp);

GSL_FUN double gsl_interp2d_eval(const gsl_interp2d * interp, const double xarr[],
                         const double yarr[], const double zarr[], const double x,
                         const double y, gsl_interp_accel * xa, gsl_interp_accel * ya);

GSL_FUN double gsl_interp2d_eval_extrap(const gsl_interp2d * interp,
                                const double xarr[], const double yarr[],
                                const double zarr[], const double x,
                                const double y, gsl_interp_accel * xa,
                                gsl_interp_accel * ya);

GSL_FUN int gsl_interp2d_eval_e(const gsl_interp2d * interp, const double xarr[],
                        const double yarr[], const double zarr[],
                        const double x, const double y, gsl_interp_accel* xa,
                        gsl_interp_accel* ya, double * z);

#ifndef GSL_DISABLE_DEPRECATED

GSL_FUN int gsl_interp2d_eval_e_extrap(const gsl_interp2d * interp,
                               const double xarr[],
                               const double yarr[],
                               const double zarr[],
                               const double x,
                               const double y,
                               gsl_interp_accel * xa,
                               gsl_interp_accel * ya,
                               double * z);

#endif /* !GSL_DISABLE_DEPRECATED */

GSL_FUN int gsl_interp2d_eval_extrap_e(const gsl_interp2d * interp,
                               const double xarr[],
                               const double yarr[],
                               const double zarr[],
                               const double x,
                               const double y,
                               gsl_interp_accel * xa,
                               gsl_interp_accel * ya,
                               double * z);

GSL_FUN double gsl_interp2d_eval_deriv_x(const gsl_interp2d * interp, const double xarr[],
                                 const double yarr[], const double zarr[],
                                 const double x, const double y, gsl_interp_accel * xa,
                                 gsl_interp_accel * ya);

GSL_FUN int gsl_interp2d_eval_deriv_x_e(const gsl_interp2d * interp, const double xarr[],
                                const double yarr[], const double zarr[],
                                const double x, const double y,
                                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z);

GSL_FUN double gsl_interp2d_eval_deriv_y(const gsl_interp2d * interp, const double xarr[],
                                 const double yarr[], const double zarr[],
                                 const double x, const double y,
                                 gsl_interp_accel* xa, gsl_interp_accel* ya);

GSL_FUN int gsl_interp2d_eval_deriv_y_e(const gsl_interp2d * interp, const double xarr[],
                                const double yarr[], const double zarr[],
                                const double x, const double y,
                                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z);

GSL_FUN double gsl_interp2d_eval_deriv_xx(const gsl_interp2d * interp, const double xarr[],
                                  const double yarr[], const double zarr[],
                                  const double x, const double y,
                                  gsl_interp_accel * xa, gsl_interp_accel * ya);

GSL_FUN int gsl_interp2d_eval_deriv_xx_e(const gsl_interp2d * interp, const double xarr[],
                                 const double yarr[], const double zarr[],
                                 const double x, const double y,
                                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z);

GSL_FUN double gsl_interp2d_eval_deriv_yy(const gsl_interp2d * interp, const double xarr[],
                                  const double yarr[], const double zarr[],
                                  const double x, const double y,
                                  gsl_interp_accel * xa, gsl_interp_accel * ya);

GSL_FUN int gsl_interp2d_eval_deriv_yy_e(const gsl_interp2d * interp, const double xarr[],
                                 const double yarr[], const double zarr[],
                                 const double x, const double y,
                                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z);

GSL_FUN double gsl_interp2d_eval_deriv_xy(const gsl_interp2d * interp, const double xarr[],
                                  const double yarr[], const double zarr[],
                                  const double x, const double y,
                                  gsl_interp_accel * xa, gsl_interp_accel * ya);

GSL_FUN int gsl_interp2d_eval_deriv_xy_e(const gsl_interp2d * interp, const double xarr[],
                                 const double yarr[], const double zarr[],
                                 const double x, const double y,
                                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z);

__END_DECLS

#endif /* __GSL_INTERP2D_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim