text:
#ifndef __MTPOWERSPECTRUM_ENGINE_H__
#define  __MTPOWERSPECTRUM_ENGINE_H__

#include <memory>
#include <vector>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include "mspass/seismic/TimeSeries.h"
#include "mspass/utility/dmatrix.h"
#include "mspass/seismic/PowerSpectrum.h"

namespace mspass::algorithms::deconvolution{
/*! \brief Multittaper power spectral estimator.

The multitaper method uses averages of spectra windowed by Slepian functions.
This class can be used to compute power spectra.  For efficiency the design
has constructors that build the Slepian functions and cache them in a
private area.  Use of the pply method that returns vector of
power spectral estimates.  Get methods can be called to get frequency
properties of the data returned.  This could have been done with a class, but
that was judged better left to a python wrapper using this class.
*/
class MTPowerSpectrumEngine
{
public:
  MTPowerSpectrumEngine();
  MTPowerSpectrumEngine(const int winsize, const double tbp, const int ntapers);
  MTPowerSpectrumEngine(const MTPowerSpectrumEngine& parent);
  ~MTPowerSpectrumEngine();
  //~MTPowerSpectrumEngine();
  MTPowerSpectrumEngine& operator=(const MTPowerSpectrumEngine& parent);
  /*! \process a TimeSeries.

  This is one of two methods for applying the multiaper algorithm to data.
  This one uses dt and data length to set the Rayleigh bin size (df).   If
  the input data vector length is not the same as the operator length an
  elog complaint is posted to parent.   Short data are processed but should
  be considered suspect unless the sizes differ by only a tiny fraction
  (e.g. and off by one error from rounding).  Long data will be truncated
  on the right (i.e. sample 0 will be the start of the window used).

  \param parent is the data to process
  \return vector containing estimated power spwecrum
  */
  mspass::seismic::PowerSpectrum apply(const mspass::seismic::TimeSeries& d);
  /*! \brief Low level processing of vector of data.

  This is lower level function that processes a raw vector of data.   Since
  it does not know the sample interval it cannot compute the rayleigh bin
  size so if callers need that feature they must do that (simple) calculation
  themselves.   Unlike the TimeSeries method this one will throw an
  exception if the input data size does not match the operator size.

  \param d is the vector of data to process.  d.size() must this->taperlen() value.
  \return vector containing estimated power spectrum (usual convention with
    0 containing 0 frequency value)
  \exception throw a MsPASSError if the size of d does not match operator length
  */
  std::vector<double> apply(const std::vector<double>& d);
  double df() const {return deltaf;};

  std::vector<double> frequencies();
  /*! Retrieve the taper length.*/
  int taper_length() const
  {
    return taperlen;
  };
  /*! Retrieve time-bandwidth product.*/
  double time_bandwidth_product()  const
  {
    return tbp;
  };
  /*! Return number of tapers used by this engine. */
  int number_tapers() const
  {
    return ntapers;
  };
  /*! \brief PUtter equivalent of df.

  The computation of the Rayleigh bin size (dt) is actually quote trivial but
  this convenience functon allows users of the  vector<double> method to
  handle the comutation easily.  It uses the internal oeprator size, however,
  to compute the df size it returns because the operator is dogmatic about
  using that size.  Users wishing to call the frequency method and the apply
  method on raw vector data need to call this function before calling frequency.

  \param dt is the data sample interval (time domain)

  \return computed df
  */
  double set_df(double dt)
  {
    deltaf=1.0/(dt*static_cast<double>(taperlen));
    return deltaf;
  };
private:
  int taperlen;
  int ntapers;
  double tbp;
  mspass::utility::dmatrix tapers;
  /* Frequency bin interval of last data processed.*/
  double deltaf;
  gsl_fft_complex_wavetable *wavetable;
  gsl_fft_complex_workspace *workspace;
};
} //namespace ed
#endif

max_stars_repo_name: Yangzhengtang/mspass

max_issues_repo_name: Yangzhengtang/mspass

max_forks_repo_name: Yangzhengtang/mspass


text:
/* ******************************************************************************
 *
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://www.apache.org/licenses/LICENSE-2.0.
 *
 *  See the NOTICE file distributed with this work for additional
 *  information regarding copyright ownership.
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
 * License for the specific language governing permissions and limitations
 * under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
// Created by agibsonccc on 2/6/16.
//

#ifndef NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H
#define NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H
#include <cblas.h>

/*
enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102 };
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113,
    AtlasConj=114};
enum CBLAS_UPLO  {CblasUpper=121, CblasLower=122};
enum CBLAS_DIAG  {CblasNonUnit=131, CblasUnit=132};
enum CBLAS_SIDE  {CblasLeft=141, CblasRight=142};
*/

