text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.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 "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_compute.h"

extern const double fk_bspline_32[32];
extern const double fk_bspline_128[128];
extern const double fk_bspline_1024[1024];

#pragma omp declare simd
double w_bspline_3d(double r,double h){
  const double A_d = 3./(2.*M_PI*h*h*h);
  double q=0.;

  if(r<0||h<=0.)
    exit(10);

  q = r/h;
  if(q<=1)
    return A_d*(2./3.-q*q + q*q*q/2.0);
  else if((1.<=q)&&(q<2.))
    return A_d*(1./6.)*(2.-q)*(2.-q)*(2.-q);
  else
    return 0.;
}

double w_bspline_3d_constant(double h){
  return 3./(2.*M_PI*h*h*h);
}

#pragma omp declare simd
double w_bspline_3d_simd(double q){
  double wq = 0.0;
  double wq1 = (0.6666666666666666 - q*q + 0.5*q*q*q);
  double wq2 = 0.16666666666666666*(2.-q)*(2.-q)*(2.-q);

  if(q<2.)
    wq = wq2;

  if(q<1.)
    wq = wq1;

  return wq;
}

#pragma omp declare simd
double dwdq_bspline_3d_simd(double q){
  double wq = 0.0;
  double wq1 = (1.5*q*q - 2.*q);
  double wq2 = -0.5*(2.-q)*(2.-q);

  if(q<2.)
    wq = wq2;

  if(q<1.)
    wq = wq1;

  return wq;
}

// https://stackoverflow.com/questions/56417115/efficient-integer-floor-function-in-c
// 341.333333333 = 1024./3.
#pragma omp declare simd
double w_bspline_3d_LUT_1024(double q){
  int kq = (int)(341.333333333*q);
  double phi = q-kq*0.002929688;
  double wq = phi*fk_bspline_1024[kq] + (1.-phi)*fk_bspline_1024[kq+1];

  return wq;
}

#pragma omp declare simd
double w_bspline_3d_LUT_128(double q){
  int kq = (int)(42.666666667*q);
  double phi = q-kq*0.0234375;
  double wq = phi*fk_bspline_128[kq] + (1.-phi)*fk_bspline_128[kq+1];

  return wq;
}

#pragma omp declare simd
double w_bspline_3d_LUT_32(double q){
  int kq = (int)(10.666666667*q);
  double phi = q-kq*0.09375;
  double wq = (1.-phi)*fk_bspline_32[kq] + phi*fk_bspline_32[kq+1];

  return wq;
}

double dwdq_bspline_3d(double r,double h){
  const double A_d = 3./(2.*M_PI*h*h*h*h);
  double q=0.;

  if(r<0||h<=0.)
    exit(10);

  q = r/h;
  if(q<=1)
    return A_d*q*(-2.+1.5*q);
  else if((1.<=q)&&(q<2.))
    return -A_d*0.5*(2.-q)*(2.-q);
  else
    return 0.;
}

double distance_2d(double xi,double yi,
                   double xj,double yj){
  double dist = 0.0;

  dist += (xi-xj)*(xi-xj);
  dist += (yi-yj)*(yi-yj);

  return sqrt(dist);
}

#pragma omp declare simd
double distance_3d(double xi,double yi,double zi,
                   double xj,double yj,double zj){
  double dist = 0.0;

  dist += (xi-xj)*(xi-xj);
  dist += (yi-yj)*(yi-yj);
  dist += (zi-zj)*(zi-zj);

  return sqrt(dist);
}

/**********************************************************************/

int compute_density_3d_chunk(int64_t node_begin, int64_t node_end,
                             int64_t nb_begin, int64_t nb_end,double h,
                             double* restrict x, double* restrict y,
                             double* restrict z, double* restrict nu,
                             double* restrict rho){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);

  #pragma omp parallel for
  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];
    double rhoii = 0.0;

    #pragma omp simd reduction(+:rhoii) nontemporal(rhoii)
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      rhoii += nu[jj]*w_bspline_3d_simd(q);
      //rhoii += nu[jj]*w_bspline_3d_LUT(q);
    }
    rho[ii] += rhoii*kernel_constant;
  }

  return 0;
}

int compute_density_3d_innerOmp(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int res;
  double dist = 0.0;
  khiter_t kbegin,kend;
  int64_t node_hash=-1,node_begin=0, node_end=0;
  int64_t nb_begin= 0, nb_end = 0;
  int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

  for (kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){

    if (kh_exist(box->hbegin, kbegin)){
      kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));
      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      for(int64_t ii=node_begin;ii<node_end;ii+=1)// this loop inside was the problem
        lsph->rho[ii] = 0.0;

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          //nb_hash  = nblist[j];
          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          compute_density_3d_chunk(node_begin,node_end,nb_begin,nb_end,h,
                                   lsph->x,lsph->y,lsph->z,lsph->nu,lsph->rho);
        }
      }
    }
  }

  return 0;
}

/**********************************************************************/

int compute_density_3d_chunk_noomp(int64_t node_begin, int64_t node_end,
                                   int64_t nb_begin, int64_t nb_end,double h,
                                   double* restrict x, double* restrict y,
                                   double* restrict z, double* restrict nu,
                                   double* restrict rho){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);

  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];
    double rhoii = 0.0;

    #pragma omp simd reduction(+:rhoii)
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      rhoii += nu[jj]*w_bspline_3d_simd(q);
      //rhoii += nu[jj]*w_bspline_3d_LUT(q);
    }
    rho[ii] += rhoii*kernel_constant;
  }

  return 0;
}

int compute_density_3d(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t pair_count = 0;
  const khint32_t num_boxes = kh_size(box->hbegin);
  const khiter_t hbegin_start = kh_begin(box->hbegin), hbegin_finish = kh_end(box->hbegin);

  //for (khint32_t kbegin = hbegin_start; kbegin != hbegin_finish; kbegin++)
  #pragma omp parallel for num_threads(24)
  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
    int res;
    int64_t node_hash=-1,node_begin=0, node_end=0;
    int64_t nb_begin= 0, nb_end = 0;
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

    if (kh_exist(box->hbegin, kbegin)){ // I have to call this!
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));

      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      for(int64_t ii=node_begin;ii<node_end;ii+=1)
        lsph->rho[ii] = 0.0;

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          //nb_hash  = nblist[j];

          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          pair_count += (node_end-node_begin)*(nb_end-nb_begin);

          compute_density_3d_chunk_noomp(node_begin,node_end,nb_begin,nb_end,h,
                                         lsph->x,lsph->y,lsph->z,lsph->nu,lsph->rho);
        }
      }
    }

    //printf("thread %d - %lf %lf\n",i,lsph->rho[node_begin],lsph->rho[node_end-1]);
  }

  printf("pair_count = %ld\n",pair_count);

  //for(int64_t i=0;i<N;i+=1000)
  //  printf("%ld - %lf\n",i,lsph->rho[i]);

  return 0;
}

/**********************************************************************/

int compute_density_3d_chunk_loopswapped(int64_t node_begin, int64_t node_end,
                                                int64_t node_hash,linkedListBox *box, double h,
                                                double* restrict x, double* restrict y,
                                                double* restrict z, double* restrict nu,
                                                double* restrict rho)
{
  int64_t pair_count;
  int res=0,nb_count=0,fused_count=0;
  //int64_t nb_begin= 0, nb_end = 0;
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);
  int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];
  int64_t nb_begin[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];
  int64_t nb_end[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

  res = neighbour_hash_3d(node_hash,nblist,box->width,box);
  for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
    if(nblist[j]>=0){
      nb_begin[nb_count] = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
      nb_end[nb_count]   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );
      nb_count+=1;
    }
    else{
      nb_begin[nb_count] = -1;
      nb_end[nb_count]   = -1;
    }
  }

  qsort(nb_begin,nb_count,sizeof(int64_t),compare_int64_t);
  qsort(nb_end  ,nb_count,sizeof(int64_t),compare_int64_t);

  fused_count = nb_count;
  {
    int64_t tmp;
    unsigned int j=0;
    while(j<nb_count-1){
      if(nb_begin[j] < 0){
        nb_begin[j] = nb_begin[nb_count-1]+1;
        nb_end[j]   = nb_end[nb_count-1]+1;
      }
      else if(nb_end[j]==nb_begin[j+1]){
        //printf("%ld:%ld , %ld:%d",nb_begin[j],nb_begin[j+1]);
        nb_end[j] = nb_end[j+1];

        tmp = nb_begin[j];
        nb_begin[j] = nb_begin[j+1];
        nb_begin[j+1] = tmp;

        tmp = nb_end[j];
        nb_end[j] = nb_end[j+1];
        nb_end[j+1] = tmp;

        nb_begin[j] = nb_begin[nb_count-1]+1;
        nb_end[j]   = nb_end[nb_count-1]+1;

        fused_count -= 1;
      }

      j+=1;
    }
  }

  qsort(nb_begin,nb_count,sizeof(int64_t),compare_int64_t);
  qsort(nb_end  ,nb_count,sizeof(int64_t),compare_int64_t);

  for(unsigned int j=0;j<fused_count;j+=1){
    pair_count += (node_end-node_begin)*(nb_end[j]-nb_begin[j]);
    for(int64_t ii=node_begin;ii<node_end;ii+=1){
      double xii = x[ii];
      double yii = y[ii];
      double zii = z[ii];

      #pragma omp simd
      for(int64_t jj=nb_begin[j];jj<nb_end[j];jj+=1){
        double q = 0.;

        double xij = xii-x[jj];
        double yij = yii-y[jj];
        double zij = zii-z[jj];

        q += xij*xij;
        q += yij*yij;
        q += zij*zij;

        q = sqrt(q)*inv_h;

        rho[ii] += nu[jj]*w_bspline_3d_simd(q);
      }
    }
  }

  for(int64_t ii=node_begin;ii<node_end;ii+=1)
    rho[ii] *= kernel_constant;

  return pair_count;
}

int compute_density_3d_loopswapped(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t pair_count = 0, ppc;
  const khint32_t num_boxes = kh_size(box->hbegin);
  const khiter_t hbegin_start = kh_begin(box->hbegin), hbegin_finish = kh_end(box->hbegin);

  //for (khint32_t kbegin = hbegin_start; kbegin != hbegin_finish; kbegin++)
  #pragma omp parallel for num_threads(24)
  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
    int res;
    int64_t node_hash=-1,node_begin=0, node_end=0;

    if (kh_exist(box->hbegin, kbegin)){ // I have to call this!
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));

      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      for(int64_t ii=node_begin;ii<node_end;ii+=1)
        lsph->rho[ii] = 0.0;

      ppc = compute_density_3d_chunk_loopswapped(node_begin,node_end,node_hash,box,h,
                                                 lsph->x,lsph->y,lsph->z,lsph->nu,lsph->rho);
      pair_count += ppc;
    }
  }

  printf("pair_count = %ld\n",pair_count);

  return 0;
}

/*******************************************************************/

int count_box_pairs(linkedListBox *box){
  int64_t pair_count = 0, particle_pair_count = 0;

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
    int res;
    int64_t node_hash=-1,node_begin=0, node_end=0;
    int64_t nb_begin= 0, nb_end = 0;
    int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

    if (kh_exist(box->hbegin, kbegin)){ // I have to call this!
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));

      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          //nb_hash  = nblist[j];

          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          pair_count += 1;
          particle_pair_count += (node_end-node_begin)*(nb_end-nb_begin);
        }
      }
    }
  }

  printf("unique ordered particle_pair_count = %ld\n",particle_pair_count);

  return pair_count;
}

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;

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
    int res;
    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)){ // I have to call this!
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));

      node_hash = kh_key(box->hbegin, kbegin);

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          //nb_hash  = nblist[j];

          node_begin[pair_count] = kh_value(box->hbegin, kbegin);
          node_end[pair_count]   = kh_value(box->hend, kend);
          nb_begin[pair_count]   = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end[pair_count]     = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          pair_count += 1;//(node_end-node_begin)*(nb_end-nb_begin);
        }
      }
    }
    //printf("thread %d - %lf %lf\n",i,lsph->rho[node_begin],lsph->rho[node_end-1]);
  }

  return pair_count;
}

int compute_density_3d_chunk_noomp_shift(int64_t node_begin, int64_t node_end,
                                         int64_t nb_begin, int64_t nb_end,double h,
                                         double* restrict x, double* restrict y,
                                         double* restrict z, double* restrict nu,
                                         double* restrict rho){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);

  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];
    double rhoii = 0.0;

    #pragma omp simd reduction(+:rhoii)
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      rhoii += nu[jj]*w_bspline_3d_simd(q);
    }
    rho[ii-node_begin] += rhoii*kernel_constant;
  }

  return 0;
}

int compute_density_3d_load_ballanced(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t *node_begin,*node_end,*nb_begin,*nb_end;
  int64_t max_box_pair_count = 0;

  max_box_pair_count = count_box_pairs(box);

  node_begin = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  node_end   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_begin   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_end     = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));

  setup_box_pairs(box,node_begin,node_end,nb_begin,nb_end);

  for(int64_t ii=0;ii<N;ii+=1)
    lsph->rho[ii] = 0.0;

  #pragma omp parallel for num_threads(24)
  for(size_t i=0;i<max_box_pair_count;i+=1){
    double local_rho[node_end[i]-node_begin[i]];

    for(size_t ii=0;ii<node_end[i]-node_begin[i];ii+=1)
      local_rho[ii] = 0.;

    compute_density_3d_chunk_noomp_shift(node_begin[i],node_end[i],nb_begin[i],nb_end[i],
                                         h,lsph->x,lsph->y,lsph->z,lsph->nu,local_rho);

    #pragma omp critical
    {
      for(size_t ii=node_begin[i];ii<node_end[i];ii+=1)
        lsph->rho[ii] += local_rho[ii - node_begin[i]];
    }
  }

  free(node_begin);
  free(node_end);
  free(nb_begin);
  free(nb_end);

  //for(int64_t i=0;i<N;i+=1000)
  //  printf("%ld - %lf\n",i,lsph->rho[i]);

  return 0;
}

/*******************************************************************/
/*******************************************************************/

/*
 Pit from hell
 */
int compute_density_3d_chunk_symmetrical(int64_t node_begin, int64_t node_end,
                                         int64_t nb_begin, int64_t nb_end,double h,
                                         double* restrict x, double* restrict y,
                                         double* restrict z, double* restrict nu,
                                         double* restrict rhoi, double* restrict rhoj){
  const double inv_h = 1./h;

  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];

    #pragma omp simd
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      double wij = w_bspline_3d_simd(q);

      rhoi[ii-node_begin] += nu[jj]*wij;
      rhoj[jj-nb_begin]   += nu[ii]*wij;
    }
  }

  return 0;
}

#define min(a,b) \
  ({ __typeof__ (a) _a = (a); \
     __typeof__ (b) _b = (b); \
   _a < _b ? _a : _b; })

int compute_density_3d_chunk_symm_nlo(int64_t node_begin, int64_t node_end,
                                      int64_t nb_begin, int64_t nb_end,double h,
                                      double* restrict x, double* restrict y,
                                      double* restrict z, double* restrict nu,
                                      double* restrict rhoi, double* restrict rhoj){
  const int64_t STRIP=512;
  const double inv_h = 1./h;

  /*
  for(int64_t ii=node_begin;ii<node_end;ii+=1){

    #pragma omp simd
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xii = x[ii];
      double yii = y[ii];
      double zii = z[ii];

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      double wij = w_bspline_3d_simd(q);

      rhoi[ii-node_begin] += nu[jj]*wij;
      rhoj[jj-nb_begin]   += nu[ii]*wij;
    }
  }*/

  for(int64_t i=0;i<(node_end-node_begin);i+=STRIP)
    for(int64_t j=0;j<( nb_end - nb_begin );j+=STRIP)
      for(int64_t ii=i;ii<min(node_end-node_begin,i+STRIP);ii+=1)
        #pragma omp simd
        for(int64_t jj=j;jj<min( nb_end - nb_begin ,j+STRIP);jj+=1){
          double q = 0.;

          double xij = x[ii+node_begin]-x[jj+nb_begin];
          double yij = y[ii+node_begin]-y[jj+nb_begin];
          double zij = z[ii+node_begin]-z[jj+nb_begin];

          q += xij*xij;
          q += yij*yij;
          q += zij*zij;

          q = sqrt(q)*inv_h;

          double wij = w_bspline_3d_simd(q);

          rhoi[ii] += nu[jj + nb_begin]*wij;
          rhoj[jj] += nu[ii+node_begin]*wij;
        }

  return 0;

}

int compute_density_3d_chunk_symm_colapse(int64_t node_begin, int64_t node_end,
                                          int64_t nb_begin, int64_t nb_end,double h,
                                          double* restrict x, double* restrict y,
                                          double* restrict z, double* restrict nu,
                                          double* restrict rhoi, double* restrict rhoj){
  const double inv_h = 1./h;

  #pragma omp simd
  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = x[ii]-x[jj];
      double yij = y[ii]-y[jj];
      double zij = z[ii]-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      double wij = w_bspline_3d_simd(q);

      rhoi[ii-node_begin] += nu[jj]*wij;
      rhoj[jj-nb_begin]   += nu[ii]*wij;
    }
  }

  return 0;
}

int cpr_int64_t(const void *a, const void *b){
  return ( *(int64_t*)a - *(int64_t*)b );
}

int unique_box_bounds(int64_t max_box_pair_count,
                      int64_t *node_begin, int64_t *node_end,
                      int64_t *nb_begin, int64_t *nb_end)
{
  int64_t box_skip=0;

  if(node_begin==NULL || node_end==NULL || nb_begin==NULL || nb_end==NULL)
    return -1;

  for(int64_t i=0;i<max_box_pair_count;i+=1){
    if(node_begin[i]>=0){
      for(int64_t j=i+1;j<max_box_pair_count;j+=1){
        if((node_begin[j]==nb_begin[i]) && (nb_begin[j]==node_begin[i]) ){
          node_begin[j] = -1; node_end[j] = -1;
          nb_begin[j]   = -1; nb_end[j]   = -1;

          //node_begin[j] -= -1; node_end[j] -= -1;
          //nb_begin[j]   -= -1; nb_end[j]   -= -1;

          box_skip += 1;
        }
      }
    }
  }

  //qsort(node_begin,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  //qsort(node_end  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  //qsort(nb_begin  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  //qsort(nb_end    ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);

  /*
  if(node_begin[i]<=nb_begin[i]){

      box_keep +=1;
    }

  qsort(node_begin,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(node_end  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_begin  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_end    ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);

  max_box_pair_count = box_keep; //max_box_pair_count - box_subtract;

  for(int64_t i=0;i<max_box_pair_count;i+=1){
    node_begin[i] *= -1; node_end[i] *= -1;
    nb_begin[i]   *= -1; nb_end[i]   *= -1;
  }

  qsort(node_begin,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(node_end  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_begin  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_end    ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);*/

  return box_skip;
}

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;

  for (khint32_t kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){
    int res;
    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)){ // I have to call this!
      khint32_t kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));

      node_hash = kh_key(box->hbegin, kbegin);

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          //nb_hash  = nblist[j];

          if(kh_value(box->hbegin, kbegin) <= kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) ))
          {
            node_begin[pair_count] = kh_value(box->hbegin, kbegin);
            node_end[pair_count]   = kh_value(box->hend, kend);
            nb_begin[pair_count]   = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
            nb_end[pair_count]     = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

            pair_count += 1;//(node_end-node_begin)*(nb_end-nb_begin);
          }
        }
      }
    }
    //printf("thread %d - %lf %lf\n",i,lsph->rho[node_begin],lsph->rho[node_end-1]);
  }

  return pair_count;
}

int compute_density_3d_symmetrical_load_ballance(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t *node_begin,*node_end,*nb_begin,*nb_end;
  int64_t max_box_pair_count = 0, particle_pair_count = 0;
  int64_t box_skip = 0;
  const double kernel_constant = w_bspline_3d_constant(h);

  max_box_pair_count = count_box_pairs(box);
  printf("max_box_pair_count = %ld\n",max_box_pair_count);

  node_begin = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  node_end   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_begin   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_end     = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));

  max_box_pair_count = setup_unique_box_pairs(box,node_begin,node_end,nb_begin,nb_end);//setup_box_pairs(box,node_begin,node_end,nb_begin,nb_end);
  printf("unique max_box_pair_count = %ld\n",max_box_pair_count);

  double avg_nb_len = 0.0;
  for(int64_t i=0;i<max_box_pair_count;i+=1)
    avg_nb_len += (nb_end[i]-nb_begin[i]);
  avg_nb_len /= max_box_pair_count;
  printf("avg_nb_len = %lf\n",avg_nb_len);

  for(int64_t i=0;i<max_box_pair_count;i+=1)
    particle_pair_count += (node_end[i]-node_begin[i])*(nb_end[i]-nb_begin[i]);
  printf("unique unordered particle_pair_count = %ld\n",particle_pair_count);

  //box_skip = unique_box_bounds(max_box_pair_count,node_begin,node_end,nb_begin,nb_end);

  for(int64_t ii=0;ii<N;ii+=1)
    lsph->rho[ii] = 0.0;

  #pragma omp parallel for schedule(dynamic,5) proc_bind(master)
  for(size_t i=0;i<max_box_pair_count;i+=1){

    double local_rhoi[node_end[i] - node_begin[i]];
    double local_rhoj[  nb_end[i] -   nb_begin[i]];

    for(size_t ii=0;ii<node_end[i]-node_begin[i];ii+=1)
      local_rhoi[ii] = 0.;

    for(size_t ii=0;ii<nb_end[i]-nb_begin[i];ii+=1)
      local_rhoj[ii] = 0.;

    compute_density_3d_chunk_symmetrical(node_begin[i],node_end[i],nb_begin[i],nb_end[i],h,
                                         lsph->x,lsph->y,lsph->z,lsph->nu,local_rhoi,local_rhoj);

    //compute_density_3d_chunk_symm_nlo(node_begin[i],node_end[i],nb_begin[i],nb_end[i],h,
    //                                   lsph->x,lsph->y,lsph->z,lsph->nu,local_rhoi,local_rhoj);

    //compute_density_3d_chunk_symm_colapse(node_begin[i],node_end[i],nb_begin[i],nb_end[i],h,
    //                                      lsph->x,lsph->y,lsph->z,lsph->nu,local_rhoi,local_rhoj);

    #pragma omp critical
    {

      for(size_t ii=node_begin[i];ii<node_end[i];ii+=1){
        lsph->rho[ii] += kernel_constant*local_rhoi[ii - node_begin[i]];
      }

      if(node_begin[i] != nb_begin[i])
        for(size_t ii=nb_begin[i];ii<nb_end[i];ii+=1){
          lsph->rho[ii] += kernel_constant*local_rhoj[ii - nb_begin[i]];
        }
    }
  }

  free(node_begin);
  free(node_end);
  free(nb_begin);
  free(nb_end);

  return 0;
}