#ifdef __cplusplus
extern "C" {
#endif
/**
 * Converts a character
 * to its proper enum
 * for row (c) or column (f) ordering
 * default is row major
 */
CBLAS_ORDER convertOrder(int from);
/**
 * Converts a character to its proper enum
 * t -> transpose
 * n -> no transpose
 * c -> conj
 */
CBLAS_TRANSPOSE convertTranspose(int from);
/**
 * Upper or lower
 * U/u -> upper
 * L/l -> lower
 *
 * Default is upper
 */
CBLAS_UPLO convertUplo(int from);

/**
 * For diagonals:
 * u/U -> unit
 * n/N -> non unit
 *
 * Default: unit
 */
CBLAS_DIAG convertDiag(int from);
/**
 * Side of a matrix, left or right
 * l /L -> left
 * r/R -> right
 * default: left
 */
CBLAS_SIDE convertSide(int from);

#ifdef __cplusplus
}
#endif

#endif  // NATIVEOPERATIONS_CBLAS_ENUM_CONVERSION_H_H

ext: h


text:
#include "stdlib.h"
#include "stdio.h"
#include "/home/lillian/work/install_fpdebug/valgrind-3.7.0/fpdebug/fpdebug.h"
#include <gsl/gsl_statistics.h>
int main(int argc, const char * argv[]) {
unsigned long int hexdouble;
int a;
a = atoi(argv[1]);
int argv_cnt = 2;
double b[5];
for(int i = 0; i < a; i++) {
sscanf(argv[argv_cnt++], "%lX", &hexdouble);
b[i] = *(double*)(&hexdouble);
}
unsigned long c = 1;
double d[5];
for(int i = 0; i < a; i++) {
sscanf(argv[argv_cnt++], "%lX", &hexdouble);
d[i] = *(double*)(&hexdouble);
}
unsigned long e = 1;
double g;
sscanf(argv[argv_cnt++], "%lX", &hexdouble);
g = *(double*)(&hexdouble);
double h;
sscanf(argv[argv_cnt++], "%lX", &hexdouble);
h = *(double*)(&hexdouble);
double result = gsl_stats_wkurtosis_m_sd(b, c, d, e, a, g, h);
//printf("%.15f\n", result);
VALGRIND_PRINT_VALUES("result", 1, &result);
return 0;
}

ext: c

max_stars_repo_name: floatfeather/FpGenetic

max_issues_repo_name: floatfeather/FpGenetic

max_forks_repo_name: floatfeather/FpGenetic


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#include <cairo.h>
#include <gtk/gtk.h>
#include <gsl/gsl_multifit.h>
#include <kplot.h>

#include "extern.h"

void
draw(GtkWidget *w, cairo_t *cr, struct curwin *cur)
{
	double		 x, y;

	cur->redraw = 0;

	/* White-out the view. */
	x = gtk_widget_get_allocated_width(w);
	y = gtk_widget_get_allocated_height(w);
	cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
	cairo_rectangle(cr, 0.0, 0.0, x, y);
	cairo_fill(cr);

	/* Draw our plot. */
	kplot_draw(cur->views[cur->view], x, y, cr);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
/*
 * SPDX-License-Identifier:  BSD-3-Clause
 *
 * sph_linked_list.c :
 *     Header containing the declarions for several
 *     cell linked list operations, including hash
 *     calculations, setup hash tables and neighbour finding.
 *
 * (C) Copyright 2021 José Hugo Elsas
 * Author: José Hugo Elsas <jhelsas@gmail.com>
 *
 */

#include <math.h>
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <stdbool.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "MZC3D64.h"
#include "MZC2D64.h"
#include "sph_data_types.h"
#include "sph_linked_list.h"

#define dbg false

int safe_free_box(linkedListBox *box){
	kh_destroy(0, box->hbegin);
	kh_destroy(1, box->hend);
	free(box);

	return 0;
}

int compare_int64_t(const void *p,const void *q){
	int64_t *data1,*data2;
	data1 = (int64_t*)p;
	data2 = (int64_t*)q;

	if(data1[0] < data2[0])         // data[0] is the hash value,
		return -1;                    // data[1] is the corresponding
	else if(data1[0] == data2[0])   // in the unsorted array
		return 0;
	else
		return 1;
}

#define safe_check_alloc(ptr,N,dtype) {                                          \
                                       (ptr) = (dtype*)malloc((N)*sizeof(dtype));\
                                       if((ptr)==NULL){                          \
                                         success=0;                              \
                                         goto finishlabel;                       \
                                       }                                         \
                                      }

#define safe_check__aligned_alloc(ptr,alignment,N,dtype) {                                        \
                                       (ptr) = (dtype*)aligned_alloc(alignment,(N)*sizeof(dtype));\
                                       if((ptr)==NULL){                                           \
                                         success=0;                                               \
                                         goto finishlabel;                                        \
                                       }                                                          \
                                      }

#define safe_free(ptr) {               \
                        if(ptr != NULL)\
                          free(ptr);   \
                       }

/*
 *  Function SPHparticle_SoA_malloc:
 *    Easily allocate the SPHparticle array
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particle struct reference
 */