int compute_density_3d_symmetrical_lb_branching(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t *node_begin,*node_end,*nb_begin,*nb_end;
  int64_t max_box_pair_count = 0, particle_pair_count = 0;
  int64_t box_skip = 0;
  const double kernel_constant = w_bspline_3d_constant(h);

  max_box_pair_count = count_box_pairs(box);
  printf("max_box_pair_count = %ld\n",max_box_pair_count);

  node_begin = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  node_end   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_begin   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_end     = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));

  max_box_pair_count = setup_box_pairs(box,node_begin,node_end,nb_begin,nb_end);

  box_skip = unique_box_bounds(max_box_pair_count,node_begin,node_end,nb_begin,nb_end);

  for(int64_t i=0;i<max_box_pair_count;i+=1)
    if(node_end[i] >= 0)
      particle_pair_count += (node_end[i]-node_begin[i])*(nb_end[i]-nb_begin[i]);
  printf("unique unordered particle_pair_count = %ld\n",particle_pair_count);

  for(int64_t ii=0;ii<N;ii+=1)
    lsph->rho[ii] = 0.0;

  #pragma omp parallel for
  for(size_t i=0;i<max_box_pair_count;i+=1){
    if(node_end[i]<0)
      continue;

    double local_rhoi[node_end[i] - node_begin[i]];
    double local_rhoj[  nb_end[i] -   nb_begin[i]];

    for(size_t ii=0;ii<node_end[i]-node_begin[i];ii+=1)
      local_rhoi[ii] = 0.;

    for(size_t ii=0;ii<nb_end[i]-nb_begin[i];ii+=1)
      local_rhoj[ii] = 0.;

    compute_density_3d_chunk_symmetrical(node_begin[i],node_end[i],nb_begin[i],nb_end[i],h,
                                         lsph->x,lsph->y,lsph->z,lsph->nu,local_rhoi,local_rhoj);

    #pragma omp critical
    {
      for(size_t ii=node_begin[i];ii<node_end[i];ii+=1)
        lsph->rho[ii] += kernel_constant*local_rhoi[ii - node_begin[i]];

      if(node_begin[i] != nb_begin[i])
        for(size_t ii=nb_begin[i];ii<nb_end[i];ii+=1)
          lsph->rho[ii] += kernel_constant*local_rhoj[ii - nb_begin[i]];

      /*
      if(node_begin[i]!=nb_begin[i])
        for(size_t ii=nb_begin[i];ii<nb_end[i];ii+=1)
          lsph->rho[ii] += kernel_constant*local_rhoj[ii - nb_begin[i]];*/
    }
  }

  free(node_begin);
  free(node_end);
  free(nb_begin);
  free(nb_end);

  return 0;
}

/*******************************************************************/
/*******************************************************************/

/*
int compute_density_3d_chunk_symmetrical(int64_t node_begin, int64_t node_end,
                                         int64_t nb_begin, int64_t nb_end,double h,
                                         double* restrict x, double* restrict y,
                                         double* restrict z, double* restrict nu,
                                         double* restrict rho)
{
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);

  double rhoi[node_end-node_begin];
  double rhoj[nb_end-nb_begin];

  for(int64_t ii=0;ii<node_end-node_begin;ii+=1)
    rhoi[ii] = 0.;

  for(int64_t jj=0;jj<nb_end-nb_begin;jj+=1)
    rhoj[jj] = 0.;

  for(int64_t ii=node_begin;ii<node_end;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];

    #pragma omp simd
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
      double q = 0.;

      double xij = xii-x[jj];
      double yij = yii-y[jj];
      double zij = zii-z[jj];

      q += xij*xij;
      q += yij*yij;
      q += zij*zij;

      q = sqrt(q)*inv_h;

      double wij = w_bspline_3d_simd(q);

      rhoi[ii-node_begin] += nu[jj]*wij;
      rhoj[jj-nb_begin]   += nu[ii]*wij;
    }
  }

    for(int64_t ii=node_begin;ii<node_end;ii+=1)
      rho[ii] += rhoi[ii-node_begin]*kernel_constant;

  if(nb_begin == node_begin)
    return 0;

    for(int64_t jj=nb_begin;jj<nb_end;jj+=1)
      rho[jj] += rhoj[jj-nb_begin]*kernel_constant;

  return 0;
}*/

/*
int cpr_int64_t(const void *a, const void *b){
  return ( *(int64_t*)a - *(int64_t*)b );
} */

/*
int unique_box_bounds(int64_t max_box_pair_count,
                      int64_t *node_begin, int64_t *node_end,
                      int64_t *nb_begin, int64_t *nb_end)
{
  int64_t box_keep=0;

  if(node_begin==NULL || node_end==NULL || nb_begin==NULL || nb_end==NULL)
    return -1;

  for(int64_t i=0;i<max_box_pair_count;i+=1)
    if(node_begin[i]<=nb_begin[i]){
      node_begin[i] *= -1; node_end[i] *= -1;
      nb_begin[i]   *= -1; nb_end[i]   *= -1;
      box_keep +=1;
    }

  qsort(node_begin,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(node_end  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_begin  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_end    ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);

  max_box_pair_count = box_keep; //max_box_pair_count - box_subtract;

  for(int64_t i=0;i<max_box_pair_count;i+=1){
    node_begin[i] *= -1; node_end[i] *= -1;
    nb_begin[i]   *= -1; nb_end[i]   *= -1;
  }

  qsort(node_begin,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(node_end  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_begin  ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);
  qsort(nb_end    ,max_box_pair_count, sizeof(int64_t),cpr_int64_t);

  return max_box_pair_count;
}

int compute_density_3d_symmetrical_lb(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t *node_begin,*node_end,*nb_begin,*nb_end;
  int64_t max_box_pair_count = 0;

  max_box_pair_count = count_box_pairs(box);

  node_begin = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  node_end   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_begin   = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));
  nb_end     = (int64_t*)malloc(max_box_pair_count*sizeof(int64_t));

  setup_box_pairs(box,node_begin,node_end,nb_begin,nb_end);

  max_box_pair_count = unique_box_bounds(max_box_pair_count,node_begin,node_end,nb_begin,nb_end);

  for(int64_t ii=0;ii<N;ii+=1)
    lsph->rho[ii] = 0.0;

  #pragma omp parallel for num_threads(24)
  for(size_t i=0;i<max_box_pair_count;i+=1){
    if(node_begin[i] > nb_begin[i])
      continue;

    compute_density_3d_chunk_symmetrical(node_begin[i],node_end[i],nb_begin[i],nb_end[i],
                                         h,lsph->x,lsph->y,lsph->z,lsph->nu,lsph->rho);
    }

  free(node_begin);
  free(node_end);
  free(nb_begin);
  free(nb_end);

  //for(int64_t i=0;i<N;i+=1000)
  //  printf("%ld - %lf\n",i,lsph->rho[i]);

  return 0;
}*/

/*******************************************************************/
/*******************************************************************/
/*******************************************************************/

#pragma omp declare simd
double pressure_from_density(double rho){
  double p = cbrt(rho);
  p = 0.5*p*rho;
  return p;
}

#pragma omp declare simd
double gamma_from_u(double ux,double uy,double uz){
  double gamma = 1.0;
  gamma += ux*ux;
  gamma += uy*uy;
  gamma += uz*uz;

  gamma = sqrt(gamma);

  return gamma;
}

/*
int compute_force_3d(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int err, res;
  double dist = 0.0;
  khiter_t kbegin,kend;
  int64_t node_hash=-1,node_begin=0, node_end=0;
  int64_t nb_hash=-1  , nb_begin= 0, nb_end = 0;
  int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

  for (kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){

    if (kh_exist(box->hbegin, kbegin)){
      kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));
      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      for(int64_t ii=node_begin;ii<node_end;ii+=1)// this loop inside was the problem
        lsph->rho[ii] = 0.0;

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          nb_hash  = nblist[j];
          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          for(int64_t ii=node_begin;ii<node_end;ii+=1){ // this loop inside was the problem
            for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
              dist = distance_3d(lsph->x[i],lsph->y[i],lsph->z[i],
                                 lsph->x[j],lsph->y[j],lsph->z[j]);
              lsph->rho[ii] += (lsph->nu[jj])*(box->w(dist,h));
            }
          }
        }
      }
    }
  }

  return 0;
} */

/*
int compute_force_3d(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int err, res;
  double dist = 0.0;
  khiter_t kbegin,kend;
  int64_t node_hash=-1,node_begin=0, node_end=0;
  int64_t nb_hash=-1  , nb_begin= 0, nb_end = 0;
  int64_t nblist[(2*box->width+1)*(2*box->width+1)*(2*box->width+1)];

  for (kbegin = kh_begin(box->hbegin); kbegin != kh_end(box->hbegin); kbegin++){

    if (kh_exist(box->hbegin, kbegin)){
      kend = kh_get(1, box->hend, kh_key(box->hbegin, kbegin));
      node_hash = kh_key(box->hbegin, kbegin);
      node_begin = kh_value(box->hbegin, kbegin);
      node_end   = kh_value(box->hend, kend);

      for(int64_t ii=node_begin;ii<node_end;ii+=1){
        lsph[ii].F.x = 0.0;
        lsph[ii].F.y = 0.0;
        lsph[ii].F.z = 0.0;
      }

      res = neighbour_hash_3d(node_hash,nblist,box->width,box);
      for(unsigned int j=0;j<(2*box->width+1)*(2*box->width+1)*(2*box->width+1);j+=1){
        if(nblist[j]>=0){
          nb_hash  = nblist[j];
          nb_begin = kh_value(box->hbegin, kh_get(0, box->hbegin, nblist[j]) );
          nb_end   = kh_value(box->hend  , kh_get(1, box->hend  , nblist[j]) );

          for(int64_t ii=node_begin;ii<node_end;ii+=1){ // this loop inside was the problem
            for(int64_t jj=nb_begin;jj<nb_end;jj+=1){
              dist = double4_distance_3d(&(lsph[ii].r),&(lsph[jj].r));
              = ();
              lsph[ii].F.x += ((lsph[ii].r.x-lsph[jj].r.x)/dist);
              lsph[ii].F.y += ;
              lsph[ii].F.z += ;
            }
          }
        }
      }
    }
  }

  return 0;
}
*/

const double fk_bspline_32[32] = {
       6.66666667e-01, 6.57754579e-01, 6.32830944e-01, 5.94614705e-01,
       5.45824802e-01, 4.89180177e-01, 4.27399774e-01, 3.63202533e-01,
       2.99307397e-01, 2.38433308e-01, 1.83299207e-01, 1.36445011e-01,
       9.83294731e-02, 6.80686561e-02, 4.47562463e-02, 2.74859298e-02,
       1.53513925e-02, 7.44632048e-03, 2.86439976e-03, 6.99316348e-04,
       4.47562463e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00};

const double fk_bspline_128[128] = {
       6.66666667e-01, 6.66115256e-01, 6.64487387e-01, 6.61822602e-01,
       6.58160445e-01, 6.53540459e-01, 6.48002188e-01, 6.41585176e-01,
       6.34328964e-01, 6.26273098e-01, 6.17457119e-01, 6.07920573e-01,
       5.97703001e-01, 5.86843948e-01, 5.75382957e-01, 5.63359570e-01,
       5.50813333e-01, 5.37783787e-01, 5.24310476e-01, 5.10432945e-01,
       4.96190735e-01, 4.81623391e-01, 4.66770456e-01, 4.51671473e-01,
       4.36365986e-01, 4.20893537e-01, 4.05293671e-01, 3.89605931e-01,
       3.73869861e-01, 3.58125002e-01, 3.42410900e-01, 3.26767097e-01,
       3.11233137e-01, 2.95848563e-01, 2.80652918e-01, 2.65685747e-01,
       2.50986591e-01, 2.36594995e-01, 2.22550503e-01, 2.08892657e-01,
       1.95661000e-01, 1.82895077e-01, 1.70634431e-01, 1.58916000e-01,
       1.47746458e-01, 1.37112949e-01, 1.27002291e-01, 1.17401303e-01,
       1.08296805e-01, 9.96756140e-02, 9.15245505e-02, 8.38304328e-02,
       7.65800797e-02, 6.97603101e-02, 6.33579430e-02, 5.73597971e-02,
       5.17526914e-02, 4.65234448e-02, 4.16588760e-02, 3.71458040e-02,
       3.29710476e-02, 2.91214257e-02, 2.55837572e-02, 2.23448610e-02,
       1.93915558e-02, 1.67106607e-02, 1.42889945e-02, 1.21133759e-02,
       1.01706240e-02, 8.44755758e-03, 6.93099549e-03, 5.60775662e-03,
       4.46465985e-03, 3.48852404e-03, 2.66616806e-03, 1.98441079e-03,
       1.43007110e-03, 9.89967859e-04, 6.50919937e-04, 3.99746206e-04,
       2.23265538e-04, 1.08296805e-04, 4.16588760e-05, 1.01706240e-05,
       6.50919937e-07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00};

const double fk_bspline_1024[1024] = {
       6.66666667e-01, 6.66658079e-01, 6.66632368e-01, 6.66589608e-01,
       6.66529876e-01, 6.66453246e-01, 6.66359796e-01, 6.66249599e-01,
       6.66122732e-01, 6.65979271e-01, 6.65819291e-01, 6.65642868e-01,
       6.65450077e-01, 6.65240994e-01, 6.65015696e-01, 6.64774257e-01,
       6.64516753e-01, 6.64243260e-01, 6.63953853e-01, 6.63648609e-01,
       6.63327602e-01, 6.62990909e-01, 6.62638605e-01, 6.62270765e-01,
       6.61887466e-01, 6.61488783e-01, 6.61074792e-01, 6.60645569e-01,
       6.60201188e-01, 6.59741726e-01, 6.59267259e-01, 6.58777861e-01,
       6.58273610e-01, 6.57754579e-01, 6.57220846e-01, 6.56672485e-01,
       6.56109573e-01, 6.55532184e-01, 6.54940396e-01, 6.54334282e-01,
       6.53713920e-01, 6.53079384e-01, 6.52430750e-01, 6.51768095e-01,
       6.51091493e-01, 6.50401020e-01, 6.49696752e-01, 6.48978765e-01,
       6.48247134e-01, 6.47501935e-01, 6.46743244e-01, 6.45971135e-01,
       6.45185686e-01, 6.44386971e-01, 6.43575066e-01, 6.42750048e-01,
       6.41911990e-01, 6.41060970e-01, 6.40197063e-01, 6.39320344e-01,
       6.38430889e-01, 6.37528774e-01, 6.36614075e-01, 6.35686866e-01,
       6.34747225e-01, 6.33795226e-01, 6.32830944e-01, 6.31854457e-01,
       6.30865839e-01, 6.29865166e-01, 6.28852514e-01, 6.27827959e-01,
       6.26791575e-01, 6.25743439e-01, 6.24683626e-01, 6.23612213e-01,
       6.22529274e-01, 6.21434885e-01, 6.20329123e-01, 6.19212062e-01,
       6.18083778e-01, 6.16944347e-01, 6.15793845e-01, 6.14632348e-01,
       6.13459930e-01, 6.12276668e-01, 6.11082637e-01, 6.09877913e-01,
       6.08662571e-01, 6.07436688e-01, 6.06200339e-01, 6.04953599e-01,
       6.03696545e-01, 6.02429251e-01, 6.01151794e-01, 5.99864249e-01,
       5.98566692e-01, 5.97259199e-01, 5.95941844e-01, 5.94614705e-01,
       5.93277856e-01, 5.91931373e-01, 5.90575332e-01, 5.89209808e-01,
       5.87834877e-01, 5.86450616e-01, 5.85057098e-01, 5.83654401e-01,
       5.82242599e-01, 5.80821769e-01, 5.79391986e-01, 5.77953326e-01,
       5.76505864e-01, 5.75049676e-01, 5.73584838e-01, 5.72111425e-01,
       5.70629514e-01, 5.69139179e-01, 5.67640496e-01, 5.66133541e-01,
       5.64618390e-01, 5.63095118e-01, 5.61563801e-01, 5.60024515e-01,
       5.58477335e-01, 5.56922337e-01, 5.55359597e-01, 5.53789190e-01,
       5.52211192e-01, 5.50625678e-01, 5.49032725e-01, 5.47432408e-01,
       5.45824802e-01, 5.44209983e-01, 5.42588027e-01, 5.40959010e-01,
       5.39323007e-01, 5.37680094e-01, 5.36030346e-01, 5.34373840e-01,
       5.32710650e-01, 5.31040853e-01, 5.29364524e-01, 5.27681738e-01,
       5.25992573e-01, 5.24297102e-01, 5.22595402e-01, 5.20887548e-01,
       5.19173617e-01, 5.17453683e-01, 5.15727823e-01, 5.13996112e-01,
       5.12258626e-01, 5.10515440e-01, 5.08766630e-01, 5.07012271e-01,
       5.05252441e-01, 5.03487213e-01, 5.01716663e-01, 4.99940869e-01,
       4.98159904e-01, 4.96373845e-01, 4.94582767e-01, 4.92786746e-01,
       4.90985857e-01, 4.89180177e-01, 4.87369781e-01, 4.85554745e-01,
       4.83735144e-01, 4.81911054e-01, 4.80082550e-01, 4.78249708e-01,
       4.76412605e-01, 4.74571315e-01, 4.72725914e-01, 4.70876478e-01,
       4.69023083e-01, 4.67165804e-01, 4.65304717e-01, 4.63439897e-01,
       4.61571420e-01, 4.59699362e-01, 4.57823799e-01, 4.55944806e-01,
       4.54062459e-01, 4.52176833e-01, 4.50288004e-01, 4.48396048e-01,
       4.46501040e-01, 4.44603057e-01, 4.42702173e-01, 4.40798465e-01,
       4.38892008e-01, 4.36982877e-01, 4.35071149e-01, 4.33156899e-01,
       4.31240203e-01, 4.29321136e-01, 4.27399774e-01, 4.25476193e-01,
       4.23550468e-01, 4.21622675e-01, 4.19692890e-01, 4.17761188e-01,
       4.15827645e-01, 4.13892336e-01, 4.11955338e-01, 4.10016726e-01,
       4.08076576e-01, 4.06134962e-01, 4.04191962e-01, 4.02247650e-01,
       4.00302103e-01, 3.98355395e-01, 3.96407603e-01, 3.94458803e-01,
       3.92509069e-01, 3.90558477e-01, 3.88607104e-01, 3.86655025e-01,
       3.84702315e-01, 3.82749050e-01, 3.80795307e-01, 3.78841159e-01,
       3.76886684e-01, 3.74931957e-01, 3.72977053e-01, 3.71022048e-01,
       3.69067018e-01, 3.67112038e-01, 3.65157185e-01, 3.63202533e-01,
       3.61248159e-01, 3.59294138e-01, 3.57340545e-01, 3.55387457e-01,
       3.53434949e-01, 3.51483097e-01, 3.49531976e-01, 3.47581662e-01,
       3.45632230e-01, 3.43683758e-01, 3.41736319e-01, 3.39789989e-01,
       3.37844845e-01, 3.35900962e-01, 3.33958416e-01, 3.32017282e-01,
       3.30077636e-01, 3.28139553e-01, 3.26203110e-01, 3.24268382e-01,
       3.22335444e-01, 3.20404373e-01, 3.18475243e-01, 3.16548131e-01,
       3.14623112e-01, 3.12700262e-01, 3.10779657e-01, 3.08861372e-01,
       3.06945483e-01, 3.05032065e-01, 3.03121194e-01, 3.01212946e-01,
       2.99307397e-01, 2.97404622e-01, 2.95504697e-01, 2.93607697e-01,
       2.91713698e-01, 2.89822776e-01, 2.87935006e-01, 2.86050465e-01,
       2.84169227e-01, 2.82291369e-01, 2.80416966e-01, 2.78546093e-01,
       2.76678827e-01, 2.74815243e-01, 2.72955417e-01, 2.71099424e-01,
       2.69247340e-01, 2.67399241e-01, 2.65555202e-01, 2.63715299e-01,
       2.61879608e-01, 2.60048204e-01, 2.58221163e-01, 2.56398561e-01,
       2.54580473e-01, 2.52766975e-01, 2.50958142e-01, 2.49154051e-01,
       2.47354777e-01, 2.45560395e-01, 2.43770982e-01, 2.41986612e-01,
       2.40207362e-01, 2.38433308e-01, 2.36664524e-01, 2.34901086e-01,
       2.33143071e-01, 2.31390554e-01, 2.29643610e-01, 2.27902316e-01,
       2.26166746e-01, 2.24436977e-01, 2.22713084e-01, 2.20995143e-01,
       2.19283229e-01, 2.17577418e-01, 2.15877786e-01, 2.14184409e-01,
       2.12497361e-01, 2.10816720e-01, 2.09142560e-01, 2.07474956e-01,
       2.05813986e-01, 2.04159724e-01, 2.02512246e-01, 2.00871628e-01,
       1.99237945e-01, 1.97611273e-01, 1.95991688e-01, 1.94379265e-01,
       1.92774080e-01, 1.91176209e-01, 1.89585728e-01, 1.88002711e-01,
       1.86427235e-01, 1.84859375e-01, 1.83299207e-01, 1.81746806e-01,
       1.80202249e-01, 1.78665611e-01, 1.77136968e-01, 1.75616394e-01,
       1.74103967e-01, 1.72599761e-01, 1.71103853e-01, 1.69616317e-01,
       1.68137230e-01, 1.66666667e-01, 1.65204687e-01, 1.63751281e-01,
       1.62306426e-01, 1.60870094e-01, 1.59442261e-01, 1.58022902e-01,
       1.56611992e-01, 1.55209505e-01, 1.53815416e-01, 1.52429700e-01,
       1.51052331e-01, 1.49683285e-01, 1.48322536e-01, 1.46970060e-01,
       1.45625830e-01, 1.44289821e-01, 1.42962009e-01, 1.41642368e-01,
       1.40330873e-01, 1.39027499e-01, 1.37732220e-01, 1.36445011e-01,
       1.35165848e-01, 1.33894704e-01, 1.32631555e-01, 1.31376375e-01,
       1.30129139e-01, 1.28889822e-01, 1.27658399e-01, 1.26434845e-01,
       1.25219133e-01, 1.24011240e-01, 1.22811140e-01, 1.21618807e-01,
       1.20434217e-01, 1.19257343e-01, 1.18088162e-01, 1.16926648e-01,
       1.15772775e-01, 1.14626518e-01, 1.13487852e-01, 1.12356752e-01,
       1.11233193e-01, 1.10117149e-01, 1.09008596e-01, 1.07907507e-01,
       1.06813859e-01, 1.05727624e-01, 1.04648779e-01, 1.03577299e-01,
       1.02513157e-01, 1.01456328e-01, 1.00406788e-01, 9.93645117e-02,
       9.83294731e-02, 9.73016473e-02, 9.62810090e-02, 9.52675330e-02,
       9.42611942e-02, 9.32619673e-02, 9.22698271e-02, 9.12847483e-02,
       9.03067058e-02, 8.93356743e-02, 8.83716286e-02, 8.74145435e-02,
       8.64643938e-02, 8.55211542e-02, 8.45847996e-02, 8.36553047e-02,
       8.27326442e-02, 8.18167931e-02, 8.09077259e-02, 8.00054177e-02,
       7.91098430e-02, 7.82209767e-02, 7.73387936e-02, 7.64632684e-02,
       7.55943760e-02, 7.47320911e-02, 7.38763885e-02, 7.30272430e-02,
       7.21846293e-02, 7.13485223e-02, 7.05188966e-02, 6.96957272e-02,
       6.88789888e-02, 6.80686561e-02, 6.72647039e-02, 6.64671071e-02,
       6.56758404e-02, 6.48908785e-02, 6.41121963e-02, 6.33397686e-02,
       6.25735701e-02, 6.18135756e-02, 6.10597598e-02, 6.03120976e-02,
       5.95705638e-02, 5.88351331e-02, 5.81057803e-02, 5.73824802e-02,
       5.66652075e-02, 5.59539371e-02, 5.52486438e-02, 5.45493022e-02,
       5.38558872e-02, 5.31683736e-02, 5.24867361e-02, 5.18109496e-02,
       5.11409888e-02, 5.04768285e-02, 4.98184434e-02, 4.91658084e-02,
       4.85188982e-02, 4.78776876e-02, 4.72421515e-02, 4.66122645e-02,
       4.59880014e-02, 4.53693371e-02, 4.47562463e-02, 4.41487038e-02,
       4.35466844e-02, 4.29501628e-02, 4.23591138e-02, 4.17735123e-02,
       4.11933330e-02, 4.06185507e-02, 4.00491401e-02, 3.94850760e-02,
       3.89263333e-02, 3.83728867e-02, 3.78247109e-02, 3.72817808e-02,
       3.67440712e-02, 3.62115568e-02, 3.56842123e-02, 3.51620127e-02,
       3.46449326e-02, 3.41329469e-02, 3.36260303e-02, 3.31241576e-02,
       3.26273035e-02, 3.21354430e-02, 3.16485507e-02, 3.11666014e-02,
       3.06895699e-02, 3.02174310e-02, 2.97501595e-02, 2.92877301e-02,
       2.88301177e-02, 2.83772970e-02, 2.79292427e-02, 2.74859298e-02,
       2.70473328e-02, 2.66134267e-02, 2.61841863e-02, 2.57595862e-02,
       2.53396013e-02, 2.49242063e-02, 2.45133761e-02, 2.41070854e-02,
       2.37053089e-02, 2.33080216e-02, 2.29151981e-02, 2.25268132e-02,
       2.21428418e-02, 2.17632586e-02, 2.13880383e-02, 2.10171558e-02,
       2.06505858e-02, 2.02883032e-02, 1.99302826e-02, 1.95764990e-02,
       1.92269270e-02, 1.88815414e-02, 1.85403171e-02, 1.82032287e-02,
       1.78702512e-02, 1.75413592e-02, 1.72165275e-02, 1.68957310e-02,
       1.65789443e-02, 1.62661424e-02, 1.59572999e-02, 1.56523917e-02,
       1.53513925e-02, 1.50542771e-02, 1.47610203e-02, 1.44715968e-02,
       1.41859815e-02, 1.39041492e-02, 1.36260745e-02, 1.33517323e-02,
       1.30810974e-02, 1.28141446e-02, 1.25508485e-02, 1.22911841e-02,
       1.20351261e-02, 1.17826493e-02, 1.15337284e-02, 1.12883382e-02,
       1.10464536e-02, 1.08080492e-02, 1.05731000e-02, 1.03415805e-02,
       1.01134657e-02, 9.88873037e-03, 9.66734920e-03, 9.44929700e-03,
       9.23454856e-03, 9.02307866e-03, 8.81486208e-03, 8.60987360e-03,
       8.40808799e-03, 8.20948005e-03, 8.01402454e-03, 7.82169626e-03,
       7.63246998e-03, 7.44632048e-03, 7.26322254e-03, 7.08315094e-03,
       6.90608047e-03, 6.73198590e-03, 6.56084202e-03, 6.39262360e-03,
       6.22730542e-03, 6.06486228e-03, 5.90526893e-03, 5.74850018e-03,
       5.59453079e-03, 5.44333554e-03, 5.29488923e-03, 5.14916663e-03,
       5.00614251e-03, 4.86579166e-03, 4.72808886e-03, 4.59300890e-03,
       4.46052654e-03, 4.33061658e-03, 4.20325378e-03, 4.07841294e-03,
       3.95606884e-03, 3.83619624e-03, 3.71876994e-03, 3.60376471e-03,
       3.49115534e-03, 3.38091660e-03, 3.27302328e-03, 3.16745016e-03,
       3.06417201e-03, 2.96316362e-03, 2.86439976e-03, 2.76785523e-03,
       2.67350479e-03, 2.58132323e-03, 2.49128533e-03, 2.40336587e-03,
       2.31753963e-03, 2.23378139e-03, 2.15206594e-03, 2.07236804e-03,
       1.99466249e-03, 1.91892406e-03, 1.84512753e-03, 1.77324769e-03,
       1.70325931e-03, 1.63513718e-03, 1.56885607e-03, 1.50439077e-03,
       1.44171605e-03, 1.38080670e-03, 1.32163749e-03, 1.26418322e-03,
       1.20841865e-03, 1.15431857e-03, 1.10185776e-03, 1.05101100e-03,
       1.00175307e-03, 9.54058747e-04, 9.07902817e-04, 8.63260059e-04,
       8.20105252e-04, 7.78413178e-04, 7.38158617e-04, 6.99316348e-04,
       6.61861154e-04, 6.25767814e-04, 5.91011108e-04, 5.57565818e-04,
       5.25406723e-04, 4.94508604e-04, 4.64846242e-04, 4.36394418e-04,
       4.09127910e-04, 3.83021501e-04, 3.58049970e-04, 3.34188099e-04,
       3.11410666e-04, 2.89692454e-04, 2.69008242e-04, 2.49332811e-04,
       2.30640942e-04, 2.12907414e-04, 1.96107009e-04, 1.80214506e-04,
       1.65204687e-04, 1.51052331e-04, 1.37732220e-04, 1.25219133e-04,
       1.13487852e-04, 1.02513157e-04, 9.22698271e-05, 8.27326442e-05,
       7.38763885e-05, 6.56758404e-05, 5.81057803e-05, 5.11409888e-05,
       4.47562463e-05, 3.89263333e-05, 3.36260303e-05, 2.88301177e-05,
       2.45133761e-05, 2.06505858e-05, 1.72165275e-05, 1.41859815e-05,
       1.15337284e-05, 9.23454856e-06, 7.26322254e-06, 5.59453079e-06,
       4.20325378e-06, 3.06417201e-06, 2.15206594e-06, 1.44171605e-06,
       9.07902817e-07, 5.25406723e-07, 2.69008242e-07, 1.13487852e-07,
       3.36260303e-08, 4.20325378e-09, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00};

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
/**
 * (c) Author:    Woongkyu Jee, woong.jee.16@ucl.ac.uk, wldndrb1@gmail.com
 * Created:   02.06.2019 ~
 *
 * University College London, Department of Chemistry
 **/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include"sp_cluster_type.h"