int SPHparticle_SoA_malloc(int N,SPHparticle **lsph){
  int success=1;
  //const int alignment = 32;
  (*lsph) = (SPHparticle*)malloc(1*sizeof(SPHparticle));
  if(lsph==NULL){
    success = 0;
    goto finishlabel;
  }
  (*lsph)->x  = NULL; (*lsph)->y  = NULL; (*lsph)->z  = NULL;
  (*lsph)->ux = NULL; (*lsph)->uy = NULL; (*lsph)->uz = NULL;
  (*lsph)->Fx = NULL; (*lsph)->Fy = NULL; (*lsph)->Fz = NULL;
  (*lsph)->nu = NULL; (*lsph)->rho= NULL;
  (*lsph)->id = NULL; (*lsph)->hash= NULL;

  safe_check_alloc((*lsph)->x   , N ,double);
  safe_check_alloc((*lsph)->y   , N ,double);
  safe_check_alloc((*lsph)->z   , N ,double);
  safe_check_alloc((*lsph)->ux  , N ,double);
  safe_check_alloc((*lsph)->uy  , N ,double);
  safe_check_alloc((*lsph)->uz  , N ,double);
  safe_check_alloc((*lsph)->Fx  , N ,double);
  safe_check_alloc((*lsph)->Fy  , N ,double);
  safe_check_alloc((*lsph)->Fz  , N ,double);
  safe_check_alloc((*lsph)->nu  , N ,double);
  safe_check_alloc((*lsph)->rho , N ,double);
  safe_check_alloc((*lsph)->id  , N ,int64_t);
  safe_check_alloc((*lsph)->hash,2*N,int64_t);

  // SPHparticle_SoA_malloc((*lsph)->x   , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->y   , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->z   , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->ux  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->uy  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->uz  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->Fx  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->Fy  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->Fz  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->nu  , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->rho , 32, N ,double);
  // SPHparticle_SoA_malloc((*lsph)->id  , 32, N ,int64_t);
  // SPHparticle_SoA_malloc((*lsph)->hash, 32,2*N,int64_t);

finishlabel:

  if(success)
    return 0;
  else{
    if(*lsph==NULL)
      return 1;

    safe_free((*lsph)->x);  safe_free((*lsph)->y);  safe_free((*lsph)->z);
    safe_free((*lsph)->ux); safe_free((*lsph)->uy); safe_free((*lsph)->uz);
    safe_free((*lsph)->Fx); safe_free((*lsph)->Fy); safe_free((*lsph)->Fz);
    safe_free((*lsph)->nu); safe_free((*lsph)->rho);
    safe_free((*lsph)->id); safe_free((*lsph)->hash);

    return 1;
  }
}

/*
 *  Function SPHparticleSOA_safe_free:
 *    Easily free the SPHparticle array
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particle struct reference
 */
int SPHparticleSOA_safe_free(int N,SPHparticle **lsph){

  if(*lsph==NULL)
    return 1;

  safe_free((*lsph)->x);  safe_free((*lsph)->y);  safe_free((*lsph)->z);
  safe_free((*lsph)->ux); safe_free((*lsph)->uy); safe_free((*lsph)->uz);
  safe_free((*lsph)->Fx); safe_free((*lsph)->Fy); safe_free((*lsph)->Fz);
  safe_free((*lsph)->nu); safe_free((*lsph)->rho);
  safe_free((*lsph)->id); safe_free((*lsph)->hash);

  free((*lsph));

  return 0;
}

/*
 *  Function gen_unif_rdn_pos:
 *    Generate particles positions uniformely in a [0,1]^3 box
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particles array
 *     seed         <int>     : seed for the PRNG
 */
int gen_unif_rdn_pos(int64_t N, int seed, SPHparticle *lsph){

	const gsl_rng_type *T=NULL;
	gsl_rng *r=NULL;

	if(lsph==NULL)
		return 1;

	gsl_rng_env_setup();                      // Initialize environment for PRNG generator

	T = gsl_rng_default;                      // Set generator type as default type
	r = gsl_rng_alloc(T);                     // allocate state for it
	gsl_rng_set(r,seed);                      // set seed accordingly

	for(int64_t i=0;i<N;i+=1){
		lsph->x[i] = gsl_rng_uniform(r);        // Generate a uniform value for X component between 0 and 1
		lsph->y[i] = gsl_rng_uniform(r);        // Generate a uniform value for Y component between 0 and 1
		lsph->z[i] = gsl_rng_uniform(r);        // Generate a uniform value for Z component between 0 and 1

		lsph->ux[i] = 0.0; lsph->Fx[i] = 0.0;   // set velocity and force components as zero
		lsph->uy[i] = 0.0; lsph->Fy[i] = 0.0;   // as they will not be needed
    lsph->uz[i] = 0.0; lsph->Fz[i] = 0.0;

		lsph->nu[i]   = 1.0/N;                  // set each particle mass as to the total mass be 1
		lsph->rho[i]  = 0.0;                    // initialize density to zero
		lsph->id[i]   = (int64_t) i;            // set the particle's id
		lsph->hash[2*i+0] = (int64_t) 0;        // initialize particle's hash value as zero
		lsph->hash[2*i+1] = (int64_t) i;        // set particle's index position by the position in current array
	}

	gsl_rng_free(r);                          // release the PRNG

	return 0;
}