#define SP_SUPPORT_TRUE 1
#define SP_SUPPORT_FALSE -1

#define DEBUG_SUPPORT

/* Get the Lowest energy state */
int sp_cluster_support_get_lowest_state( gsl_vector* v )
{
    int Return = 0;
    double min;

    min = gsl_vector_get(v,0);
    for(int i=0;i<4;i++)
    {
        if( gsl_vector_get(v,i) < min  )
        {
            min = gsl_vector_get(v,i);
            Return = i;
        }
    }
    return Return;
    // Returns the index of element in gsl_vector* eval
}

void sp_cluster_support_sign_gs_eigenvector( void* sp_sys_void )
{
    sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    int low_idx;

    for(int i=0;i<sp_sys->number_of_sp_ion;i++)
    {
        low_idx = sp_cluster_support_get_lowest_state(sp_sys->sp_ion[i].eigen_value);
        if( gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,0,low_idx) < 0. )          // eigen_vector : matrix(4,4) // eigen_vector_gs : vector(4)
        {
            gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,0,low_idx,-1.*gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,0,low_idx));
            gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,1,low_idx,-1.*gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,1,low_idx));
            gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,2,low_idx,-1.*gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,2,low_idx));
            gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,3,low_idx,-1.*gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,3,low_idx));
        }
    }
    return;
}

void sp_cluster_support_load_gs_eigenvector( void* sp_sys_void )
{
    sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    int low_idx;
    for(int i=0;i<sp_sys->number_of_sp_ion;i++)
    {   low_idx = sp_cluster_support_get_lowest_state(sp_sys->sp_ion[i].eigen_value);

        gsl_vector_set(sp_sys->sp_ion[i].eigen_vector_gs,0,gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,0,low_idx));
        gsl_vector_set(sp_sys->sp_ion[i].eigen_vector_gs,1,gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,1,low_idx));
        gsl_vector_set(sp_sys->sp_ion[i].eigen_vector_gs,2,gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,2,low_idx));
        gsl_vector_set(sp_sys->sp_ion[i].eigen_vector_gs,3,gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,3,low_idx));
    }
    return;
}

/* Matrix Viewer: Print out a N x N matrix on console */
void sp_cluster_support_matrix_view( const gsl_matrix* m )
{
    if( m != NULL )
    {   for(int i=0;i<m->size1;i++)
        {   for(int j=0;j<m->size2;j++)
            {   //printf("%s%12.4lf",gsl_matrix_get(m,i,j)>0.?"+":"",gsl_matrix_get(m,i,j));
                printf("%12.4e",gsl_matrix_get(m,i,j));
            }
        puts("");
        }
    }
    else
        puts("sp_cluster_support_matrix_view input 'm' (gsl_matrix*) is a null pointer ... in SP_Support.c or SP_Support.h");

    return;
}

/* Vector Viewer: Print out a vector on console */
void sp_cluster_support_vector_view( const gsl_vector* v )
{
    if( v != NULL )
    {   for(int i=0;i<v->size;i++)
        {       //printf("%s%12.4f",gsl_vector_get(v,i)>0.?"+":"",gsl_vector_get(v,i));
                printf("%12.4e",gsl_vector_get(v,i));
        }
        puts("");
    }
    else
        puts("sp_cluster_support_vector_view input 'v' (gsl_vector*) is a null pointer ... in SP_Support.c or SP_Support.h");

    return;
}

/* Matrix Viewer: Print out a N x N matrix on file */
void sp_cluster_support_matrix_view_f( FILE* fp, const gsl_matrix* m )
{
    if( m != NULL )
    {   for(int i=0;i<m->size1;i++)
        {   for(int j=0;j<m->size2;j++)
                fprintf(fp,"%s%.18lf\t",gsl_matrix_get(m,i,j)>0.?"+":"",gsl_matrix_get(m,i,j));
        fprintf(fp,"\n");
        }
    }
    else
        fputs("sp_cluster_support_matrix_view input 'm' (gsl_matrix*) is a null pointer ... in SP_Support.c or SP_Support.h",fp);

    fflush(fp);

    return;
}