/*
 *  Function gen_unif_rdn_pos_box:
 *    Generate particles positions uniformely in a
 *    [Xmin,Xmax]x[Ymin,Ymax]x[Zmin,Zmax] box
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particles array
 *     seed         <int>     : seed for the PRNG
 *      box  <linkedListBox*> : Cell Linked List box
 */
int gen_unif_rdn_pos_box(int64_t N, int seed, linkedListBox *box,SPHparticle *lsph){

  const gsl_rng_type *T=NULL;
  gsl_rng *r=NULL;

  if(lsph==NULL)
    return 1;

	gsl_rng_env_setup();                      // Initialize environment for PRNG generator

	T = gsl_rng_default;                      // Set generator type as default type
	r = gsl_rng_alloc(T);                     // allocate state for it
	gsl_rng_set(r,seed);                      // set seed accordingly

  for(int64_t i=0;i<N;i+=1){
    lsph->x[i] = gsl_rng_uniform(r)*(box->Xmax-box->Xmin) // Generate a uniform value for X
                     + box->Xmin;                         // component between Xmin and Xmax
    lsph->y[i] = gsl_rng_uniform(r)*(box->Ymax-box->Ymin) // Generate a uniform value for Y
                     + box->Ymin;                         // component between Ymin and Ymax
    lsph->z[i] = gsl_rng_uniform(r)*(box->Zmax-box->Zmin) // Generate a uniform value for Z
                     + box->Zmin;                         // component between Zmin and Zmax

		lsph->ux[i] = 0.0; lsph->Fx[i] = 0.0;   // set velocity and force components as zero
		lsph->uy[i] = 0.0; lsph->Fy[i] = 0.0;   // as they will not be needed
    lsph->uz[i] = 0.0; lsph->Fz[i] = 0.0;

		lsph->nu[i]   = 1.0/N;                  // set each particle mass as to the total mass be 1
		lsph->rho[i]  = 0.0;                    // initialize density to zero
		lsph->id[i]   = (int64_t) i;            // set the particle's id
		lsph->hash[2*i+0] = (int64_t) 0;        // initialize particle's hash value as zero
		lsph->hash[2*i+1] = (int64_t) i;        // set particle's index position by the position in current array
  }

  gsl_rng_free(r);

  return 0;
}

/*
 *  Function compute_hash_MC3D:
 *    Compute the Morton Z hashes for all particles based on
 *    the parameters given by the Cell Linked List box
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particles array
 *     box  <linkedListBox*>  : Box of cell linked lists
 */
int compute_hash_MC3D(int64_t N, SPHparticle *lsph, linkedListBox *box){

	if(lsph==NULL)
		return 1;

	if(box==NULL)
		return 1;

	const double etax = box->Nx/(box->Xmax - box->Xmin);
	const double etay = box->Ny/(box->Ymax - box->Ymin);
	const double etaz = box->Nz/(box->Zmax - box->Zmin);

	for(int64_t i=0;i<N;i+=1){
		uint32_t kx,ky,kz;
		kx = (uint32_t)(etax*(lsph->x[i] - box->Xmin));      // Compute cell index in the X direction
		ky = (uint32_t)(etay*(lsph->y[i] - box->Ymin));      // Compute cell index in the Y direction
		kz = (uint32_t)(etaz*(lsph->z[i] - box->Zmin));      // Compute cell index in the Z direction

		if((kx<0)||(ky<0)||(kz<0))                           // The can't be negative indexes
			return 1;
		else if((kx>=box->Nx)||(ky>=box->Nx)||(kz>=box->Nx)) // Nor indexes greater than the upper bounds
			return 1;
		else{
			lsph->hash[2*i+0] = ullMC3Dencode(kx,ky,kz);         // If ok, compute the corresponding Morton Z hash
      lsph->hash[2*i+1] = i;
		}
	}

	return 0;
}

/*
 * To swap two arrays, copy data from the original array into
 * a swap buffer in the right order, and then overrite the original array.
 */
#define swap_loop(N,lsph,temp_swap,member,type) for(int64_t i=0;i<(N);i+=1)                            \
																	 	              (temp_swap)[i] = (lsph)->member[(lsph)->hash[2*i+1]];\
																		            memcpy((lsph)->member,temp_swap,(N)*sizeof(type))

/*
 *  Function reorder_lsph_SoA:
 *    reorder the array according to the indexes reordered according to
 *    their hashes.
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particles array
 *     swap_arr  <void*>      : buffer for reordering data in the array
 */