/* Vector Viewer: Print out a vector on file */
void sp_cluster_support_vector_view_f( FILE* fp, const gsl_vector* v )
{
    if( v != NULL )
    {   for(int i=0;i<v->size;i++)
            fprintf(fp,"%s%.18f\t",gsl_vector_get(v,i)>0.?"+":"",gsl_vector_get(v,i));
        fprintf(fp,"\n");
    }
    else
        fputs("sp_cluster_support_vector_view input 'v' (gsl_vector*) is a null pointer ... in SP_Support.c or SP_Support.h",fp);

    fflush(fp);

    return;
}

/* Transformation Matrix calculator */
/*
 * this function takes a vector, and the vector is transformed along transformed z-axis
 * the final return is the transformation matrix (rank 2 tensor), and the this data type is
 * gsl_matrix*
 */

gsl_matrix* sp_cluster_support_get_transformation_matrix( const gsl_vector* v )
{
    // gsl_vector_get(v,0) == 'x'
    // gsl_vector_get(v,1) == 'y'
    // gsl_vector_get(v,2) == 'z'
    gsl_matrix* pReturn = NULL;

    const double rxy = sqrt(pow(gsl_vector_get(v,0),2.)+pow(gsl_vector_get(v,1),2.));
    const double R   = sqrt(pow(gsl_vector_get(v,0),2.)+pow(gsl_vector_get(v,1),2.)+pow(gsl_vector_get(v,2),2.));
    double n1, n2, tmp; // dummy variables for workspace
    pReturn = gsl_matrix_calloc(4,4);   // Normally a rotation matrix is 3x3
                                        // For the 1st column or row has an element of '1' at T_11.

    if( pReturn != NULL )
    {
        gsl_matrix_set(pReturn,0,0,1.);

        if( gsl_vector_get(v,0) == 0. && gsl_vector_get(v,1) == 0. && gsl_vector_get(v,2) > 0. ) // if vector 'v' is on z-axis
        {   gsl_matrix_set(pReturn,1,1,1.); gsl_matrix_set(pReturn,2,2,1.); gsl_matrix_set(pReturn,3,3,1.);
            // set the matrix as I
        }
        else if( gsl_vector_get(v,0) == 0. && gsl_vector_get(v,1) == 0. && gsl_vector_get(v,2) < 0. )   // if vector 'v' is on negative z-axis
        {   gsl_matrix_set(pReturn,1,1,1.); gsl_matrix_set(pReturn,2,2,1.); gsl_matrix_set(pReturn,3,3,-1.);
            // set the matrix has xy-plane reflection
        }
        else
        {   gsl_matrix_set(pReturn,3,1,gsl_vector_get(v,0)/R);
            gsl_matrix_set(pReturn,3,2,gsl_vector_get(v,1)/R);
            gsl_matrix_set(pReturn,3,3,gsl_vector_get(v,2)/R);  // set k' in the transformed (local) symmetry

            gsl_matrix_set(pReturn,2,1,gsl_vector_get(v,2)*gsl_vector_get(v,0)/rxy);
            gsl_matrix_set(pReturn,2,2,gsl_vector_get(v,2)*gsl_vector_get(v,1)/rxy);
            gsl_matrix_set(pReturn,2,3,-R*sqrt(1.-gsl_vector_get(v,2)*gsl_vector_get(v,2)/R/R));

            n1 = 1./sqrt(pow(gsl_matrix_get(pReturn,2,1),2.)
               + pow(gsl_matrix_get(pReturn,2,2),2.)
               + pow(gsl_matrix_get(pReturn,2,3),2.));

            for(int i=0;i<3;i++)
            {   tmp = gsl_matrix_get(pReturn,2,i+1);
                tmp = tmp*n1;
                gsl_matrix_set(pReturn,2,i+1,tmp);  // set j' in the transformed (local) symmetry
            }

            gsl_matrix_set(pReturn,1,1,n1/R*(pow(gsl_vector_get(v,2),2.)*gsl_vector_get(v,1)/rxy
                        + R*gsl_vector_get(v,1)*sqrt(1.-pow(gsl_vector_get(v,2),2.)/R/R)));
            gsl_matrix_set(pReturn,1,2,n1/R*(-pow(gsl_vector_get(v,2),2.)*gsl_vector_get(v,0)/rxy
                        - R*gsl_vector_get(v,0)*sqrt(1.-pow(gsl_vector_get(v,2),2.)/R/R)));
            gsl_matrix_set(pReturn,1,3,0.);

            n2 = 1./sqrt(pow(gsl_matrix_get(pReturn,1,1),2.)
                    + pow(gsl_matrix_get(pReturn,1,2),2.)
                    + pow(gsl_matrix_get(pReturn,1,3),2.));

            for(int i=0;i<3;i++)
            {   tmp = gsl_matrix_get(pReturn,1,i+1);
                tmp = tmp/n2;
                gsl_matrix_set(pReturn,1,i+1,tmp);  // set i' in the transformed (local) symmetry
            }
        }
    }
    else
        puts("sp_cluster_get_trans_mat 'pReturn' alloc error ... in SP_support.c or SP_support.h");

    return pReturn;
}   /* SEE THE DETAILS IN THE 'N1' RING BINDER */

double sp_cluster_support_kronecker_delta( int a, int b )
{
    double Return = 0.;

    if( a == b )
        Return = 1.;

    return Return;
}

// norm of vector
double sp_cluster_support_get_norm( double v1, double v2, double v3 )
{   return pow(v1*v1+v2*v2+v3*v3,0.5);
}

// SPLINER

double** sp_cluster_support_get_spline( const double** data, const int knot_stride )
{
    double** pReturn = (double**)malloc((knot_stride-1)*sizeof(double*));
    for(int i=0;i<knot_stride-1;i++)
        pReturn[i] = (double*)calloc(4,sizeof(double));

    // Work Space
    const int n = knot_stride;
    double f[n];  double t[n];
    double a[n];  double b[n];  double c[n];  double d[n];
    double z[n];  double L[n];  double h[n];
    double alpha[n];  double mu[n];

    // Measuring slopes
    const double rslope = 0.;
    const double lslope = (data[3][1]-data[0][1])/(data[3][0]-data[0][0]);

    // Memset
    memset(f,0,n*sizeof(double));   memset(t,0,n*sizeof(double));
    memset(a,0,n*sizeof(double));   memset(b,0,n*sizeof(double));   memset(c,0,n*sizeof(double));   memset(d,0,n*sizeof(double));
    memset(z,0,n*sizeof(double));   memset(L,0,n*sizeof(double));   memset(h,0,n*sizeof(double));
    memset(alpha,0,n*sizeof(double));   memset(mu,0,n*sizeof(double));

    // Make Spline

    for(int i=0;i<n;i++)
    {   f[i] = data[i][1];
        t[i] = data[i][0];      }

    h[0] = t[1] - t[0];
    alpha[0] = 3.*(f[1]-f[0])/h[0]-3.*lslope;
    L[0] = 2.*h[0];
    mu[0] = 0.5;
    z[0] = alpha[0]/L[0];
    b[0] = lslope;
    // End Left Init

    for(int k=1;k<n-1;k++)
    {
        h[k] = t[k+1] - t[k];

        alpha[k] = (3./(h[k]*h[k-1]))*(f[k+1]*h[k-1]-f[k]*(h[k]+h[k-1])+f[k-1]*h[k]);

        L[k] = 2.*(h[k]+h[k-1])-h[k-1]*mu[k-1];

        mu[k] = h[k]/L[k];

        z[k] = (alpha[k]-h[k-1]*z[k-1])/L[k];
    }

    // Right End init
    alpha[n-1] = 3.*rslope - 3.*(f[n-1]-f[n-2])/h[n-2];
    L[n-1] = h[n-2]*(2.-mu[n-2]);
    z[n-1] = (alpha[n-1]-h[n-2]*z[n-2])/L[n-1];
    c[n-1] = z[n-1];

    for(int j=n-2;j>=0;j--)
    {
        c[j] = z[j]-mu[j]*c[j+1];

        b[j] = (f[j+1]-f[j])/h[j] - h[j]*(c[j+1]+2.*c[j])/3.;

        d[j] = (c[j+1]-c[j])/(3.*h[j]);

        a[j] = f[j];
    }
    // R-E END

    // SAVE DATA
    for(int i=0;i<knot_stride-1;i++)
    {   pReturn[i][0] = d[i];
        pReturn[i][1] = c[i] - 3.*t[i]*d[i];
        pReturn[i][2] = b[i] - 2.*c[i]*t[i] + 3.*t[i]*t[i]*d[i];
        pReturn[i][3] = a[i] - b[i]*t[i] + c[i]*t[i]*t[i] - d[i]*t[i]*t[i]*t[i];
    }

    return pReturn;
}

//int sp_cluster_support_get_atom_number(

void sp_cluster_support_print_xyz( void* sp_sys_void, const double cur_energy, const int rank, const int numtasks )
{
	int offset;
	// cast the type into "sp_cluster_system*"
	sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;

	int atom_number=0;
	for(int i=0;i<sp_sys->number_of_classic_ion;i++)
	{	if( sp_sys->classic_ion[i].if_shell == SP_SUPPORT_FALSE )
			atom_number++;
		// CNT ONLY WHEN ITS CORE
	}

	if( cur_energy != 0. )
	{
		//printf("\t%d\n",sp_sys->number_of_classic_ion+sp_sys->number_of_sp_ion);
		printf("\t%d\n",atom_number+sp_sys->number_of_sp_ion);
		printf(" SCF DONE %.6lf\n",cur_energy);
	}

	// CORE POSITION ONLY ... "xyz" Format Compatible
	for(int n=0;n<sp_sys->number_of_sp_ion+sp_sys->number_of_classic_ion;n++)
	{
    	    offset = n-sp_sys->number_of_classic_ion;
	    if( n < sp_sys->number_of_classic_ion )
	    {
		if( sp_sys->classic_ion[n].if_shell == SP_SUPPORT_FALSE ) // if it is core
		{	fprintf(stdout,"%3s%12.6lf%12.6lf%12.6lf\n", sp_sys->classic_ion[n].atom_name,
			gsl_vector_get(sp_sys->classic_ion[n].core_position,0), gsl_vector_get(sp_sys->classic_ion[n].core_position,1), gsl_vector_get(sp_sys->classic_ion[n].core_position,2));
		}
	    }
	    else // this is for printing sp-ions
	    {
		    offset = n - sp_sys->number_of_classic_ion;
		    fprintf(stdout,"%3s%12.6lf%12.6lf%12.6lf\n", sp_sys->sp_ion[offset].atom_name,
			    gsl_vector_get(sp_sys->sp_ion[offset].core_position,0), gsl_vector_get(sp_sys->sp_ion[offset].core_position,1), gsl_vector_get(sp_sys->sp_ion[offset].core_position,2));
	    }
	}

	printf("\n");
	printf(" CONFIGURATION_XYZ_SC_INFO\n");
	printf(" %d\t%d\n",sp_sys->number_of_classic_ion,sp_sys->number_of_sp_ion);

	// SHOW SHELL CORE POSITION BOTH
	for(int n=0;n<sp_sys->number_of_sp_ion+sp_sys->number_of_classic_ion;n++)
	{
    	    offset = n-sp_sys->number_of_classic_ion;
	    if( n < sp_sys->number_of_classic_ion )
	    {
		if( sp_sys->classic_ion[n].if_shell == SP_SUPPORT_FALSE )
		{
			fprintf(stdout,"%3s%3s%12.6lf%12.6lf%12.6lf\n", sp_sys->classic_ion[n].atom_name,"c",
			gsl_vector_get(sp_sys->classic_ion[n].core_position,0), gsl_vector_get(sp_sys->classic_ion[n].core_position,1), gsl_vector_get(sp_sys->classic_ion[n].core_position,2));
		}
		else if( sp_sys->classic_ion[n].if_shell == SP_SUPPORT_TRUE )
		{
			fprintf(stdout,"%3s%3s%12.6lf%12.6lf%12.6lf\n", sp_sys->classic_ion[n].atom_name,"s",
			gsl_vector_get(sp_sys->classic_ion[n].core_position,0), gsl_vector_get(sp_sys->classic_ion[n].core_position,1), gsl_vector_get(sp_sys->classic_ion[n].core_position,2));
		}
	    }
	    else // this is for printing sp-ions
	    {
		    offset = n - sp_sys->number_of_classic_ion;
		    fprintf(stdout,"%3s%15.6lf%12.6lf%12.6lf\n", sp_sys->sp_ion[offset].atom_name,
			    gsl_vector_get(sp_sys->sp_ion[offset].core_position,0), gsl_vector_get(sp_sys->sp_ion[offset].core_position,1), gsl_vector_get(sp_sys->sp_ion[offset].core_position,2));
	    }
	}

	return;
}