int reorder_lsph_SoA(int N, SPHparticle *lsph, void *swap_arr){

	int64_t *int64_temp_swap = (int64_t *)swap_arr;
	swap_loop(N,lsph,int64_temp_swap,id ,int64_t);
	double *double_temp_swap = (double *)swap_arr;
	swap_loop(N,lsph,double_temp_swap,nu ,double);
	swap_loop(N,lsph,double_temp_swap,rho,double);
	swap_loop(N,lsph,double_temp_swap,x  ,double);
	swap_loop(N,lsph,double_temp_swap,y  ,double);
	swap_loop(N,lsph,double_temp_swap,z  ,double);
	swap_loop(N,lsph,double_temp_swap,ux ,double);
	swap_loop(N,lsph,double_temp_swap,uy ,double);
	swap_loop(N,lsph,double_temp_swap,uz ,double);
	swap_loop(N,lsph,double_temp_swap,Fx ,double);
	swap_loop(N,lsph,double_temp_swap,Fy ,double);
	swap_loop(N,lsph,double_temp_swap,Fz ,double);

	return 0;
}

/*
 *  Function setup_interval_hashtables:
 *    Setup the information where each cell begins and ends in the SPH
 *    Array in a pair of hash tables for quick consultation.
 *    This effectively completes the creation of the implicit
 *    cell linked list structure initiated by sorting the array
 *    according to hashes.
 *
 *    Arguments:
 *       N          <int>     : Number of Particles
 *     lsph <SPHparticle**>   : SoA Particles array
 *     box  <linkedListBox*>  : Box of cell linked lists
 */
int setup_interval_hashtables(int64_t N,SPHparticle *lsph,linkedListBox *box){

	int ret;
	int64_t hash0 = lsph->hash[2*0];                             // Store the first hash in hash0
	khiter_t kbegin,kend;

	if(lsph==NULL)
		return 1;

	if(box==NULL)
		return 1;

	kbegin = kh_put(0, box->hbegin, lsph->hash[2*0], &ret);      // Insert start hash as first value
	kh_value(box->hbegin, kbegin) = (int64_t)0;                  // Set the start index of first cell as 0
	for(int64_t i=0;i<N;i+=1){
		lsph->hash[i] = lsph->hash[2*i];                           // Clean the first half of the array only to contain hashes
		if(lsph->hash[i] == hash0)                                 // If the the particle hash is the same as the previous
			continue;                                                // skip because they are in the same cell
		hash0 = lsph->hash[i];                                     // otherwise, store the new hash in hash0

		kend   = kh_put(1, box->hend  , lsph->hash[i-1], &ret);    // If so, insert the previous hash as an cell end
		kh_value(box->hend  , kend)   = i;                         // and set this index as the upper index of that cell

		kbegin = kh_put(0, box->hbegin, lsph->hash[i  ], &ret);    // and then insert this new hash as the begining
		kh_value(box->hbegin, kbegin) = i;                         // of a new cell in with the lower index as this index
	}
	kend   = kh_put(1, box->hend  , lsph->hash[2*(N-1)], &ret);  // finally set the last hash for the end of the last cell,
	kh_value(box->hend  , kend)   = N;                           // and the number of particles as the index of the last cell

	return 0;
}

/*
 *  Function neighbour_hash_3d:
 *    Find the hashes of all valid cells that neighbor cell with
 *    int hash.
 *
 *    Arguments:
 *       hash    <int64_t>    : Hash of the central cell
 *       nblist <int64_t*>    : Array for hashes of the neighbor cells
 *       width     <int>      : Width of the neighborhood
 *     box  <linkedListBox*>  : Box of cell linked lists
 */
int neighbour_hash_3d(int64_t hash,int64_t *nblist,int width, linkedListBox *box){
	int idx=0,kx=0,ky=0,kz=0;

	kx = ullMC3DdecodeX(hash);   // To find a given cell's neighbors, start by first
	ky = ullMC3DdecodeY(hash);   // recovering the integer index of that cell in
	kz = ullMC3DdecodeZ(hash);   // the overall box domain

	for(int ix=-width;ix<=width;ix+=1){            // iterate over all cells surrounding the original cell
		for(int iy=-width;iy<=width;iy+=1){          // as long as they are within width of the original cell
			for(int iz=-width;iz<=width;iz+=1){
				if((kx+ix<0)||(ky+iy<0)||(kz+iz<0))                              // if the cell is off bounds for having lower index then 0
					nblist[idx++] = -1;                                            // annotate the cell as invalid
				else if( (kx+ix>=box->Nx)||(ky+iy>=box->Ny)||(kz+iz>=box->Nz) )  // or if it is off bounds for having higher index then
					nblist[idx++] = -1;                                            // the maximum, also annotate as invalid
				else if( kh_get(0, box->hbegin, ullMC3Dencode(kx+ix,ky+iy,kz+iz)) == kh_end(box->hbegin) ) // Also, if no corresponding hash
					nblist[idx++] = -1;                                            // is not found in the hash table, also annotate as invalid
				else                                                             // at last, if it has no problems
					nblist[idx++] = ullMC3Dencode(kx+ix,ky+iy,kz+iz);              // annotate the hash of the corresponding cell
			}
		}
	}

	return 0;
}

/*
 *  Function count_box_pairs:
 *    Count the number of valid cell pairs
 *
 *    Arguments:
 *     box  <linkedListBox*>  : Box of cell linked lists
 */
int count_box_pairs(linkedListBox *box){
  int64_t pair_count = 0;                                                     // initialize the number of pairs as zero

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){ // iterate over the cells
    int64_t node_hash=-1;
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

    if (kh_exist(box->hbegin, kbegin)){                                       // check if the cell is valid
      node_hash = kh_key(box->hbegin, kbegin);                                // get the hash of that cell

      neighbour_hash_3d(node_hash,nblist,box->width,box);                     // fetch the list of neighbors
      for(int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){ // and iterate over those neighbors
        if(nblist[j]>=0){                                                     // if the neighbor is valid
          pair_count += 1;                                                    // count it as a valid for a pair
        }
      }
    }
  }

  return pair_count;
}

/*
 *  Function setup_box_pairs:
 *    Pre-compute indexes for all valid cell pairs.
 *
 *  Arguments:
 *     box  <linkedListBox*>  : Box of cell linked lists
 *  Returns:
 *     node_begin <int64_t*>  : Array for receiver cell begin indexes
 *     node_end   <int64_t*>  : Array for receiver cell end indexes
 *     nb_begin   <int64_t*>  : Array for sender cell begin indexes
 *     nb_end     <int64_t*>  : Array for sender cell end indexes
 */
int setup_box_pairs(linkedListBox *box,
                    int64_t *node_begin,int64_t *node_end,
                    int64_t *nb_begin,int64_t *nb_end)
{
  int64_t pair_count = 0,particle_pair_count = 0;

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){    // iterate over the cells
    int64_t node_hash=-1;
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

    if (kh_exist(box->hbegin, kbegin)){                                                       // check if the cell is valid
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));                     // get the iterator for the end index

      node_hash = kh_key(box->hbegin, kbegin);                                                // and also get the corresponding hash

      neighbour_hash_3d(node_hash,nblist,box->width,box);                                     // fetch the list of neighbors
      for(int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){                 // and iterate over those neighbors
        if(nblist[j]>=0){                                                                     // check if the neighbor is valid
          node_begin[pair_count] = kh_value(box->hbegin, kbegin);                             // get the cell's begin index
          node_end[pair_count]   = kh_value(box->hend, kend);                                 // get the cell's end index
          nb_begin[pair_count]   = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) ); // get the neighbor's begin index
          nb_end[pair_count]     = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) ); // get the neighbor's end index

          particle_pair_count += (nb_end[pair_count]-nb_begin[pair_count])*                   // just for bookeeping/profiling,
                                       (node_end[pair_count]-node_begin[pair_count]);         // annotate the number of particle pairs
          pair_count += 1;                                                                    // that are part of the computation
        }
      }
    }
  }

  if(dbg)
    printf("particle_pair_count = %ld\n",particle_pair_count);

  return pair_count;
}

/*
 *  Function setup_unique_box_pairs:
 *    Pre-compute indexes for all valid unique cell pairs.
 *    by unique means that if [node_begin,node_end)x[nb_begin,nb_end)
 *    is in the list, [nb_begin,nb_end)x[node_begin,node_end) is not.
 *
 *  Arguments:
 *     box  <linkedListBox*>  : Box of cell linked lists
 *  Returns:
 *     node_begin <int64_t*>  : Array for receiver cell begin indexes
 *     node_end   <int64_t*>  : Array for receiver cell end indexes
 *     nb_begin   <int64_t*>  : Array for sender cell begin indexes
 *     nb_end     <int64_t*>  : Array for sender cell end indexes
 */
int setup_unique_box_pairs(linkedListBox *box,
                           int64_t *node_begin,int64_t *node_end,
                           int64_t *nb_begin,int64_t *nb_end)
{
  int64_t pair_count = 0, particle_pair_count = 0;

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){      // iterate over the cells
    int64_t node_hash=-1;
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

    if (kh_exist(box->hbegin, kbegin)){                                                         // check if the cell is valid
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));                       // get the iterator for the end index

      node_hash = kh_key(box->hbegin, kbegin);                                                  // and also get the corresponding hash

      neighbour_hash_3d(node_hash,nblist,box->width,box);                                       // fetch the list of neighbors
      for(int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){                   // and iterate over those neighbors
        if(nblist[j]>=0){                                                                       // check if the neighbor is valid
          if(kh_value(box->hbegin, kbegin) <=                                                   // and if the node_index is smaller
                     kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) ))                 // then the neighbor index
          {
            node_begin[pair_count] = kh_value(box->hbegin, kbegin);                             // get the cell's begin index
            node_end[pair_count]   = kh_value(box->hend, kend);                                 // get the cell's end index
            nb_begin[pair_count]   = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) ); // get the neighbor's begin index
            nb_end[pair_count]     = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) ); // get the neighbor's end index

            particle_pair_count += (nb_end[pair_count]-nb_begin[pair_count])*                   // just for bookeeping/profiling,
                                         (node_end[pair_count]-node_begin[pair_count]);         // annotate the number of particle pairs
            pair_count += 1;                                                                    // that are part of the computation
          }
        }
      }
    }
  }

  if(dbg)
    printf("particle_pair_count = %ld\n",particle_pair_count);

  return pair_count;
}