/// DIIS_SUPPORT TOOLS

/*
int if_diis;
int diis_max_depth;
int diis_cur_depth;

gsl_vector** diis_error_vector;                         // SAVE AS CIRCULAR QUEUE
gsl_vector*  diis_coefficient_vector;
gsl_permutation* diis_ws_p;
gsl_matrix*  diis_error_matrix;
gsl_matrix*  diis_error_matrix_inv;
*/

void sp_cluster_support_diis_error_vector_queue_init( void* sp_sys_void )
{
    sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    sp_sys->diis_error_vector_queue_front = 0;
    sp_sys->diis_error_vector_queue_rear  = 0;
    return;
}

int sp_cluster_support_diis_error_vector_queue_isfull( void* sp_sys_void )
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    if( (sp_sys->diis_error_vector_queue_rear+1)%sp_sys->diis_max_depth == sp_sys->diis_error_vector_queue_front )
        return SP_SUPPORT_TRUE;
    else
        return SP_SUPPORT_FALSE;
}

int sp_cluster_support_diis_error_vector_queue_isempty( void* sp_sys_void )
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    if( sp_sys->diis_error_vector_queue_front == sp_sys->diis_error_vector_queue_rear )
        return SP_SUPPORT_TRUE;
    else
        return SP_SUPPORT_FALSE;
}

void sp_cluster_support_diis_error_vector_enqueue( void* sp_sys_void )
{
    sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    double delta_evec[4];
    int low_idx;

    double delta_evec_sum = 0.;

    //if( !((sp_sys->diis_error_vector_queue_rear+1)%sp_sys->diis_max_depth == sp_sys->diis_error_vector_queue_front) )   // check if queue is full
    if( sp_cluster_support_diis_error_vector_queue_isfull( sp_sys ) == SP_SUPPORT_FALSE )
    {
        sp_sys->diis_error_vector_queue_rear = (sp_sys->diis_error_vector_queue_rear+1)%sp_sys->diis_max_depth;         // rear index ++
        sp_sys->diis_cur_depth++;

        for(int i=0;i<sp_sys->number_of_sp_ion;i++)
        {
            low_idx = sp_cluster_support_get_lowest_state(sp_sys->sp_ion[i].eigen_value);

            // current eigenvector - (loaded) previous eigenvector_gs
            delta_evec[0] = gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,0,low_idx) - gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,0);
            delta_evec[1] = gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,1,low_idx) - gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,1);
            delta_evec[2] = gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,2,low_idx) - gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,2);
            delta_evec[3] = gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,3,low_idx) - gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,3);

            delta_evec_sum += delta_evec[0]*delta_evec[0] + delta_evec[1]*delta_evec[1] + delta_evec[2]*delta_evec[2] + delta_evec[3]*delta_evec[3];

            gsl_vector_set(sp_sys->diis_error_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+0,delta_evec[0]);
            gsl_vector_set(sp_sys->diis_error_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+1,delta_evec[0]);
            gsl_vector_set(sp_sys->diis_error_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+2,delta_evec[0]);
            gsl_vector_set(sp_sys->diis_error_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+3,delta_evec[0]);

            // loading .. previous eigen_vector_gs (note that the function 'sp_cluster_support_load_gs_eigenvector()' must be called before enqueue)
            gsl_vector_set(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+0, gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,0));
            gsl_vector_set(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+1, gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,1));
            gsl_vector_set(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+2, gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,2));
            gsl_vector_set(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ],i*4+3, gsl_vector_get(sp_sys->sp_ion[i].eigen_vector_gs,3));

            #ifdef DEBUG_SUPPORT
            printf("\ndelta evec\n");
            printf("%12.6lf%12.6lf%12.6lf%12.6lf\n",delta_evec[0],delta_evec[1],delta_evec[2],delta_evec[3]);
            printf("%12.6lf%12.6lf%12.6lf%12.6lf\n",
            gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,0,low_idx),
            gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,1,low_idx),
            gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,2,low_idx),
            gsl_matrix_get(sp_sys->sp_ion[i].eigen_vector,3,low_idx));
            #endif
/*
            printf("%12.6lf%12.6lf%12.6lf%12.6lf\n",gsl_vector_get(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ], i*4 + 0 ),
            gsl_vector_get(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ], i*4 + 1 ),
            gsl_vector_get(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ], i*4 + 2 ),
            gsl_vector_get(sp_sys->diis_prev_eigen_vector[ sp_sys->diis_error_vector_queue_rear ], i*4 + 3 ));
*/
        /* Note that the both,

            Error vector ( sp_sys->diis_error_vector )
            Prev  evecs  ( sp_sys->diis_prev_eigen_vector )

            Follows the form of (gsl_vector)v[i][j]

            v[i]    ... 'i' -> index in the queue // queue is circular type, max dept is 'sp_sys->diis_max_depth'
            v[i][j] ... 'j' -> DATA, s_sp1, px_sp1, py_sp1, pz_sp1 // s_sp2, px_sp2, py_sp2, pz_sp2 // ....
                                i.e., length of number_of_sp_ions * 4 (each lone pair ground-state MO saved with stride of 4)
        */
        }
        printf("delta_evec_sum: %12.8lf\n",sqrt(delta_evec_sum)/sp_sys->number_of_sp_ion);

    }
    return;     // if queue is full do nothing
}

void sp_cluster_support_diis_error_vector_dequeue( void* sp_sys_void )
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    if( sp_cluster_support_diis_error_vector_queue_isempty( sp_sys ) == SP_SUPPORT_FALSE )
    {   sp_sys->diis_error_vector_queue_front = (sp_sys->diis_error_vector_queue_front+1)%sp_sys->diis_max_depth;
        sp_sys->diis_cur_depth--;
    }
    #ifdef DEBUG_SUPPORT
    const int front = sp_sys->diis_error_vector_queue_front;
    const int rear  = sp_sys->diis_error_vector_queue_rear;
    printf("%d\t%d\n",front,rear);
    #endif

    return;     // if queue is empty do nothing
}

double sp_cluster_support_diis_get_error_vector_gnorm( void* sp_sys_void )      // returns latest .. i.e., the one at rear
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    double ret;
    if( sp_cluster_support_diis_error_vector_queue_isempty( sp_sys ) == SP_SUPPORT_FALSE )
    {   ret = gsl_blas_dnrm2( sp_sys->diis_error_vector[ sp_sys->diis_error_vector_queue_rear ] );
        return ret/sp_sys->number_of_sp_ion;
    }
    return 0.;
}

void sp_cluster_support_diis_solve_least_square_problem( void* sp_sys_void )
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;

    // int count = front > rear ? (MAX - front + rear) : (rear - front);
    const int queue_capacity = sp_sys->diis_max_depth;
    const int front = sp_sys->diis_error_vector_queue_front;
    const int rear  = sp_sys->diis_error_vector_queue_rear;
    const int cur_size = front > rear ? (queue_capacity - front + rear) : (rear - front);       // number of elements in the queue (error_vector)

/*
gsl_vector*  diis_coefficient_vector;
gsl_permutation* diis_ws_p;
gsl_matrix*  diis_error_matrix;
gsl_matrix*  diis_error_matrix_inv;
*/
    size_t len = cur_size + 1;

printf("\n#### IN LSP SOLVER / len : %zu\n",len);
printf("ws size : ErrorMat / %zu %zu\n",sp_sys->diis_error_matrix->size1,sp_sys->diis_error_matrix->size2);
printf("cursize/max_depth/cur_depth : %4d%4d%4d\n", cur_size, sp_sys->diis_max_depth - 1, sp_sys->diis_cur_depth );
    //if( !(cur_size == (sp_sys->diis_max_depth - 1)) )
    //if( !(cur_size == (sp_sys->diis_max_depth)) )
    //if( !(len == (sp_sys->diis_max_depth)) )
    if( !(len == sp_sys->diis_error_matrix->size1) )
    {
printf("Resizing!!\n");
        // if current size is not same with actual max_depth in use ... i.e., max_depth - 1, since circular queue in use
        // have to resize the workspaces ... diis_ws_p (gsl_permutation*) / diis_error_matrix (gsl_matrix*) / diis_error_matrix_inv (gsl_matrix*)

        gsl_vector_free(sp_sys->diis_coefficient_vector);
        gsl_vector_free(sp_sys->diis_least_square_condition);
        gsl_permutation_free(sp_sys->diis_ws_p);
        gsl_matrix_free(sp_sys->diis_error_matrix);
        gsl_matrix_free(sp_sys->diis_error_matrix_inv);

        sp_sys->diis_coefficient_vector = gsl_vector_calloc(len);
        sp_sys->diis_least_square_condition = gsl_vector_calloc(len);
        sp_sys->diis_ws_p = gsl_permutation_calloc(len);
        sp_sys->diis_error_matrix = gsl_matrix_calloc(len,len);
        sp_sys->diis_error_matrix_inv = gsl_matrix_calloc(len,len);      // len = cur_size (error_vector length) + 1 // 1 is to represent least square sense of coeff
    }

    int ii=0; int jj=0;     // tags to get actual matrix indices !
    double elem_tmp;
    for(int i=(front+1)%queue_capacity; i!=(rear+1)%queue_capacity ; i=(i+1)%queue_capacity)    // for loop circulate over queue elements
    {
        for(int j=(front+1)%queue_capacity; j!=(rear+1)%queue_capacity ; j=(j+1)%queue_capacity)
        {
            #ifdef DEBUG_SUPPORT
            printf("i,j / tag_i,tag_j : %d\t%d / %d\t%d\n",i,j,ii,jj);
            //printf("size1 / size2 : %d\t%d\n", sp_sys->diis_error_matrix->size1, sp_sys->diis_error_matrix->size2);
            #endif

            gsl_blas_ddot(sp_sys->diis_error_vector[j],sp_sys->diis_error_vector[i],&elem_tmp);
            gsl_matrix_set(sp_sys->diis_error_matrix,ii,jj,elem_tmp);
            jj++;
        }
        jj=0;
        ii++;
    }
    printf("ERROR - 0 ?\n");
    ///// SET RESTS
     for(int i=0;i<len;i++)
    {   gsl_matrix_set(sp_sys->diis_error_matrix,i,len-1,-1.);
        gsl_matrix_set(sp_sys->diis_error_matrix,len-1,i,-1.);        }
    printf("ERROR?\n");
    gsl_matrix_set(sp_sys->diis_error_matrix,len-1,len-1,0.);
    printf("ERROR?\n");

    ///// SET least Square Condition vector
    gsl_vector_set(sp_sys->diis_least_square_condition,len-1,-1.);

    //// calculate inverse
    int signum; // variable for LU_decomp

    #ifdef DEBUG_SUPPORT
    printf("ErrorMatrix\n");
    sp_cluster_support_matrix_view( sp_sys->diis_error_matrix );
    #endif

    /*  Say,
        Error matrix                    E
        Coefficient Vector              C
        RHS Least Square Condition      L

        Need to find 'C' vector

        Solve : EC = L, i.e., C = E_Inverse  L
    */

    gsl_linalg_LU_decomp(sp_sys->diis_error_matrix,sp_sys->diis_ws_p,&signum);
    gsl_linalg_LU_invert(sp_sys->diis_error_matrix,sp_sys->diis_ws_p,sp_sys->diis_error_matrix_inv);    // invert is saved in 'sp_sys->diis_error_matrix_inv'
    // Inverse Found

    // int gsl_blas_dgemv(CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix *A, const gsl_vector *x, double beta, gsl_vector *y)
    gsl_blas_dgemv(CblasNoTrans,1.,sp_sys->diis_error_matrix_inv,sp_sys->diis_least_square_condition,0.,sp_sys->diis_coefficient_vector);
    // Solve LeaseSquare Problem Answer saved in 'sp_sys->diis_coefficient_vector'

    #ifdef DEBUG_SUPPORT
    //sp_cluster_support_vector_view( sp_sys->diis_least_square_condition );
    printf("coefficient vector, ignore the last dummy lambda value\n");
    sp_cluster_support_vector_view( sp_sys->diis_coefficient_vector );
    //sp_cluster_support_matrix_view( sp_sys->diis_error_matrix );
    printf("InverseMatrix\n");
    sp_cluster_support_matrix_view( sp_sys->diis_error_matrix_inv );
    printf("cursize/max_depth/cur_depth : %d\t%d\t%d\n", cur_size, sp_sys->diis_max_depth - 1, sp_sys->diis_cur_depth );

    printf(" ---- FINALISE LEAST SQUARE SOLVER \n\n");
    #endif
    return;
}