/*
 *  Function print_sph_particles_density:
 *    Write in disk the density values for the computed SPH particles.
 *
 *  Arguments:
 *     prefix  <const char*>  : prefix for the file name
 *     is_cll    <bool>       : boolean defined naive or cell linked list case
 *       N      <int64_t>     : Number of SPH particles
 *       h      <double>      : smoothing length
 *     seed    <long int>     : seed utilized for the PRNG
 *     runs     <int>         : Number of runs executed
 *     lsph  <SPHparticle*>   : Array of SPH particles
 *     box  <linkedListBox*>  : Box of cell linked lists
 *  Returns:
 *     0  <int>   : Error code
 *    fp <FILE*>  : Written file with particle densities
 */
int print_sph_particles_density(const char *prefix, bool is_cll, int64_t N, double h,
																long int seed, int runs, SPHparticle *lsph, linkedListBox *box){
	FILE *fp;
	char filename[2048+1];

	if(is_cll){
		sprintf(filename,
						"data/cd3d(%s,runs=%d)-P(seed=%ld,N=%ld,h=%lg)-B(Nx=%d,Ny=%d,Nz=%d)-D(%lg,%lg,%lg,%lg,%lg,%lg).csv",
						prefix,runs,seed,N,h,box->Nx,box->Ny,box->Nz,box->Xmin,box->Ymin,box->Zmin,box->Xmax,box->Ymax,box->Zmax);

		fp = fopen(filename,"w");
		fprintf(fp,"id,x,y,z,rho\n");
		for(int64_t i=0;i<N;i+=1)
			fprintf(fp,"%ld,%lf,%lf,%lf,%lf\n",lsph->id[i],lsph->x[i],lsph->y[i],lsph->z[i],lsph->rho[i]);
		fclose(fp);
	}
	else{
		sprintf(filename,
						"data/cd3d(%s,runs=%d)-P(seed=%ld,N=%ld,h=%lg)-B(Nx=%d,Ny=%d,Nz=%d)-D(%lg,%lg,%lg,%lg,%lg,%lg).csv",
						prefix,runs,seed,N,h,box->Nx,box->Ny,box->Nz,box->Xmin,box->Ymin,box->Zmin,box->Xmax,box->Ymax,box->Zmax);

		fp = fopen(filename,"w");
		fprintf(fp,"id,x,y,z,rho\n");
		for(int64_t i=0;i<N;i+=1)
			fprintf(fp,"%ld,%lf,%lf,%lf,%lf\n",lsph->id[i],lsph->x[i],lsph->y[i],lsph->z[i],lsph->rho[i]);
		fclose(fp);
	}

	return 0;
}

/*
 *  Function print_sph_particles_density:
 *    Write in disk the density values for the computed SPH particles.
 *
 *  Arguments:
 *     prefix  <const char*>  : prefix for the file name
 *     is_cll    <bool>       : boolean defined naive or cell linked list case
 *       N      <int64_t>     : Number of SPH particles
 *       h      <double>      : smoothing length
 *     seed    <long int>     : seed utilized for the PRNG
 *     runs     <int>         : Number of runs executed
 *     lsph  <SPHparticle*>   : Array of SPH particles
 *     box  <linkedListBox*>  : Box of cell linked lists
 *    times   <double*>       : timings for each of the processing steps
 *  Returns:
 *     0  <int>   : Error code
 *    fp <FILE*>  : Written file with particle densities
 */
int print_time_stats(const char *prefix, bool is_cll, int64_t N, double h,
										 long int seed, int runs, SPHparticle *lsph, linkedListBox *box,double *times){
  FILE *fp;
	char filename[1024+1];

	if(is_cll){
    const int COMPUTE_BLOCKS = 5;
    double t[COMPUTE_BLOCKS], dt[COMPUTE_BLOCKS], total_time, dtotal_time;
  	sprintf(filename,
						"data/times-(%s,runs=%d)-P(seed=%ld,N=%ld,h=%lg)-B(Nx=%d,Ny=%d,Nz=%d)-D(%lg,%lg,%lg,%lg,%lg,%lg).csv",
						prefix,runs,seed,N,h,box->Nx,box->Ny,box->Nz,box->Xmin,box->Ymin,box->Zmin,box->Xmax,box->Ymax,box->Zmax);

  	fp = fopen(filename,"w");
		fprintf(fp,"id,compute_hash_MC3D,sorting,reorder_lsph_SoA,setup_interval_hashtables,compute_density\n");
		for(int run=0;run<runs;run+=1)
			fprintf(fp,"%d,%lf,%lf,%lf,%lf,%lf\n",run,times[COMPUTE_BLOCKS*run+0],times[COMPUTE_BLOCKS*run+1],times[COMPUTE_BLOCKS*run+2],
                                                times[COMPUTE_BLOCKS*run+3],times[COMPUTE_BLOCKS*run+4]);
		fclose(fp);

  	total_time = 0.;
  	for(int k=0;k<COMPUTE_BLOCKS;k+=1){
    	t[k]=0.; dt[k]=0.;
    	for(int run=0;run<runs;run+=1)
      	t[k] += times[COMPUTE_BLOCKS*run+k];
    	t[k] /= runs;
    	for(int run=0;run<runs;run+=1)
      	dt[k] += (times[COMPUTE_BLOCKS*run+k]-t[k])*(times[COMPUTE_BLOCKS*run+k]-t[k]);
    	dt[k] /= runs;
    	dt[k] = sqrt(dt[k]);

    	total_time += t[k];
  	}

  	dtotal_time = 0.;
  	for(int run=0;run<runs;run+=1){
	    double rgm = 0.;
  	  for(int k=0;k<COMPUTE_BLOCKS;k+=1)
    	  rgm += times[COMPUTE_BLOCKS*run+k];

    	dtotal_time += (rgm-total_time)*(rgm-total_time);
  	}
  	dtotal_time /= runs;
  	dtotal_time = sqrt(dtotal_time);

  	if(dbg){
      printf("compute_hash_MC3D calc time                 : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[0],dt[0],100*t[0]/total_time,100*dt[0]/total_time);
  	  printf("qsort calc time                             : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[1],dt[1],100*t[1]/total_time,100*dt[1]/total_time);
  	  printf("reorder_lsph_SoA calc time                  : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[2],dt[2],100*t[2]/total_time,100*dt[2]/total_time);
  	  printf("setup_interval_hashtables calc time         : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[3],dt[3],100*t[3]/total_time,100*dt[3]/total_time);
  	  printf("compute_density_3d load balanced calc time  : %.5lf +- %.6lf s : %.3lg%% +- %.3lg%%\n",t[4],dt[4],100*t[4]/total_time,100*dt[4]/total_time);
  	  printf("compute_density_3d load balanced total time : %.5lf +- %.6lf s : %.3lf%%\n",total_time,dtotal_time,100.);
    }
	}
	else{
    const int COMPUTE_BLOCKS = 1;
    double t[COMPUTE_BLOCKS], dt[COMPUTE_BLOCKS], total_time, dtotal_time;
		sprintf(filename,
						"data/times-(%s,runs=%d)-P(seed=%ld,N=%ld,h=%lg)-B(Nx=%d,Ny=%d,Nz=%d)-D(%lg,%lg,%lg,%lg,%lg,%lg).csv",
						prefix,runs,seed,N,h,box->Nx,box->Ny,box->Nz,box->Xmin,box->Ymin,box->Zmin,box->Xmax,box->Ymax,box->Zmax);

  	fp = fopen(filename,"w");
		fprintf(fp,"id,compute_density\n");
		for(int run=0;run<runs;run+=1)
			fprintf(fp,"%d,%lf\n",run,times[COMPUTE_BLOCKS*run+0]);
		fclose(fp);

  	total_time = 0.;
  	for(int k=0;k<COMPUTE_BLOCKS;k+=1){
    	t[k]=0.; dt[k]=0.;
    	for(int run=0;run<runs;run+=1)
      	t[k] += times[COMPUTE_BLOCKS*run+k];
    	t[k] /= runs;
    	for(int run=0;run<runs;run+=1)
      	dt[k] += (times[COMPUTE_BLOCKS*run+k]-t[k])*(times[COMPUTE_BLOCKS*run+k]-t[k]);
    	dt[k] /= runs;
    	dt[k] = sqrt(dt[k]);

    	total_time += t[k];
  	}

  	dtotal_time = 0.;
  	for(int run=0;run<runs;run+=1){
	    double rgm = 0.;
  	  for(int k=0;k<1;k+=1)
    	  rgm += times[COMPUTE_BLOCKS*run+k];

    	dtotal_time += (rgm-total_time)*(rgm-total_time);
  	}
  	dtotal_time /= runs;
  	dtotal_time = sqrt(dtotal_time);

  	if(dbg)
      printf("compute_density_3d naive %s : %.5lf +- %.6lf s : %.3lf%%\n",prefix,total_time,dtotal_time,100.);
	}

  return 0;
}

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


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>
#include <gsl/gsl_blas_types.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_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_short * dest, const gsl_matrix_short * src);

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 int gsl_matrix_short_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, 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_equal (const gsl_matrix_short * a, const gsl_matrix_short * b);

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_scale_rows (gsl_matrix_short * a, const gsl_vector_short * x);
GSL_FUN int gsl_matrix_short_scale_columns (gsl_matrix_short * a, const gsl_vector_short * 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__ */