void sp_cluster_support_diis_least_square_result_update( void* sp_sys_void )
{   sp_cluster_system* sp_sys = (sp_cluster_system*)sp_sys_void;
    int low_idx;
    int jj = 0;
    double cs, cx, cy, cz, norm_factor;

    const int queue_capacity = sp_sys->diis_max_depth;
    const int front = sp_sys->diis_error_vector_queue_front;
    const int rear  = sp_sys->diis_error_vector_queue_rear;

    for(int i=0;i<sp_sys->number_of_sp_ion;i++)
    {
        low_idx = sp_cluster_support_get_lowest_state(sp_sys->sp_ion[i].eigen_value);

        // results are in 'sp_sys->diis_coefficient_vector' ... data structure sp_sys->diis_coefficient_vector->size or 'stride'?
        // size can also be obtained by 'sp_sys->diis_cur_depth'

        jj = 0;
        cs = 0.; cx = 0.; cy = 0.; cz = 0.;

        for(int j=(front+1)%queue_capacity; j!=(rear+1)%queue_capacity ; j=(j+1)%queue_capacity)
        {
            cs += gsl_vector_get( sp_sys->diis_coefficient_vector,jj) * gsl_vector_get(sp_sys->diis_prev_eigen_vector[j],i*4+0);
            cx += gsl_vector_get( sp_sys->diis_coefficient_vector,jj) * gsl_vector_get(sp_sys->diis_prev_eigen_vector[j],i*4+1);
            cy += gsl_vector_get( sp_sys->diis_coefficient_vector,jj) * gsl_vector_get(sp_sys->diis_prev_eigen_vector[j],i*4+2);
            cz += gsl_vector_get( sp_sys->diis_coefficient_vector,jj) * gsl_vector_get(sp_sys->diis_prev_eigen_vector[j],i*4+3);
            jj++;

            printf("index j / tag jj: %d\t%d\n",j,jj);
        }

#ifdef DEBUG_SUPPORT
printf("%dth\t%12.4lf%12.4lf%12.4lf%12.4lf%20.4lf\n",i+1,cs,cx,cy,cz,sqrt(cs*cs+cx*cx+cy*cy+cz*cz));
#endif
        norm_factor = sqrt(cs*cs+cx*cx+cy*cy+cz*cz);
        norm_factor = 1.;

        gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,0,low_idx,cs/norm_factor);
        gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,1,low_idx,cx/norm_factor);
        gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,2,low_idx,cy/norm_factor);
        gsl_matrix_set(sp_sys->sp_ion[i].eigen_vector,3,low_idx,cz/norm_factor);
        // set eigenvector ... updated!!!

        // finish up the rest ...
    }

    return;
}

ext: c


text:
/** @file header.h
 *
 */

#ifndef HEADER_H
#define HEADER_H

#define _GNU_SOURCE

#include <time.h>
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <omp.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_sf_expint.h>
#include <ctype.h>
#include "cuba.h"

#define PSC  		101L
#define ST 		    	102L
#define TK 		    	103L

#define GROWTH		104L
#define DERGROWTH		105L

#define NONLINEAR	   	106L
#define LINEAR		107L

#define GAUSSIAN		114L
#define NONGAUSSIAN	115L

#define INIT  		116L
#define LOCAL		117L
#define EQUILATERAL	118L
#define ORTHOGONAL	119L
#define QSF			120L
#define HS			121L
#define NGLOOP		122L
#define derNGLOOP		123L

#define QUADRATIC		124L
#define TIDE		125L
#define GAMMA		126L

#define LPOWER		127L
#define NLPOWER		128L
#define CDM 		90L
#define BA 			91L
#define TOT 		92L

#define CO10  		131L
#define CO21  		132L
#define CO32  		133L
#define CO43  		134L
#define CO54  		135L
#define CO65  		136L
#define CII 		137L

#define MATTER		138L
#define LINEMATTER	139L
#define LINE 		140L

#define DST			141L
#define GFILTER		142L
#define BSPLINE		143L

#define TREE		144L
#define LOOP		145L
#define WIR			146L
#define NOIR		147L

#define HALO 		148L

#define MEAN            149L
#define HMshot          150L
#define HMclust		151L

#define PS_KMIN   	1.e-5
#define PS_KMAX   	300.1
#define PS_ZMAX   	14.

#define CLEANUP        1
#define DO_NOT_EVALUATE -1.0
#define NPARS     6
#define MAXL 2000

extern struct globals gb;

/**
 * List of limHaloPT header files
 */
/// \cond DO_NOT_DOCUMENT
#include "../Class/include/class.h"
/// \endcond
#include "global_structs.h"
#include "hcubature.h"
#include "utilities.h"
#include "setup_teardown.h"
#include "cosmology.h"
#include "ir_res.h"
#include "line_ingredients.h"
#include "ps_halo_1loop.h"
#include "ps_line_hm.h"
#include "ps_line_pt.h"
#include "survey_specs.h"
#include "wnw_split.h"

/**
 * A structure passed to the integrators to hold the parameters fixed in the integration
 */
struct integrand_parameters
{
	double p1;
	double p2;
	double p3;
	double p4;
	double p5;
	double p6;
	double p7;
	double p8;
	double p9;
	double p10;
	double p11;
	long 	p12;
	long 	p13;
};

/**
 * Another structure passed to the integrators to hold the parameters fixed in the integration
 */
struct integrand_parameters2
{

	struct Cosmology *p1;
	struct Cosmology *p2;
	struct Cosmology *p3;

	double p4;
	double p5;
	double p6;
	double p7;
	double p8;
	double p9;
	double p10;
	double p11;
	double p12;

	long p13;
	long p14;
	long p15;
	long p16;
	long p17;
	long p18;

	int p19;

	double *p20;
	size_t p22;
};

#endif

ext: h

max_stars_repo_name: amoradinejad/limHaloPT

max_issues_repo_name: amoradinejad/limHaloPT

max_forks_repo_name: amoradinejad/limHaloPT


text:
#include "cones.h"
#ifdef LAPACK_LIB_FOUND
#include <cblas.h>
#include <lapacke.h>
#endif
void projectsdc(double * X, int n, Work * w);

/* in place projection (with branches) */
void projCone(double *x, Cone * k, Work * w)
{
	int i;
	int count;

	/* project onto positive orthant */
	for(i = k->f; i < k->f+k->l; ++i)
	{
		if(x[i] < 0.0) x[i] = 0.0;
		//x[i] = (x[i] < 0.0) ? 0.0 : x[i];
	}
	count = k->l+k->f;
	/* project onto SOC */
	for(i = 0; i < k->qsize; ++i)
	{
		double v1 = x[count];
		double s = calcNorm(&(x[count+1]),k->q[i]-1);
		double alpha = (s + v1)/2.0;

		if(s <= v1) { /* do nothing */ }
		else if (s <= - v1) {
			memset(&(x[count]), 0, k->q[i]*sizeof(double));
		} else {
			x[count] = alpha;
			scaleArray(&(x[count+1]), alpha/s, k->q[i]-1);
      //cblas_dscal(k->q[i]-1, alpha/s, &(x[count+1]),1);
		}
		count += k->q[i];
	}
#ifdef LAPACK_LIB_FOUND
	/* project onto PSD cone */
	for (i=0; i < k->ssize; ++i){
		projectsdc(&(x[count]),k->s[i],w);
		count += (k->s[i])*(k->s[i]);
	}
#else
  if(k->ssize > 0){
    coneOS_printf("WARNING: solving SDP, no lapack library specified in makefile!\n");
    coneOS_printf("ConeOS will return a wrong answer!\n");
}
#endif
	/* project onto OTHER cones */
}

#ifdef LAPACK_LIB_FOUND
void projectsdc(double *X, int n, Work * w)
{ /* project onto the positive semi-definite cone */
  if (n == 1) {
    if(X[0] < 0.0) X[0] = 0.0;
    return;
  }

  int i, j, m=0;
  double * Xs = w->Xs;
  double * Z = w->Z;
  double * e = w->e;
  memcpy(Xs,X,n*n*sizeof(double));

  // Xs = X + X', save div by 2 for eigen-recomp
  for (i = 0; i < n; ++i){
    cblas_daxpy(n, 1, &(X[i]), n, &(Xs[i*n]), 1);
    //b_daxpy(n, 1, &(X[i]), n, &(Xs[i*n]), 1);
  }

  double EIG_TOL = 1e-8;
  double vupper = calcNorm(Xs,n*n);
  LAPACKE_dsyevr( LAPACK_COL_MAJOR, 'V', 'V', 'U', n, Xs, n, 0.0, vupper, -1, -1, EIG_TOL, &m, e, Z, n , NULL);

  memset(X, 0, n*n*sizeof(double));
  for (i = 0; i < m; ++i) {
    cblas_dsyr(CblasColMajor, CblasLower, n, e[i]/2, &(Z[i*n]), 1, X, n);
    //b_dsyr('L', n, -e[i]/2, &(Z[i*n]), 1, Xs, n);
  }
  // fill in upper half
  for (i = 0; i < n; ++i){
    for (j = i+1; j < n; ++j){
      X[i + j*n] = X[j + i*n];
    }
  }
}
#endif

max_stars_repo_name: cvxgrp/coneos

max_issues_repo_name: cvxgrp/coneos

max_forks_repo_name: cvxgrp/coneos