text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <omp.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <sys/time.h>

// ***************************************************************************************
// * Pre-Port Entwicklung und Testen eines OMP-parallisierten, adaptiven Simpson-integrators für
// * die fiesen doppelintgrale in imd_colrad.c
// ***************************************************************************************
int num_threads;
int simpson_error;
int funcevals;
int iter;

const double tolmax=1e-30;
#define INITIAL_STACK_SIZE 128; //128   /* initial size of new stacks */

const double alpha=0.816496580927726;
const double beta=0.447213595499958;
static const double xgkq[12] =
{
  0.0,
  -0.942882415695480,
  -0.816496580927726,
  -0.641853342345781,
  -0.447213595499958,
  -0.236383199662150,
  0.0,
  0.236383199662150,
  0.447213595499958,
  0.641853342345781,
  0.816496580927726,
  0.942882415695480
};

/* the stack structure */
struct stack_s{
  int el_count;            /* count of elements on stack */
  int el_size;             /* size of an element */
  int mem_reserve;         /* allocated memory for stack */
  void* elements;          /* pointer to begin of stack */
};

struct my_f_params { double ne; double T;double mu; double E;double DeltaE; int allowed;};

typedef struct _work_t{
  double a;
  double b;
  double tol;
  double S;
  double fa;
  double fb;
  double fm;
  double rec;
  int iter;
  struct my_f_params * p; //pointer auf params
} work_t;

typedef struct stack_s* stack_t;

typedef struct _work_t_gkq{
  double a;
  double b;
  double toler;
  double I_13;
  double fa;
  double fb;
  struct my_f_params * p; //pointer auf params
  double (*f)(double, struct my_f_params*);
  stack_t stack_inner;
  short int is_parent;
  short int subtask_nr; //inner task nr
  short int task_nr; //outer task nr
  short int subtasks_left; //erst pop'n wenn aller inner tasks (intergals) vollständig!
  double inner_integrals[5];
} work_gkq;

// **************************************** FUNCTION DEFS *************************************************
void create_stack(stack_t* stack, int element_size);
int empty_stack(stack_t stack);
void push_stack(stack_t stack, void* element);
void pop_stack(stack_t stack, void* element);

double
integral2(     double (*f)(double, struct my_f_params*), /* function to integrate */
     stack_t stack);

double gkq_adapt_double(stack_t stack);
double gkq_adapt_single(stack_t stack); //for single integral, first integral

double gkq_double(double (*f)(double, struct my_f_params*), double (*finner)(double, struct my_f_params*),
                  double a, double b,
                  double TOL, struct my_f_params* p,stack_t stack);

double gkq_adapt_serial(double (*f)(double, struct my_f_params*), double a, double b, double fa,double fb, double toler,double I_13, struct my_f_params* p);

work_gkq gkq_create_inner_task(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p);
double gkq_single(double (*f)(double, struct my_f_params*), double a, double b,
                       double TOL, struct my_f_params* p,stack_t stack);

int terminate_serial;
int terminate_gkq;

// static double myfun(double x,struct my_f_params* p)
// {
//   double T=p->T;
//   double fermi=1/(1+exp(-x/T));
//   double fermi2=1- 1/(1+exp(-x/T));
//   double sigma=x/T*log(x/T*2)/pow(1/x,2.0); //einfach nur ärger machen

//   // return  exp(-x*x)/p->T+log(x*pow(p->T,2.0))/p->ne;
//   return fermi*fermi2*sigma;
// }

// static double myfun2(double x,void* pv)
// {
//   struct my_f_params *p= (struct my_f_params*) pv;

//   double T=p->T;
//   double fermi=1/(1+exp(-x/T));
//   double fermi2=1- 1/(1+exp(-x/T));
//   double sigma=x/T*log(x/T*2)/pow(1/x,2.0); //einfach nur ärger machen

//   // return  exp(-x*x)/p->T+log(x*pow(p->T,2.0))/p->ne;
//   return fermi*fermi2*sigma;
// }

static double inner_integrand(double x,void *pv)
{
  struct my_f_params *p= (struct my_f_params*) pv;
  double T=p->T;
  double fermi=1/(1+exp(-x/T));
  double fermi2=1- 1/(1+exp(-x/T));
  return exp(-x*x); //fermi*fermi2;
}
static double inner_integrand2(double x,struct my_f_params* p)
{
  double T=p->T;
  double fermi=1/(1+exp(-x/T));
  double fermi2=1- 1/(1+exp(-x/T));
  return fermi*fermi2;
}
void get_integ_bounds_inner(double *integ_bnd, double x, struct my_f_params *p)
{
  integ_bnd[0]=-4;
  integ_bnd[1]=4;
}

gsl_integration_workspace * gsinner;
gsl_integration_workspace * gsinner2;

static double myfun(double x,struct my_f_params* p)
{
  double T=p->T;
  double fermi=1/(1+exp(-x/T));
  double fermi2=1- 1/(1+exp(-x/T));
  double sigma=x/T*log(x/T*2)/pow(1/x,2.0); //einfach nur ärger machen

  // double result, err;
  // gsl_function gslfun;
  // gslfun.function=&inner_integrand;
  // gslfun.params=p;

  // gsl_integration_qag(&gslfun, x, 300, 0.0, 1e-3, 1000,1,
  //                      gsinner2, &result, &err);

  // gsl_integration_workspace_free(gsinner2);
  // stack_t stack2;
  // create_stack(&stack2, sizeof(work_gkq));

  //double result=gkq(inner_integrand2, x, 600, 1e-3, p);//, stack);
  // free(stack2->elements);

  return sigma;
  //return  exp(-x*x);
};

static double myfun_gsl(double x,void* pv)
{
  struct my_f_params *p= (struct my_f_params*) pv;

  double T=p->T;
  double fermi=1/(1+exp(-x/T));
  double fermi2=1- 1/(1+exp(-x/T));
  double sigma=x/T*log(x/T*2)/pow(1/x,2.0); //einfach nur ärger machen

  gsl_function gslfun;
  gslfun.function=&inner_integrand;
  gslfun.params=pv;
  double result, err;
  gsl_integration_qag(&gslfun, x, 600, 0.0, 1e-3, 1000,1,
                       gsinner, &result, &err);

  return sigma*result;
}

// static double myfun2(double x,void* pv)
// {
//   struct my_f_params *p= (struct my_f_params*) pv;
// static float sinfc(float x) { return sinf(x); }
// static double frand48c(double x) { funcevals++; return (double) drand48(); }

// *************************************************************************************************************
int main(int argc,char** argv)
{
  iter=0;
  num_threads=omp_get_num_threads();
  funcevals=0;
  simpson_error=0;

  struct timeval start, end;

  double gsltime, simpsontime;
  // ************************************************
  gsinner=   gsl_integration_workspace_alloc (1000);

  // double gslinteg=0.0;
  // double integ_err;

  // gsl_function gslfun;
  // gslfun.function=&myfun;

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

  double pi=3.141592653589793;
  double xmin, xmax;
  double answer = 0.0;

  int i;

  xmin=2;
  xmax=100;

  // ***********************************************
  double T0=1;
  double ne0=1e5;

  struct my_f_params ptest;
  ptest.T=T0;
  ptest.ne=ne0;

  // *************************************************
  int method=GSL_INTEG_GAUSS61;
  gsl_function gslfun;
  gslfun.function=&myfun_gsl;
  gslfun.params=&ptest;
  gsl_integration_workspace * gswork=NULL;
  gswork=   gsl_integration_workspace_alloc (1000);
  double gslinteg=0.0;
  double integ_err;

  gettimeofday(&start, NULL);
  gsl_integration_qag(&gslfun, xmin, xmax, 0.0, 1e-8, 1000,1,
                       gswork, &gslinteg, &integ_err);
  gettimeofday(&end, NULL);
  gsltime = ((end.tv_sec  - start.tv_sec) * 1000000u +
              end.tv_usec - start.tv_usec) /1.e6;

  printf("gslres:%.12e, gsltime:%.4e\n", gslinteg, gsltime);
  // *************************************************

  stack_t stack;  //global stack
  create_stack(&stack, sizeof(work_gkq));

  gettimeofday(&start, NULL);
  double integ=gkq_double(myfun,inner_integrand2, xmin, xmax, 1e-2, &ptest,stack);

  gettimeofday(&end, NULL);
  free(stack->elements);
  free(stack);

  double partime = ((end.tv_sec  - start.tv_sec) * 1000000u +
                     end.tv_usec - start.tv_usec) /1.e6;

  printf("my:%.12e, partime:%.4e\n", integ,partime);

  return 0;
}

/******************************************
 * create new stack
 ******************************************/
void create_stack(
             stack_t* stack,     /* stack to create */
             int element_size)   /* size of a stack element */
{
  int initial_size = INITIAL_STACK_SIZE;

  /* allocate memory for new stack struct */
  (*stack) = (stack_t) malloc(sizeof(struct stack_s));
  if (!(*stack)){
    fprintf(stderr, "error: could not allocate memory for stack.. Abort.\n");
    exit(1);
  }

  /* allocate memory for stack elements */
  (*stack)->elements = (void*) malloc(element_size * initial_size);
  (*stack)->mem_reserve = initial_size;
  if (!(*stack)->elements){
    fprintf(stderr, "error: could not allocate memory for stack.. Abort.\n");
    exit(1);
  }

  (*stack)->el_size = element_size;
  (*stack)->el_count = 0;

}

/*****************************************
 * check if the stack is empty
 *****************************************/
int empty_stack(stack_t stack)
{
  //MYMOD
  if(stack==NULL)
    return 1;
  //ENDOF MYMOD
  return stack->el_count <= 0;
}

/*****************************************
 * push a element on stack
 *****************************************/
void push_stack(
           stack_t stack,    /* target stack */
           void* element)    /* element to push */
{
  int i, new_reserve;
  int log2_count;

  /* check if we need more memory for stack */
  if (stack->el_count >= stack->mem_reserve)
  {

      /* calculate new size for the stack
         it should be a power of two */
      // for (i = stack->el_count, log2_count = 0; i > 0;  i>>1, log2_count++);
      //   new_reserve = 1 << log2_count;

      log2_count=0;
      for(i=stack->el_count;i>0; i=i*0.5) //i>>1)
      {
        log2_count++;
      }
      new_reserve= 1 << log2_count;

      /* reallocate memory for phase thread tables
         and nullify new values */
      stack->elements = (void *) realloc(stack->elements,
                      stack->el_size * new_reserve);
      if (!stack->elements){
        fprintf(stderr, "error: can't reallocate stack.. Aborting\n");
        exit(1);
      }

      stack->mem_reserve = new_reserve;
  }

  /* now push the element on top of the stack */
  memcpy((char*)stack->elements + stack->el_count*stack->el_size,
            element, stack->el_size);
  stack->el_count++;

}

/*****************************************
 * pop an element from stack
 *****************************************/
void pop_stack(
          stack_t stack,    /* target stack */
          void* element)    /* where poped el. should be stored */
{
  if (stack->el_count <= 0){
    fprintf(stderr, "error: trying to pop from empty stack.\n");
    exit(2);
  }

  stack->el_count--;
  memcpy(element,
         (char*)stack->elements + stack->el_count*stack->el_size,
         stack->el_size);

}

// ***************************************************************************
// *      Gauss-kronard quadrature
// ***************************************************************************
work_gkq gkq_create_inner_task(double (*f)(double, struct my_f_params*), double a, double b, double TOL, struct my_f_params* p)
{
  work_gkq work;

  work.f=f;
  struct my_f_params* pinner;
  pinner=(struct my_f_params*) malloc(sizeof(struct my_f_params)); //ACHTUNG: Muss wieder gefreed werden
  pinner->mu= p->mu;
  pinner->T=p->T;
  pinner->ne=p->ne;
  work.p=pinner;
  work.a=a;
  work.b=b;
  work.subtasks_left=0;//beliebig!

  double m=0.5*(a+b);
  double h=0.5*(b-a);
  double y[13];
  double fa=y[0]=f(a,p);
  double fb=y[12]=f(b,p);
  int i;
  for(i=1;i<12;i++)
    y[i]=f(m+xgkq[i]*h,p);

  double I_4= (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));  // 4-point gauss-lobatto
  double I_7= (h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+  // 7-point kronrod
               625.0*(y[4]+y[8])+672.0*y[6]);

  double I_13= h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
                  0.188821573960182*(y[3]+y[9])+0.199773405226859*(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
                  0.242611071901408*y[6]);   //13-point Kronrod

  double Err1=fabs(I_7-I_13);
  double Err2=fabs(I_4-I_13);

  double r=(Err2 != 0.0) ? Err1/Err2 : 1.0;
  double toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;

  if(I_13 == 0)
    I_13=b-a;
  I_13=fabs(I_13);

// printf("create inner task: a:%f, b:%f, I4:%f,I7:%f,I13:%f\n",a,b,I_4,I_7,I_13);
  work.toler = toler;
  work.I_13=I_13;
  work.fa=fa;
  work.fb=fb;

  work.is_parent=0;
  for(i=0;i<5;i++)
    work.inner_integrals[i]=0.0;
  return work;

}

double gkq_double(double (*f)(double, struct my_f_params*), double (*finner)(double, struct my_f_params*),
                  double a, double b,
                  double TOL, struct my_f_params* p,stack_t stack)
{
  //1st integration
  //create parent task

  double result=0.0;
// *********************************************

  double m=0.5*(a+b);
  double h=0.5*(b-a);
  int i;

  double y[13];

  stack_t st;
  create_stack(&st, sizeof(work_gkq));

  double integ_bnd[2];

  get_integ_bounds_inner(integ_bnd, a, p);
  y[0]=f(a,p)* gkq_single(finner,integ_bnd[0],integ_bnd[1],1e-3, p, st);

  for(i=1;i<12;i++)
  {
    get_integ_bounds_inner(integ_bnd, m+xgkq[i]*h, p);
    y[i]=f(m+xgkq[i]*h,p)*gkq_single(finner,integ_bnd[0],integ_bnd[1],1e-3,p,st);
  }

  get_integ_bounds_inner(integ_bnd, b, p);
  y[12]=f(b,p)*gkq_single(finner,integ_bnd[0],integ_bnd[1],1e-3, p, st);

  free(st->elements);
  free(st);

  double fa=y[0];
  double fb=y[12];

  double I_4= (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));  // 4-point gauss-lobatto
  double I_7= (h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+  // 7-point kronrod
               625.0*(y[4]+y[8])+672.0*y[6]);

  double I_13= h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
                  0.188821573960182*(y[3]+y[9])+0.199773405226859*(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
                  0.242611071901408*y[6]);   //13-point Kronrod

  double Err1=fabs(I_7-I_13);
  double Err2=fabs(I_4-I_13);

  double r=(Err2 != 0.0) ? Err1/Err2 : 1.0;
  double toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;

  if(I_13 == 0)
    I_13=b-a;
  I_13=fabs(I_13);

  printf("First approx I_13:%.4e\n",I_13);

  //Prepare work and push onto stack
  work_gkq work;
  work.a = a;
  work.b = b;
  work.toler = toler;
  work.I_13=I_13;
  work.fa=fa;
  work.fb=fb;

  work.p=p;
  work.f=f;
  work.is_parent=1;
  work.task_nr=0; //fange bei 0 das zählen an

  for(i=0;i<5;i++)
    work.inner_integrals[i]=0.0;

  // ALLOC INNER WS FOR INTEGR.
  // gsl_integration_workspace *gsinner2=  gsl_integration_workspace_alloc (1000);

  stack_t stack_inner;
  create_stack(&stack_inner, sizeof(work_gkq));
  work.stack_inner=stack_inner;

  //push work on inner stack
  //for comp of 5 inner intgrals
  work_gkq winner;
  double m_inner, h_inner, mll_inner,ml_inner, mr_inner, mrr_inner;

  m_inner=  (work.a+work.b)/2;
  h_inner=  (work.a+work.b)/2;
  mll_inner=(m_inner-alpha*h_inner);
  ml_inner= (m_inner-beta*h_inner);
  mr_inner= (m_inner+beta*h_inner);
  mrr_inner=(m_inner+alpha*h_inner);

  //1st subtask
  get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
  winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0], integ_bnd[1],1e-3, work.p);
  winner.subtask_nr=0;
  winner.task_nr=work.task_nr;

  push_stack(work.stack_inner,&winner);

  //2nd subtask
  get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
  winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0], integ_bnd[1],1e-3, work.p);
  winner.subtask_nr=1;
  winner.task_nr=work.task_nr;
  push_stack(work.stack_inner,&winner);

  //3rd subtask
  get_integ_bounds_inner(integ_bnd, m_inner, work.p);
  winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0], integ_bnd[1],1e-3, work.p);
  winner.subtask_nr=2;
  winner.task_nr=work.task_nr;
  push_stack(work.stack_inner,&winner);

  //4th subtask
  get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
  winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0], integ_bnd[1],1e-3, work.p);
  winner.subtask_nr=3;
  winner.task_nr=work.task_nr;
  push_stack(work.stack_inner,&winner);

  //5th subtask
  get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
  winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0], integ_bnd[1],1e-3, work.p);
  winner.subtask_nr=4;
  winner.task_nr=work.task_nr;
  push_stack(work.stack_inner,&winner);

  //push outer integtand work on outer stack
  work.subtasks_left=5;
  push_stack(stack, &work);

/*
  int elemsinner=work.stack_inner->el_count;
  int elout=stack->el_count;

  work_gkq* top=(work_gkq*) stack->elements +(stack->el_count-1); //get top element

  int inner2=0;
  inner2=top->stack_inner->el_count;
*/
  // printf("elemsinner:%d,ou:%d,in2:%d\n",elemsinner,elout, inner2);

  // int in2=  top->stack_inner->el_count;

  result=gkq_adapt_double(stack);//,stack);

  // result=gkq_adapt_serial(f,a,b,fa,fb,toler,I_13, p);
  // gsl_integration_workspace_free(gsinner2);

  // free(stack->elements);

  return result;
}

double gkq_single(double (*f)(double, struct my_f_params*), double a, double b,
                       double TOL, struct my_f_params* p,stack_t stack)
{
  double result=0.0;
// *********************************************

  double m=0.5*(a+b);
  double h=0.5*(b-a);
  int i;

  double y[13];
  double fa=y[0]=f(a,p);
  double fb=y[12]=f(b,p);

  for(i=1;i<12;i++)           //hier müssen direkt 13 integrale berechnet werden!!!
    y[i]=f(m+xgkq[i]*h,p);    //auch das sollte parallel erfolgen (pragma parfor?)

  double I_4= (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));  // 4-point gauss-lobatto
  double I_7= (h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+  // 7-point kronrod
               625.0*(y[4]+y[8])+672.0*y[6]);

  double I_13= h*(0.0158271919734802*(y[0]+y[12])+0.0942738402188500*(y[1]+y[11])+0.155071987336585*(y[2]+y[10])+
                  0.188821573960182*(y[3]+y[9])+0.199773405226859*(y[4]+y[8])+0.224926465333340*(y[5]+y[7])+
                  0.242611071901408*y[6]);   //13-point Kronrod

  double Err1=fabs(I_7-I_13);
  double Err2=fabs(I_4-I_13);

  double r=(Err2 != 0.0) ? Err1/Err2 : 1.0;
  double toler=(r > 0.0 && r < 1.0) ? TOL/r : TOL;

  if(I_13 == 0)
    I_13=b-a;
  I_13=fabs(I_13);

  //Prepare work and push onto stack
  work_gkq work;
  work.a = a;
  work.b = b;
  work.toler = toler;
  work.I_13=I_13;
  work.fa=fa;
  work.fb=fb;
  work.p=p;
  work.a=a;
  work.f=f;

  push_stack(stack, &work);
  result=gkq_adapt_single(stack);//,stack);

  return result;
}

double gkq_adapt_serial(double (*f)(double, struct my_f_params*), double a, double b,
                        double fa,double fb, double toler,double I_13, struct my_f_params* p)
{
  double m = (a+b)/2;
  double h = (b-a)/2;
  double mll=m-alpha*h;
  double ml=m-beta*h;
  double mr=m+beta*h;
  double mrr=m+alpha*h;

  double fmll=f(mll,p);
  double fml=f(ml,p);
  double fm=f(m,p);
  double fmr=f(mr,p);
  double fmrr=f(mrr,p);
  double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
  double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

  if (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr)
  {
    if ((mll <= a || b <= mrr) && !terminate_serial) //Error
     {
      // out_of_tolerance=true; // Interval contains no more machine numbers
      printf("OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e\n", mll,b,b,mrr);
      terminate_serial=1;
     }

// printf("me ok:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f\n", omp_get_thread_num(), a,b,toler,I_4,I_7);

    return I_7;
  }
  else
  {

    // printf("me NOOOO:%d, a:%f,b:%f, tler:%.5e,I_4:%f,I_7:%f\n", omp_get_thread_num(), a,b,toler,I_4,I_7);

    return  gkq_adapt_serial(f, a,mll,fa,fmll,toler,I_13,p) +
            gkq_adapt_serial(f, mll,ml,fmll,fml,toler,I_13,p) +
            gkq_adapt_serial(f, ml,m,fml,fm,toler,I_13,p) +
            gkq_adapt_serial(f, m,mr,fm,fmr,toler,I_13,p) +
            gkq_adapt_serial(f, mr,mrr,fmr,fmrr,toler,I_13,p) +
            gkq_adapt_serial(f, mrr,b,fmrr,fb,toler,I_13,p);
  }

}

double gkq_adapt_single(stack_t stack)
{
  work_gkq work;
  // work.iter=0;
  int ready, idle, busy;
  double integral_result = 0.0;

  busy = 0;
  terminate_gkq=0;

#pragma omp parallel default(none) \
    shared(stack, integral_result,busy) \
    private(work, idle, ready)
  {
// printf("me:%d, err:%d\n",omp_get_thread_num(),simpson_error);

    ready = 0;
    idle = 1;

    while(!ready) // && !terminate_gkq)//  && !simpson_error) //<-- so NICHT!
    {
      #pragma omp critical (stack)
      {
        if (!empty_stack(stack))
        {
          /* we have new work */
          pop_stack(stack, &work);
          if (idle)
          {
            /* say others i'm busy */
            busy += 1;
            idle = 0;
          }
        }
        else
        {
          /* no new work on stack */
          if (!idle){
            busy -= 1;
            idle = 1;
          }

          /* nobody has anything to do; let us leave the loop */
          if (busy == 0)
          {
            ready = 1;
          }
        }
      } /* end critical(stack) */

      if (idle)
        continue; //if ready==1 --> leave loop
// double I_prev=work.I_prev;

double (*f)(double, struct my_f_params*)=work.f;

      double a = work.a;
      double b = work.b;
      double toler = work.toler;
      double I_13=work.I_13;
      double fa=work.fa;
      double fb=work.fb;
      // int iter=work.iter;
      // double *y= work.y; // brauch ich nicht!
      struct my_f_params * p = work.p;

      double m = (a+b)/2;
      double h = (b -a)/2;
      double mll=m-alpha*h;
      double ml=m-beta*h;
      double mr=m+beta*h;
      double mrr=m+alpha*h;

      double fmll=f(mll,p);
      double fml=f(ml,p);
      double fm=f(m,p);
      double fmr=f(mr,p);
      double fmrr=f(mrr,p);
      double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
      double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

// if(myid==1)
//   printf("I_7:%.4e, I_13:%.4e,I_4:%.4e, minus:%.4e, to:%.4e\n",I_7,I_13,I_4,I_7-I_4, toler*I_13);
// int maxiter=50; //max. subdivisions
// double abstol=1e-30;
// work.I_prev=I_7; // für abstolcheck in nächster recursion

      if (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr) // || iter > maxiter || fabs(I_7-I_prev) < abstol )
      {
        if ((mll <= a || b <= mrr)) //Error
         {
           // out_of_tolerance=true; // Interval contains no more machine numbers
           // printf("OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e,I_7-I_4:%.4e, tol:%.4e,I_13:%.4e\n",
           //        mll,b,b,mrr,I_7-I_4, toler*I_13,I_13);

         }
        #pragma omp critical (integral_result)
        {
          integral_result += I_7;      //Terminate recursion.
        }
      }
      else  //subdivide interval and push new work on stack
      {
        #pragma omp critical (stack)
        {
          // work.iter=iter+1;

          work.a=a;
          work.b=mll;
          work.fa=fa;
          work.fb=fmll;
          push_stack(stack, &work);

          work.a=mll;
          work.b=ml;
          work.fa=fmll;
          work.fb=fml;
          push_stack(stack, &work);

          work.a=ml;
          work.b=m;
          work.fa=fml;
          work.fb=fm;
          push_stack(stack, &work);

          work.a=m;
          work.b=mr;
          work.fa=fm;
          work.fb=fmr;
          push_stack(stack, &work);

          work.a=mr;
          work.b=mrr;
          work.fa=fmr;
          work.fb=fmrr;
          push_stack(stack, &work);

          work.a=mrr;
          work.b=b;
          work.fa=fmrr;
          work.fb=fb;

          push_stack(stack, &work);

        } // pragma critical stack
      }   // else ..non-acceptable error
    } // while
  } /* end omp parallel */
  return integral_result;
}

double gkq_adapt_double(stack_t stack)
{
  work_gkq work;
  work_gkq* pwork_outer;

  int ready, idle, busy;
  double integral_result = 0.0;
  busy = 0;
  int myid;
  int elcnt; //nr of outer task elements on stack

#pragma omp parallel default(none)  \
    shared(stack, integral_result,busy,iter,elcnt) \
    private(work, idle, ready,myid,pwork_outer)
  {

    myid=omp_get_thread_num();

    ready = 0;
    idle = 1;

    while(!ready) // && !terminate_gkq)//  && !simpson_error) //<-- so NICHT!
    {

// printf("\n NEXT, curapprox:%.4e\n",integral_result);
      #pragma omp critical (stack)
      {
        //pointer to outer work element
        //work_gkq* pwork=(work_gkq*) stack->elements + stack->el_count*stack->el_size;

        elcnt=stack->el_count;

// if(elcnt>1)
//  getchar();

        stack_t stack_inner;
        if(elcnt>0)
        {
          pwork_outer=(work_gkq*) stack->elements +(elcnt-1); //get top element

          //IST WOHL SUBOPTIMAL
            while(pwork_outer->task_nr > 0 && pwork_outer->subtasks_left==0) //this task is complete, goto next outer task
            {
              printf("myid:%d, outer task nr:%d, complete:subs:%d, decrement\n",myid,pwork_outer->task_nr,pwork_outer->subtasks_left);
              pwork_outer--;
            }
            stack_inner=pwork_outer->stack_inner;
        }
        else
        {
          stack_inner=NULL;
          pwork_outer=NULL;
        }

        printf("myid:%d,elcnt:%d,sinner:%p,pwout:%p\n",myid,elcnt,stack_inner,pwork_outer);

        { //stackinner gehört zu work, und work ist private!
          if(!empty_stack(stack_inner))
          {
            printf("myid:%d,iter:%d, inner stack not empty,pop..\n",myid,iter);
            pop_stack(stack_inner, &work);
//if letztes inner-work elem, blockiere outer work elem bis ich das integral hab

            iter++;

            if (idle)
            {
              busy += 1;
              idle = 0;
            }
          }
          else //inner stack is empty, pop from outer
          {
            printf("myid:%d,iter:%d, inner stack is empty?\n",myid,iter);
            if (!empty_stack(stack)) //work elem ist blockiert solange inner integs nicht vollständig
            {
              // printf("myid:%d,iter:%d, inner stack empty. work on outer with cnt=%d\n",myid,iter,pwork_outer->subtasks_left);
              if(pwork_outer->subtasks_left==0)
              {
                printf("myid:%d,iter:%d, outer stack not empty\n",myid,iter);
                pop_stack(stack, &work);
                iter++;
                if (idle)
                {
                  busy += 1;
                  idle = 0;
                }
              }
              else
              {
                printf("myid:%d,iter:%d, inner integs not complete:%d\n",myid,iter,pwork_outer->subtasks_left);
                if (!idle)
                {
                  busy -= 1;
                  idle = 1;
                }

              }
            }
            else //auch outer stack ist leer
            {
              printf("myid:%d,iter:%d, outer stack empty too\n",myid,iter);
              if (!idle)
              {
                busy -= 1;
                idle = 1;
              }
              if (busy == 0)
              {
                ready = 1;
              }
            }
          }
        } // critical inner stack
      } //critical outer stack

      if (idle)
      {
        printf("myid:%d,iter:%d, noth8ing to do\n",myid,iter);
        continue; //if ready==1 --> leave loop
      }

      //work on inner tasks first,  if available
      if(work.is_parent==0)
      {
        printf("myid:%d,iter:%d,tasknr:%d,innertask:%d, left subs:%d\n",
                myid,iter,work.task_nr,work.subtask_nr,pwork_outer->subtasks_left); //subtasks nr nicht immer aktuell!

// getchar();

        double (*f)(double, struct my_f_params*) = work.f;
        double a = work.a;
        double b = work.b;
        double toler = work.toler;
        double I_13=work.I_13;
        double fa=work.fa;
        double fb=work.fb;
        // double *y= work.y; // brauch ich nicht!
        struct my_f_params * p = work.p;

        double m = (a+b)/2;
        double h = (b -a)/2;
        double mll=m-alpha*h;
        double ml=m-beta*h;
        double mr=m+beta*h;
        double mrr=m+alpha*h;

        double fmll=f(mll,p);
        double fml=f(ml,p);
        double fm=f(m,p);
        double fmr=f(mr,p);
        double fmrr=f(mrr,p);
        double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
        double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

// printf("myid:%d,iter:%d non parent checkpoint 1 passed\n",myid,iter);

        if (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr)
        {
          if ((mll <= a || b <= mrr))// && !terminate_gkq) //Error
           {
             printf("OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e\n", mll,a,b,mrr);
           }

// printf("myid:%d,iter:%d non parent checkpoint 2.1 passed\n",myid,iter);

          int tasknr=work.subtask_nr;

          printf("myid:%d,iter:%d, winner error acceptable: I_7:%.4e, I_4:%.4e, I_13:%.4e,toler*I_13:%.4e,toler:%.4e\n",
                  myid,iter,I_7,I_4, I_13,toler*I_13,toler);

          #pragma omp critical (innerinteg) // workers greifen auf selbes array zu.nur an anderer Stelle..evtl. besser 6 versch.doubles statt array?
          {
            double *inner_integrals=pwork_outer->inner_integrals;
            inner_integrals[tasknr]+=I_7;
            pwork_outer->subtasks_left=pwork_outer->subtasks_left-1; //hier sollte subtasks_left aktuell sein

            printf("myid:%d,iter:%d, inner_integral[%d]:%.4e added:%.4e,subtask left:%d\n",
                    myid,iter,tasknr,inner_integrals[tasknr],I_7,pwork_outer->subtasks_left);
          }
        }
        else  //subdivide interval and push new work on stack
        {

          stack_t stack_inner;
          #pragma omp critical //(stack) //(stack_inner)
          {
            stack_inner =  pwork_outer->stack_inner;
            pwork_outer->subtasks_left=pwork_outer->subtasks_left-1;//weil dieser stack durch 6 andere ersetzt wird

            work.task_nr=pwork_outer->task_nr; //inner task soll auch wissen wer sein outer task is (bisher nur für debug)
            //subtasknr bleibt dieselbe (bezogen auf parent task)
            work.a=a;
            work.b=mll;
            work.fa=fa;
            work.fb=fmll;

            push_stack(stack_inner, &work);

            work.a=mll;
            work.b=ml;
            work.fa=fmll;
            work.fb=fml;

            push_stack(stack_inner, &work);

            work.a=ml;
            work.b=m;
            work.fa=fml;
            work.fb=fm;

            push_stack(stack_inner, &work);

            work.a=m;
            work.b=mr;
            work.fa=fm;
            work.fb=fmr;

            push_stack(stack_inner, &work);

            work.a=mr;
            work.b=mrr;
            work.fa=fmr;
            work.fb=fmrr;

            push_stack(stack_inner, &work);

            work.a=mrr;
            work.b=b;
            work.fa=fmrr;
            work.fb=fb;

            push_stack(stack_inner, &work);

            pwork_outer->subtasks_left=pwork_outer->subtasks_left+6;

            printf("myid:%d,iter:%d, inner_integral approx not accurate..subdivided: I_4:%f,I_7:%f,I_13:%f,left subs:%d\n",
                    myid,iter,I_4,I_7,I_13,pwork_outer->subtasks_left);

          } // pragma critical stack_inner
        }   // else ..non-acceptable error
      }// !isparent
      else  //parent task without any inner tasks
      {

printf("myid:%d,iter:%d, work is parent,subs left:%d\n",myid,iter,work.subtasks_left);

        // stack_t stack_inner = work.stack_inner; //brauch ich hier nicht
        // free(stack_inner->elements);
        // free(stack_inner);

        double (*f)(double, struct my_f_params*) = work.f;
        double a = work.a;
        double b = work.b;
        double toler = work.toler;
        double I_13=work.I_13;
        double fa=work.fa;
        double fb=work.fb;
        // double *y= work.y; // brauch ich nicht!
        struct my_f_params * p = work.p;

        double m = (a+b)/2;
        double h = (b -a)/2;
        double mll=m-alpha*h;
        double ml=m-beta*h;
        double mr=m+beta*h;
        double mrr=m+alpha*h;

        //the inner intergrals are pre-calculated and saved in an array
        //the function only calculates a prefactor...
        double fmll=f(mll,p)*work.inner_integrals[0];
        double fml=  f(ml,p)*work.inner_integrals[1];
        double fm=   f(m,p)*work.inner_integrals[2];
        double fmr=  f(mr,p)*work.inner_integrals[3];
        double fmrr= f(mrr,p)*work.inner_integrals[4];

        double I_4=h/6.0*(fa+fb+5.0*(fml+fmr));   // 4-point Gauss-Lobatto formula.
        double I_7=h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);

        if (fabs(I_7-I_4) <= toler*I_13 || mll <= a || b <= mrr)
        {
          if ((mll <= a || b <= mrr))// && !terminate_gkq) //Error
           {
             // out_of_tolerance=true; // Interval contains no more machine numbers
             printf("outer task OUT OF TOLERANCE !!!, mll:%.4e, a:%.4e, b:%.4e, mrr:%.4e\n", mll,b,b,mrr);
           }
           //#pragma omp critical(integal_result)
           {
              #pragma omp atomic
              integral_result += I_7;
              printf("myid:%d,iter:%d, I7 added, integs: 0=%f,1=%f,2=%f,3=%f,4=%f\n", myid,iter,
                    work.inner_integrals[0],work.inner_integrals[1],work.inner_integrals[2],
                    work.inner_integrals[3],work.inner_integrals[4]);
           }
        }
        else  //subdivide interval and push new work on stack
        {

printf("myid:%d,iter:%d, iouter_integral approx not accurate: I_7=%.4e,I_4=%.4e,I_13=%.4e, toler*I_13:%.4e\n"
        "integs: 0=%f,1=%f,2=%f,3=%f,4=%f..subdivide\n",
        myid,iter,I_7,I_4,I_13,I_13*toler,work.inner_integrals[0],work.inner_integrals[1],work.inner_integrals[2],
                    work.inner_integrals[3],work.inner_integrals[4] );

          stack_t stack_inner = work.stack_inner;
          #pragma omp critical (stack)
          {
            //create sub-stack for inner-integrals

            work.subtasks_left=5; //for all new outer work elements

            create_stack(&stack_inner, sizeof(work_gkq));
            work.is_parent=1;

            work_gkq winner;

            //outer task elem
            work.a=a;
            work.b=mll;
            work.fa=fa;
            work.fb=fmll;

            double m_inner=  (work.a+work.b)/2;
            double h_inner=  (work.a+work.b)/2;
            double mll_inner=(m_inner-alpha*h_inner);
            double ml_inner= (m_inner-beta*h_inner);
            double mr_inner= (m_inner+beta*h_inner);
            double mrr_inner=(m_inner+alpha*h_inner);

            //Jedes subinterval bekommt 5 inner tasks für die inneren integrale!

            int curcnt=stack->el_count;
            work.task_nr=curcnt; //weil ich bei 0 anfange zu zählen und push_stack gleich el_count incrementiert

            double integ_bnd[2];

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.task_nr=work.task_nr;
            winner.subtask_nr=0;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            push_stack(stack, &work);

            // *************************
            // *   2nd subinterval
            // **************************
            work.a=mll;
            work.b=ml;
            work.fa=fmll;
            work.fb=fml;
            work.task_nr++;

            m_inner=  (work.a+work.b)/2;
            h_inner=  (work.a+work.b)/2;
            mll_inner=(m_inner-alpha*h_inner);
            ml_inner= (m_inner-beta*h_inner);
            mr_inner= (m_inner+beta*h_inner);
            mrr_inner=(m_inner+alpha*h_inner);

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=0;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            // #pragma omp critical (stack)
            // {
              push_stack(stack, &work);
            // }

            // *************************
            // *   3rd subinterval
            // **************************
            work.a=ml;
            work.b=m;
            work.fa=fml;
            work.fb=fm;
            work.task_nr++;

            m_inner=  (work.a+work.b)/2;
            h_inner=  (work.a+work.b)/2;
            mll_inner=(m_inner-alpha*h_inner);
            ml_inner= (m_inner-beta*h_inner);
            mr_inner= (m_inner+beta*h_inner);
            mrr_inner=(m_inner+alpha*h_inner);

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=0;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2, integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            // #pragma omp critical (stack)
            // {
              push_stack(stack, &work);
            // }

            // *************************
            // *   4th subinterval
            // **************************
            work.a=m;
            work.b=mr;
            work.fa=fm;
            work.fb=fmr;
            work.task_nr++;

            m_inner=  (work.a+work.b)/2;
            h_inner=  (work.a+work.b)/2;
            mll_inner=(m_inner-alpha*h_inner);
            ml_inner= (m_inner-beta*h_inner);
            mr_inner= (m_inner+beta*h_inner);
            mrr_inner=(m_inner+alpha*h_inner);

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=0;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            // #pragma omp critical (stack)
            {
              push_stack(stack, &work);
            }

            // *************************
            // *   5th subinterval
            // **************************
            work.a=mr;
            work.b=mrr;
            work.fa=fmr;
            work.fb=fmrr;
            work.task_nr++;

            m_inner=  (work.a+work.b)/2;
            h_inner=  (work.a+work.b)/2;
            mll_inner=(m_inner-alpha*h_inner);
            ml_inner= (m_inner-beta*h_inner);
            mr_inner= (m_inner+beta*h_inner);
            mrr_inner=(m_inner+alpha*h_inner);

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=0;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            // #pragma omp critical (stack)
            {
              push_stack(stack, &work);
            }

            // *************************
            // *   6th subinterval
            // **************************
            work.a=mrr;
            work.b=b;
            work.fa=fmrr;
            work.fb=fb;
            work.task_nr++;

            m_inner=  (work.a+work.b)/2;
            h_inner=  (work.a+work.b)/2;
            mll_inner=(m_inner-alpha*h_inner);
            ml_inner= (m_inner-beta*h_inner);
            mr_inner= (m_inner+beta*h_inner);
            mrr_inner=(m_inner+alpha*h_inner);

            //1st subtask
            get_integ_bounds_inner(integ_bnd, mll_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=0;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //2nd subtask
            get_integ_bounds_inner(integ_bnd, ml_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=1;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //3rd subtask
            get_integ_bounds_inner(integ_bnd, m_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=2;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //4th subtask
            get_integ_bounds_inner(integ_bnd, mr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=3;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            //5th subtask
            get_integ_bounds_inner(integ_bnd, mrr_inner, work.p);
            winner=gkq_create_inner_task(inner_integrand2,integ_bnd[0],integ_bnd[1],1e-3, work.p);
            winner.subtask_nr=4;
            winner.task_nr=work.task_nr;
            push_stack(work.stack_inner,&winner);

            // #pragma omp critical (stack)
            {
              push_stack(stack, &work);
            }

            printf("myid:%d,iter:%d, 6 subinvtervals with 5 subtasks each pushed to stack outer\n",myid,iter);

          } // pragma critical
        }   // else ..non-acceptable error

      }

    } // while
  } /* end omp parallel */

  return integral_result;
}

//SIMPSON

  // *****************************************************************************************************************************************
// double integral2(double (*f)(double, struct my_f_params*), stack_t stack)
// {
//   work_t work;
//   int ready, idle, busy;
//   double integral_result = 0.0;

//   busy = 0;

// #pragma omp parallel default(none) \
//     shared(stack, integral_result,f,busy,simpson_error) \
//     private(work, idle, ready)
//   {

// // printf("me:%d, err:%d\n",omp_get_thread_num(),simpson_error);

//     ready = 0;
//     idle = 1;

//     while(!ready  && !simpson_error) //<-- so NICHT!
//     {
// #pragma omp critical (stack)
//       {
//         if (!empty_stack(stack))
//         {
//           /* we have new work */
//           pop_stack(stack, &work);
//           if (idle)
//           {
//             /* say others i'm busy */
//             busy += 1;
//             idle = 0;
//           }
//         }
//         else{
//           /* no new work on stack */
//           if (!idle){
//             busy -= 1;
//             idle = 1;
//           }

//           /* nobody has anything to do; let us leave the loop */
//           if (busy == 0)
//           {
//             ready = 1;
//           }
//         }
//       } /* end critical(stack) */

//       if (idle)
//         continue; //if ready==1 --> leave loop

//       double b = work.b;
//       double a = work.a;
//       double tol = work.tol;
//       double S=work.S; // previous TOTAL integral
//       double fa=work.fa;
//       double fb=work.fb;
//       double fm=work.fm;
//       int    rec=work.rec;

//       double h = (b - a)/2;
//       double mid = (a+b)/2;
//       double lm=(a+mid)/2;
//       double rm=(mid+b)/2;

//       double flm=f(lm,work.p);
//       double frm=f(rm,work.p);

//       double Sl=h/6*(fa+4*flm+fm);
//       double Sr=h/6*(fm+4*frm+fb);

//       double delta=Sl+Sr-S;

//       // serious numerical trouble: it won't converge
//       if ((tol/2 == tol) || fabs(a-lm) <= tol) // || tol < tolmax)
//       {
//              simpson_error = 1;
// #pragma omp critical (integral_result)
//              integral_result = S;
//       }
//       //need something against spurious convergence
//       //for now: iter > 5 (c.f. numerical recipes)
//       if( work.iter > 5 &&  (rec <= 0 || fabs(delta) <= 15*tol))  //error acceptable
//       {
// #pragma omp critical (integral_result)
//         integral_result += Sl+Sr+delta/15;
//       }
//       else // error not acceptable
//       {
//         //push new subintervals to stack
//         work.a = a;
//         work.b = mid;
//         work.tol = tol/2;
//         work.S = Sl;
//         work.fa=fa;
//         work.fb=fm;
//         work.fm=flm;
//         work.rec=rec-1;
//         work.iter=work.iter+1;

//   #pragma omp critical (stack)
//         {
//           //LEFT
//           push_stack(stack, &work);
//           //prepare RIGHT side and push to stack
//           work.a = mid;
//           work.b = b;
//           work.tol = tol/2;
//           work.S=Sr;
//           work.fa=fm;
//           work.fb=fb;
//           work.fm=frm;
//           work.rec=rec-1;
//           push_stack(stack, &work);

//         }
//       }
//     } /* while */
//   } /* end omp parallel */

//   return integral_result;
// }

max_stars_repo_name: fmqeisfeld/IMD

max_issues_repo_name: fmqeisfeld/IMD

max_forks_repo_name: fmqeisfeld/IMD


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>

#include <omp.h>

#include "util.h"

int main(int argc, char *argv[]) {
    char* input_fileName1 = argv[1];
    char* input_fileName2 = argv[2];
    int N_doc_bg = atoi(argv[3]);
    int N_kw = atoi(argv[4]);
    int N_obs = atoi(argv[5]);
    int N_iter = atoi(argv[6]);
    int m = atoi(argv[7]);
    double p = 0.89;
    double q = atof(argv[8]);
    char* output_fileName = argv[9];
    int N_doc = 480000/2;

    int* matrix = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* matrix_bg = (int*) malloc(sizeof(int) * N_kw * N_kw);
    int* matrix_padded = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* true_index = (int*) malloc(sizeof(int) * N_kw);
    int* permutation = (int*) malloc(sizeof(int) * N_obs);
    gsl_matrix* matrix_obs;

    // Setup
    for (int round = 0; round < 10; round++)
    {
        char input_fileName1_extend[40];
        char input_fileName2_extend[40];
        sprintf(input_fileName1_extend, "%s%d", input_fileName1, round);
        sprintf(input_fileName2_extend, "%s%d", input_fileName2, round);

        struct timeval tv1,tv2;
        gettimeofday(&tv1, NULL);
        read_matrix(&true_index, &matrix_bg, 1.0*N_doc/N_doc_bg, N_kw, input_fileName2_extend);
        read_matrix(&true_index, &matrix, 1.0, N_obs, input_fileName1_extend);
        gettimeofday(&tv2, NULL);
        printf("Reading done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
        fflush(stdout);

        for (int iter = 0; iter < 10; iter++)
        {
            printf("Run %d\n", iter);
            matrix_obs = gsl_matrix_alloc(N_obs, N_obs);

            gettimeofday(&tv1, NULL);
            pad_matrix(&matrix_padded, &matrix, m, p, q, N_obs, N_doc);
            gettimeofday(&tv2, NULL);
            printf("Padding done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            gettimeofday(&tv1, NULL);
            observe_matrix(matrix_obs, &matrix_padded, N_obs);
            gettimeofday(&tv2, NULL);
            printf("Observed matrix generated: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            // Permute observed matrix randomly and attack

            gettimeofday(&tv1, NULL);
            attack(matrix_obs, &matrix_bg, &permutation, m, p, q, N_kw, N_obs, N_doc, N_iter);
            gettimeofday(&tv2, NULL);
            printf("Main attack done: %f.\n", (double) ((tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec)) / N_iter);
            fflush(stdout);

            //for (int ii = 0; ii < 20; ii++)
            //   printf("%d, %d, %d, %d, %d\n", permutation[ii], true_index[ii], matrix[ii*N_obs+ii], matrix_bg[true_index[ii]*N_kw+true_index[ii]], matrix_bg[permutation[ii]*N_kw+permutation[ii]]);

            char output_fileName_full[40];
            sprintf(output_fileName_full, "%s%d-%d", output_fileName, round, iter);
            print_result(output_fileName_full, &permutation, &true_index, N_obs);
            //sprintf(output_fileName_full, "%s-%d-full", output_fileName, iter);
            //print_full_result(output_fileName_full, permutation_matrix, N_kw);
        }
    }

    free(matrix);
    free(matrix_bg);
    free(matrix_padded);
    gsl_matrix_free(matrix_obs);
    return(0);
}

double log_score(int idx1, int idx2, gsl_matrix* matrix_obs, int** matrix, int** permutation, int m, double p, double q, int N_kw, int N_doc)
{
    int idx1_m = (*permutation)[idx1];
    int idx2_m = (*permutation)[idx2];

    double score = 0.0;
    if (idx1 == idx2)
    {
        int N1 = m * (*matrix)[idx1_m*N_kw + idx2_m];
        double N1_mean = p * N1;
        double N1_var = p * (1-p) * N1 + (m * N_doc - N1) * N1 * 1.0 / N_doc / m;

        int N2 = m * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);
        double N2_mean = q * N2;
        double N2_var = q * (1-q) * N2 + (m * N_doc - N2) * N2 * 1.0 / N_doc / m;

        double N3_var = 1.0 * m / N_doc * (*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);

        int k = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

        score = gsl_cdf_gaussian_P(floor(k - N1_mean - N2_mean) + 0.5, sqrt(N1_var + N2_var + N3_var));
        score -= gsl_cdf_gaussian_P(floor(k - N1_mean - N2_mean) - 0.5, sqrt(N1_var + N2_var + N3_var));
    }

    else
    {
        int N1 = m * (*matrix)[idx1_m*N_kw + idx2_m];
        double N1_mean = p * p * N1;
        double N1_var  = p * p * (1-p) * (1-p)* N1;

        int N2 = m * ((*matrix)[idx1_m*N_kw + idx1_m] + (*matrix)[idx2_m*N_kw + idx2_m] - 2*(*matrix)[idx1_m*N_kw + idx2_m]);
        double N2_mean = p * q * N2;
        double N2_var  = p * q * (1-p) * (1-q) * N2;

        int N3 = m * N_doc - N1 - N2;
        double N3_mean = q * q * N3;
        double N3_var  = q * q * (1-q) * (1-q) * N3;

        double N4_var = 1.0 * m / N_doc * (*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);

        int M = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

        score = gsl_cdf_gaussian_P(floor(M - N1_mean - N2_mean - N3_mean) + 0.5, sqrt(N1_var + N2_var + N3_var + N4_var));
        score -= gsl_cdf_gaussian_P(floor(M - N1_mean - N2_mean - N3_mean) - 0.5, sqrt(N1_var + N2_var + N3_var + N4_var));
    }

    if (score == 0)
        return(-500.0);
    return(log(score));
}

void attack(gsl_matrix* matrix_obs, int** matrix, int** permutation, int m, double p, double q, int N_kw, int N_obs, int N_doc, int N_iter)
{
    // Initialise data structures
    double* score_matrix = (double*) malloc(sizeof(double) * N_kw * N_kw);
    double* score_row1 = (double*) malloc(sizeof(double) * N_kw);
    double* score_row2 = (double*) malloc(sizeof(double) * N_kw);

    int* permutation_tmp = (int*) malloc(sizeof(int) * N_obs);
    int* permutation_inv = (int*) malloc(sizeof(int) * N_kw);

    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    struct timeval tv1,tv2;
    gettimeofday(&tv1, NULL);

    solution_initial(permutation, matrix_obs, matrix, m, p, q, N_kw, N_obs, N_doc);

    for (int ii = 0; ii < N_obs; ii++)
        permutation_tmp[ii] = (*permutation)[ii];
    for (int ii = 0; ii < N_kw; ii++)
        permutation_inv[ii] = -1;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_inv[permutation_tmp[ii]] = ii;

    // Compute initial score
    #pragma omp parallel for shared(score_matrix, matrix_obs, matrix)
    for (int ii = 0; ii < N_obs * N_obs; ii++)
        score_matrix[ii] = log_score((int) (ii / N_obs), ii % N_obs, matrix_obs, matrix, permutation, m, p, q, N_kw, N_doc);
    gettimeofday(&tv2, NULL);
    printf("Initial score computed: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));

    // Iterations of simulated annealing
    double temp = (double) N_kw;
    int N_stuck = 0;
    for (int iter = 0; iter < N_iter; iter++)
    {
        /* Status code */
	    if (iter % (N_iter / 10) == 0)
        {
            gettimeofday(&tv1, NULL);
            printf("Iteration: %d, %d, %d.\n", iter, N_stuck, (int) (tv1.tv_sec - tv2.tv_sec));
            fflush(stdout);
            gettimeofday(&tv2, NULL);
        }

	    if (N_stuck >= (N_iter / 20))
            iter = N_iter;

        int idx1, idx2;
        permutation_generation(&idx1, &idx2, &permutation_tmp, permutation, &permutation_inv, matrix_obs, matrix, m, p, q, N_kw, N_obs, N_doc);
        if (idx1 == idx2)
        {
            N_stuck++;
            continue;
        }

        int ii = 0;
        #pragma omp parallel for shared(score_row1)
        for (ii = 0; ii < N_obs; ii++)
            score_row1[ii] = log_score(idx1, ii, matrix_obs, matrix, &permutation_tmp, m, p, q, N_kw, N_doc);

        if (idx2 >= 0)
            #pragma omp parallel for shared(score_row2)
            for (ii = 0; ii < N_obs; ii++)
                score_row2[ii] = log_score(idx2, ii, matrix_obs, matrix, &permutation_tmp, m, p, q, N_kw, N_doc);

        double score_diff = 0;
        for (int ii = 0; ii < N_obs; ii++)
            score_diff += score_row1[ii];
        for (int ii = 0; ii < N_obs; ii++)
            score_diff -= score_matrix[idx1*N_obs + ii];
        if (idx2 >= 0)
        {
            for (int ii = 0; ii < N_obs; ii++)
                score_diff += score_row2[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_diff -= score_matrix[idx2*N_obs + ii];
        }

        // compute difference in score, with exponentiation
        score_diff = score_diff / temp;

        if (score_diff < -40)
            score_diff = 0;
        else if (score_diff > 0)
            score_diff = 1.01;
        else
            score_diff = exp(score_diff);

        if (gsl_ran_flat(r, 0, 1) < score_diff)
        {
            // Update the scores
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[idx1*N_obs + ii] = score_row1[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[ii*N_obs + idx1] = score_row1[ii];
            if (idx2 >= 0)
            {
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[idx2*N_obs + ii] = score_row2[ii];
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[ii*N_obs + idx2] = score_row2[ii];
            }

            // Update the permutation
            permutation_inv[(*permutation)[idx1]] = -1;
            (*permutation)[idx1] = permutation_tmp[idx1];
            permutation_inv[permutation_tmp[idx1]] = idx1;
            if (idx2 >= 0)
            {
                (*permutation)[idx2] = permutation_tmp[idx2];
                permutation_inv[permutation_tmp[idx2]] = idx2;
            }
            N_stuck = 0;
        }
        else
        {
            // Update the permutation
            permutation_tmp[idx1] = (*permutation)[idx1];
            if (idx2 >= 0)
                permutation_tmp[idx2] = (*permutation)[idx2];
            N_stuck += 1;
        }

        temp *= 0.995;
    }

    free(score_matrix);
    free(score_row1);
    free(score_row2);
    gsl_rng_free(r);
}

void print_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    int count2 = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    fprintf(fp, "%d\n", count);
    fclose(fp);
    printf("Success: %d/%d.\n", count, N_obs);
}

void print_full_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    fprintf(fp, "%d\n", count);
    fclose(fp);
}

ext: c

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
/*solver/fixed_point.c
* Solves a fixed point equation using
* sequence acceleration.
*
* Author: Benjamin Vatter j.
* email : benjaminvatterj@gmail.com
* date: 15 August 2015
*/

#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include "fixed_point.h"

/*
* Finds the fixed point of a function f
* Translated from SciPy's fixed_point function
*/
void fixed_point(int n, ff_function *f, double *x0,
                 double xtol, int maxiter, double *out)
{
	//out = malloc(sizeof(double)*n);
	double *p = malloc(sizeof(double)*n);
	double *p1, *p2;
	double d;
	int pass;
	int i, j;

	for (i=0; i<n; i++)
		out[i] = x0[i];

	for (i=0; i<maxiter; i++)
	{
		p1 = malloc(sizeof(double)*n);
		p2 = malloc(sizeof(double)*n);
		f->function(out, f->params, p1);
		f->function(p1, f->params, p2);

		pass = 1;
		for (j=0; j<n; j++)
		{
			d = p2[j] - 2.0*p1[j] + out[j];
			if (d==0) {
				p[j] = p2[j];
			} else {
				p[j] = out[j] - pow(p1[j] - out[j], 2.0) / d;
			}

			if (out[j] == 0 && fabs(p[j]) > xtol){
				pass = 0;
			}
			if (out[j] != 0 && (p[j] - out[j])/out[j] > xtol){
				pass = 0;
			}
			out[j] = p[j];
		}
		if (pass == 1) {
			free(p1);
			free(p2);
			break;
		}
		if(i+1 >= maxiter){
			free(p1);
			free(p2);
			pass = -1;
			break;
		}
		free(p1);
		free(p2);
	}

	free(p);
	// Free garbage result in case of no convergence
	if (pass == -1) {
		free(out);
	}
}


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_4_3.c
 * \brief Source file to optimize Runge-Kutta 4 steps 3rd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_4_3.h"

#define DEBUG_RK_4_3 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 4 steps 3rd order Runge-Kutta
 * method.
 */
int
rk_tb_4_3 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b32 (tb) = r[4];
  b43 (tb) = r[5];
  b42 (tb) = ((1.L / 3.L - b43 (tb) * sqr (t3 (tb)))
              - t1 (tb) * (0.5L - b43 (tb) * t3 (tb)))
    / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = (0.5L - b42 (tb) * t2 (tb) - b43 (tb) * t3 (tb)) / t1 (tb);
  if (isnan (b41 (tb)))
    return 0;
  b31 (tb) = ((1.L / 6.L - b42 (tb) * b21 (tb) * t1 (tb)) / b43 (tb)
              - b32 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (b31 (tb)))
    return 0;
  rk_b_4 (tb);
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 3rd order, 4th order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_4_3t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b32 (tb) = r[4];
  b43 (tb) = (0.25L - 1.L / 3.L * t1 (tb)
              - (1.L / 3.L - 0.5L * t1 (tb)) * t2 (tb))
    / (t3 (tb) * (t3 (tb) - t2 (tb)) * (t3 (tb) - t1 (tb)));
  b42 (tb) = ((1.L / 3.L - b43 (tb) * sqr (t3 (tb)))
              - t1 (tb) * (0.5L - b43 (tb) * t3 (tb)))
    / (t2 (tb) * (t2 (tb) - t1 (tb)));
  b41 (tb) = (0.5L - b42 (tb) * t2 (tb) - b43 (tb) * t3 (tb)) / t1 (tb);
  b31 (tb) = ((1.L / 6.L - b42 (tb) * b21 (tb) * t1 (tb)) / b43 (tb)
              - b32 (tb) * t2 (tb)) / t1 (tb);
  rk_b_4 (tb);
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3t: end\n");
#endif
  if (isnan (b31 (tb)) || isnan (b41 (tb)) || isnan (b42 (tb))
      || isnan (b43 (tb)))
    return 0;
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 2nd-3rd order Runge-Kutta
 * pair.
 */
int
rk_tb_4_3p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3p: start\n");
#endif
  if (!rk_tb_4_3 (optimize))
    return 0;
  tb = optimize->coefficient;
  e41 (tb) = 0.5L / t1 (tb);
  e42 (tb) = 0.L;
  rk_e_4 (tb);
#if DEBUG_RK_4_3
  rk_print_e (optimize, "rk_tb_4_3p", stderr);
  fprintf (stderr, "rk_tb_4_3p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 2nd-3rd order, 3rd-4th order
 * in equations depending only in time, Runge-Kutta pair.
 */
int
rk_tb_4_3tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_tb_4_3tp: start\n");
#endif
  if (!rk_tb_4_3t (optimize))
    return 0;
  tb = optimize->coefficient;
  e42 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  e41 (tb) = (0.5L - e42 (tb) * t2 (tb)) / t1 (tb);
  rk_e_4 (tb);
#if DEBUG_RK_4_3
  rk_print_e (optimize, "rk_tb_4_3tp", stderr);
  fprintf (stderr, "rk_tb_4_3tp: end\n");
#endif
  if (isnan (e42 (tb)) || isnan (e41 (tb)))
    return 0;
  return 1;
}

/**
 * Function to calculate the objective function of a 4 steps 3rd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_3 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_3: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 4 steps 3rd order, 4th
 * oder in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_3t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3t: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_3t: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 4 steps 2nd-3rd order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_3p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3p: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (e40 (tb) < 0.L)
    o += e40 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_3p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 4 steps 2nd-3rd order,
 * 3rd-4th oder in equations depending only in time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_3tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3tp: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (e40 (tb) < 0.L)
    o += e40 (tb);
  if (e41 (tb) < 0.L)
    o += e41 (tb);
  if (e42 (tb) < 0.L)
    o += e42 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_3
  fprintf (stderr, "rk_objective_tb_4_3tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_3tp: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#include <gbpLib.h>
#include <gbpMisc.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_spline.h>

void compute_triaxiality(double *x_in,
                         double *y_in,
                         double *z_in,
                         double  x_cen,
                         double  y_cen,
                         double  z_cen,
                         double  box_size,
                         int     n_particles,
                         size_t *sort_index,
                         double  return_values[3],
                         double  return_vectors[3][3]) {
    double                     s, s_new;
    double                     q, q_new;
    double                     a_new;
    double                     b_new;
    double                     c_new;
    int                        i, j;
    double                     M_tmp;
    double                     x_tmp;
    double                     y_tmp;
    double                     z_tmp;
    double *                   m_p;
    double *                   x;
    double *                   y;
    double *                   z;
    double                     inv_q2;
    double                     inv_s2;
    double                     inv_r2, r2;
    double                     M[9];
    int                        n_iterations;
    int                        continue_flag;
    double                     convergence      = 0.0001;
    int                        n_iterations_max = 200;
    gsl_vector *               eigen_vector;
    gsl_matrix *               eigen_matrix;
    gsl_eigen_symmv_workspace *w;
    gsl_matrix_view            m;

    // Initialize a bunch of stuff
    a_new                = 1.;
    b_new                = 1.;
    c_new                = 1.;
    q_new                = 1.;
    s_new                = 1.;
    return_vectors[0][0] = 1.;
    return_vectors[1][0] = 0.;
    return_vectors[2][0] = 0.;
    return_vectors[0][1] = 0.;
    return_vectors[1][1] = 1.;
    return_vectors[2][1] = 0.;
    return_vectors[0][2] = 0.;
    return_vectors[1][2] = 0.;
    return_vectors[2][2] = 1.;

    if(n_particles > 0)
        continue_flag = GBP_TRUE;
    else
        continue_flag = GBP_FALSE;
    eigen_vector = gsl_vector_alloc(3);
    eigen_matrix = gsl_matrix_alloc(3, 3);
    w            = gsl_eigen_symmv_alloc(3);

    m_p = (double *)SID_malloc(sizeof(double) * n_particles);
    x   = (double *)SID_malloc(sizeof(double) * n_particles);
    y   = (double *)SID_malloc(sizeof(double) * n_particles);
    z   = (double *)SID_malloc(sizeof(double) * n_particles);
    if(sort_index != NULL) {
        for(i = 0; i < n_particles; i++) {
            m_p[i] = 1.;
            x[i]   = d_periodic((double)x_in[sort_index[i]] - x_cen, box_size);
            y[i]   = d_periodic((double)y_in[sort_index[i]] - y_cen, box_size);
            z[i]   = d_periodic((double)z_in[sort_index[i]] - z_cen, box_size);
        }
    } else {
        for(i = 0; i < n_particles; i++) {
            m_p[i] = 1.;
            x[i]   = d_periodic((double)x_in[i] - x_cen, box_size);
            y[i]   = d_periodic((double)y_in[i] - y_cen, box_size);
            z[i]   = d_periodic((double)z_in[i] - z_cen, box_size);
        }
    }

    // Iterate until convergence
    n_iterations = 0;
    while(continue_flag) {
        q      = q_new;
        s      = s_new;
        inv_q2 = 1. / (q * q);
        inv_s2 = 1. / (s * s);

        // Construct the moment of inertia tensor
        for(i = 0; i < 9; i++)
            M[i] = 0.;
        for(i = 0; i < n_particles; i++) {
            x_tmp = x[i];
            y_tmp = y[i];
            z_tmp = z[i];
            M_tmp = m_p[i];
            r2    = x_tmp * x_tmp + y_tmp * y_tmp * inv_q2 + z_tmp * z_tmp * inv_s2;
            if(r2 > 0.) {
                inv_r2 = 1. / r2;
                M[0] += x_tmp * x_tmp * M_tmp * inv_r2;
                M[1] += y_tmp * x_tmp * M_tmp * inv_r2;
                M[2] += z_tmp * x_tmp * M_tmp * inv_r2;
                M[3] += x_tmp * y_tmp * M_tmp * inv_r2;
                M[4] += y_tmp * y_tmp * M_tmp * inv_r2;
                M[5] += z_tmp * y_tmp * M_tmp * inv_r2;
                M[6] += x_tmp * z_tmp * M_tmp * inv_r2;
                M[7] += y_tmp * z_tmp * M_tmp * inv_r2;
                M[8] += z_tmp * z_tmp * M_tmp * inv_r2;
            }
        }

        // Solve for (and sort) the eigen values and eigen vectors
        m = gsl_matrix_view_array(M, 3, 3);
        gsl_eigen_symmv(&m.matrix, eigen_vector, eigen_matrix, w);
        gsl_eigen_symmv_sort(eigen_vector, eigen_matrix, GSL_EIGEN_SORT_ABS_DESC);

        // Convert gsl vectors and such into something simpler to use
        for(i = 0; i < 3; i++) {
            return_values[i]     = sqrt(gsl_vector_get(eigen_vector, i));
            return_vectors[i][0] = gsl_matrix_get(eigen_matrix, 0, i);
            return_vectors[i][1] = gsl_matrix_get(eigen_matrix, 1, i);
            return_vectors[i][2] = gsl_matrix_get(eigen_matrix, 2, i);
        }
        q_new = return_values[1] / return_values[0];
        s_new = return_values[2] / return_values[0];

        // Check for convergence
        n_iterations++;
        if(n_iterations >= n_iterations_max)
            continue_flag = GBP_FALSE;
        if((double)fabs((float)((q_new - q) / q)) < convergence && (double)fabs((float)((s_new - s) / s)) < convergence)
            continue_flag = GBP_FALSE;
    }

    // Clean-up
    gsl_eigen_symmv_free(w);
    gsl_vector_free(eigen_vector);
    gsl_matrix_free(eigen_matrix);
    SID_free((void **)(&x));
    SID_free((void **)(&y));
    SID_free((void **)(&z));
    SID_free((void **)(&m_p));
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


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

/* Bulirsch-Stoer Implicit */

/* Author:  G. Jungman
 */
/* Bader-Deuflhard implicit extrapolative stepper.
 * [Numer. Math., 41, 373 (1983)]
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

#define SEQUENCE_COUNT 8
#define SEQUENCE_MAX   7

/* Bader-Deuflhard extrapolation sequence */
static const int bd_sequence[SEQUENCE_COUNT] =
  { 2, 6, 10, 14, 22, 34, 50, 70 };

typedef struct
{
  gsl_matrix *d;                /* workspace for extrapolation         */
  gsl_matrix *a_mat;            /* workspace for linear system matrix  */
  gsl_permutation *p_vec;       /* workspace for LU permutation        */

  double x[SEQUENCE_MAX];       /* workspace for extrapolation */

  /* state info */
  size_t k_current;
  size_t k_choice;
  double h_next;
  double eps;

  /* workspace for extrapolation step */
  double *yp;
  double *y_save;
  double *yerr_save;
  double *y_extrap_save;
  double *y_extrap_sequence;
  double *extrap_work;
  double *dfdt;
  double *y_temp;
  double *delta_temp;
  double *weight;
  gsl_matrix *dfdy;

  /* workspace for the basic stepper */
  double *rhs_temp;
  double *delta;

  /* order of last step */
  size_t order;
}
bsimp_state_t;

/* Compute weighting factor */

static void
compute_weights (const double y[], double w[], size_t dim)
{
  size_t i;
  for (i = 0; i < dim; i++)
    {
      double u = fabs(y[i]);
      w[i] = (u > 0.0) ? u : 1.0;
    }
}

/* Calculate a choice for the "order" of the method, using the
 * Deuflhard criteria.
 */

static size_t
bsimp_deuf_kchoice (double eps, size_t dimension)
{
  const double safety_f = 0.25;
  const double small_eps = safety_f * eps;

  double a_work[SEQUENCE_COUNT];
  double alpha[SEQUENCE_MAX][SEQUENCE_MAX];

  int i, k;

  a_work[0] = bd_sequence[0] + 1.0;

  for (k = 0; k < SEQUENCE_MAX; k++)
    {
      a_work[k + 1] = a_work[k] + bd_sequence[k + 1];
    }

  for (i = 0; i < SEQUENCE_MAX; i++)
    {
      alpha[i][i] = 1.0;
      for (k = 0; k < i; k++)
        {
          const double tmp1 = a_work[k + 1] - a_work[i + 1];
          const double tmp2 = (a_work[i + 1] - a_work[0] + 1.0) * (2 * k + 1);
          alpha[k][i] = pow (small_eps, tmp1 / tmp2);
        }
    }

  a_work[0] += dimension;

  for (k = 0; k < SEQUENCE_MAX; k++)
    {
      a_work[k + 1] = a_work[k] + bd_sequence[k + 1];
    }

  for (k = 0; k < SEQUENCE_MAX - 1; k++)
    {
      if (a_work[k + 2] > a_work[k + 1] * alpha[k][k + 1])
        break;
    }

  return k;
}

static void
poly_extrap (gsl_matrix * d,
             const double x[],
             const unsigned int i_step,
             const double x_i,
             const double y_i[],
             double y_0[], double y_0_err[], double work[], const size_t dim)
{
  size_t j, k;

  DBL_MEMCPY (y_0_err, y_i, dim);
  DBL_MEMCPY (y_0, y_i, dim);

  if (i_step == 0)
    {
      for (j = 0; j < dim; j++)
        {
          gsl_matrix_set (d, 0, j, y_i[j]);
        }
    }
  else
    {
      DBL_MEMCPY (work, y_i, dim);

      for (k = 0; k < i_step; k++)
        {
          double delta = 1.0 / (x[i_step - k - 1] - x_i);
          const double f1 = delta * x_i;
          const double f2 = delta * x[i_step - k - 1];

          for (j = 0; j < dim; j++)
            {
              const double q_kj = gsl_matrix_get (d, k, j);
              gsl_matrix_set (d, k, j, y_0_err[j]);
              delta = work[j] - q_kj;
              y_0_err[j] = f1 * delta;
              work[j] = f2 * delta;
              y_0[j] += y_0_err[j];
            }
        }

      for (j = 0; j < dim; j++)
        {
          gsl_matrix_set (d, i_step, j, y_0_err[j]);
        }
    }
}

/* Basic implicit Bulirsch-Stoer step.  Divide the step h_total into
 * n_step smaller steps and do the Bader-Deuflhard semi-implicit
 * iteration.  */

static int
bsimp_step_local (void *vstate,
                  size_t dim,
                  const double t0,
                  const double h_total,
                  const unsigned int n_step,
                  const double y[],
                  const double yp[],
                  const double dfdt[],
                  const gsl_matrix * dfdy,
                  double y_out[],
                  const gsl_odeiv_system * sys)
{
  bsimp_state_t *state = (bsimp_state_t *) vstate;

  gsl_matrix *const a_mat = state->a_mat;
  gsl_permutation *const p_vec = state->p_vec;

  double *const delta = state->delta;
  double *const y_temp = state->y_temp;
  double *const delta_temp = state->delta_temp;
  double *const rhs_temp = state->rhs_temp;
  double *const w = state->weight;

  gsl_vector_view y_temp_vec = gsl_vector_view_array (y_temp, dim);
  gsl_vector_view delta_temp_vec = gsl_vector_view_array (delta_temp, dim);
  gsl_vector_view rhs_temp_vec = gsl_vector_view_array (rhs_temp, dim);

  const double h = h_total / n_step;
  double t = t0 + h;

  double sum;

  /* This is the factor sigma referred to in equation 3.4 of the
     paper.  A relative change in y exceeding sigma indicates a
     runaway behavior. According to the authors suitable values for
     sigma are >>1.  I have chosen a value of 100*dim. BJG */

  const double max_sum = 100.0 * dim;

  int signum, status;
  size_t i, j;
  size_t n_inter;

  /* Calculate the matrix for the linear system. */
  for (i = 0; i < dim; i++)
    {
      for (j = 0; j < dim; j++)
        {
          gsl_matrix_set (a_mat, i, j, -h * gsl_matrix_get (dfdy, i, j));
        }
      gsl_matrix_set (a_mat, i, i, gsl_matrix_get (a_mat, i, i) + 1.0);
    }

  /* LU decomposition for the linear system. */

  gsl_linalg_LU_decomp (a_mat, p_vec, &signum);

  /* Compute weighting factors */

  compute_weights (y, w, dim);

  /* Initial step. */

  for (i = 0; i < dim; i++)
    {
      y_temp[i] = h * (yp[i] + h * dfdt[i]);
    }

  gsl_linalg_LU_solve (a_mat, p_vec, &y_temp_vec.vector, &delta_temp_vec.vector);

  sum = 0.0;

  for (i = 0; i < dim; i++)
    {
      const double di = delta_temp[i];
      delta[i] = di;
      y_temp[i] = y[i] + di;
      sum += fabs(di) / w[i];
    }

  if (sum > max_sum)
    {
      return GSL_EFAILED ;
    }

  /* Intermediate steps. */

  status = GSL_ODEIV_FN_EVAL (sys, t, y_temp, y_out);

  if (status)
    {
      return status;
    }

  for (n_inter = 1; n_inter < n_step; n_inter++)
    {
      for (i = 0; i < dim; i++)
        {
          rhs_temp[i] = h * y_out[i] - delta[i];
        }

      gsl_linalg_LU_solve (a_mat, p_vec, &rhs_temp_vec.vector, &delta_temp_vec.vector);

      sum = 0.0;

      for (i = 0; i < dim; i++)
        {
          delta[i] += 2.0 * delta_temp[i];
          y_temp[i] += delta[i];
          sum += fabs(delta[i]) / w[i];
        }

      if (sum > max_sum)
        {
          return GSL_EFAILED ;
        }

      t += h;

      status = GSL_ODEIV_FN_EVAL (sys, t, y_temp, y_out);

      if (status)
        {
          return status;
        }
    }

  /* Final step. */

  for (i = 0; i < dim; i++)
    {
      rhs_temp[i] = h * y_out[i] - delta[i];
    }

  gsl_linalg_LU_solve (a_mat, p_vec, &rhs_temp_vec.vector, &delta_temp_vec.vector);

  sum = 0.0;

  for (i = 0; i < dim; i++)
    {
      y_out[i] = y_temp[i] + delta_temp[i];
      sum += fabs(delta_temp[i]) / w[i];
    }

  if (sum > max_sum)
    {
      return GSL_EFAILED ;
    }

  return GSL_SUCCESS;
}

static void *
bsimp_alloc (size_t dim)
{
  bsimp_state_t *state = (bsimp_state_t *) malloc (sizeof (bsimp_state_t));

  state->d = gsl_matrix_alloc (SEQUENCE_MAX, dim);
  state->a_mat = gsl_matrix_alloc (dim, dim);
  state->p_vec = gsl_permutation_alloc (dim);

  state->yp = (double *) malloc (dim * sizeof (double));
  state->y_save = (double *) malloc (dim * sizeof (double));
  state->yerr_save = (double *) malloc (dim * sizeof (double));
  state->y_extrap_save = (double *) malloc (dim * sizeof (double));
  state->y_extrap_sequence = (double *) malloc (dim * sizeof (double));
  state->extrap_work = (double *) malloc (dim * sizeof (double));
  state->dfdt = (double *) malloc (dim * sizeof (double));
  state->y_temp = (double *) malloc (dim * sizeof (double));
  state->delta_temp = (double *) malloc (dim * sizeof(double));
  state->weight = (double *) malloc (dim * sizeof(double));

  state->dfdy = gsl_matrix_alloc (dim, dim);

  state->rhs_temp = (double *) malloc (dim * sizeof(double));
  state->delta = (double *) malloc (dim * sizeof (double));

  {
    size_t k_choice = bsimp_deuf_kchoice (GSL_SQRT_DBL_EPSILON, dim); /*FIXME: choice of epsilon? */
    state->k_choice = k_choice;
    state->k_current = k_choice;
    state->order = 2 * k_choice;
  }

  state->h_next = -GSL_SQRT_DBL_MAX;

  return state;
}

/* Perform the basic semi-implicit extrapolation
 * step, of size h, at a Deuflhard determined order.
 */
static int
bsimp_apply (void *vstate,
             size_t dim,
             double t,
             double h,
             double y[],
             double yerr[],
             const double dydt_in[],
             double dydt_out[],
             const gsl_odeiv_system * sys)
{
  bsimp_state_t *state = (bsimp_state_t *) vstate;

  double *const x = state->x;
  double *const yp = state->yp;
  double *const y_save = state->y_save;
  double *const yerr_save = state->yerr_save;
  double *const y_extrap_sequence = state->y_extrap_sequence;
  double *const y_extrap_save = state->y_extrap_save;
  double *const extrap_work = state->extrap_work;
  double *const dfdt = state->dfdt;
  gsl_matrix *d = state->d;
  gsl_matrix *dfdy = state->dfdy;

  const double t_local = t;
  size_t i, k;

  if (h + t_local == t_local)
    {
      return GSL_EUNDRFLW;      /* FIXME: error condition */
    }

  DBL_MEMCPY (y_extrap_save, y, dim);

  /* Save inputs */
  DBL_MEMCPY (y_save, y, dim);
  DBL_MEMCPY (yerr_save, yerr, dim);

  /* Evaluate the derivative. */
  if (dydt_in != NULL)
    {
      DBL_MEMCPY (yp, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t_local, y, yp);

      if (s != GSL_SUCCESS)
	{
          return s;
	}
    }

  /* Evaluate the Jacobian for the system. */
  {
    int s = GSL_ODEIV_JA_EVAL (sys, t_local, y, dfdy->data, dfdt);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Make a series of refined extrapolations,
   * up to the specified maximum order, which
   * was calculated based on the Deuflhard
   * criterion upon state initialization.  */
  for (k = 0; k <= state->k_current; k++)
    {
      const unsigned int N = bd_sequence[k];
      const double r = (h / N);
      const double x_k = r * r;

      int status = bsimp_step_local (state,
                                     dim, t_local, h, N,
                                     y_extrap_save, yp,
                                     dfdt, dfdy,
                                     y_extrap_sequence,
                                     sys);

      if (status == GSL_EFAILED)
        {
          /* If the local step fails, set the error to infinity in
             order to force a reduction in the step size */

          for (i = 0; i < dim; i++)
            {
              yerr[i] = GSL_POSINF;
            }

          break;
        }

      else if (status != GSL_SUCCESS)
	{
	  return status;
	}

      x[k] = x_k;

      poly_extrap (d, x, k, x_k, y_extrap_sequence, y, yerr, extrap_work, dim);
    }

  /* Evaluate dydt_out[]. */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          DBL_MEMCPY (y, y_save, dim);
          DBL_MEMCPY (yerr, yerr_save, dim);
          return s;
        }
    }

  return GSL_SUCCESS;
}

static unsigned int
bsimp_order (void *vstate)
{
  bsimp_state_t *state = (bsimp_state_t *) vstate;
  return state->order;
}

static int
bsimp_reset (void *vstate, size_t dim)
{
  bsimp_state_t *state = (bsimp_state_t *) vstate;

  state->h_next = 0;

  DBL_ZERO_MEMSET (state->yp, dim);

  return GSL_SUCCESS;
}

static void
bsimp_free (void * vstate)
{
  bsimp_state_t *state = (bsimp_state_t *) vstate;

  free (state->delta);
  free (state->rhs_temp);

  gsl_matrix_free (state->dfdy);

  free (state->weight);
  free (state->delta_temp);
  free (state->y_temp);
  free (state->dfdt);
  free (state->extrap_work);
  free (state->y_extrap_sequence);
  free (state->y_extrap_save);
  free (state->y_save);
  free (state->yerr_save);
  free (state->yp);

  gsl_permutation_free (state->p_vec);
  gsl_matrix_free (state->a_mat);
  gsl_matrix_free (state->d);
  free (state);
}

static const gsl_odeiv_step_type bsimp_type = {
  "bsimp",                      /* name */
  1,                            /* can use dydt_in */
  1,                            /* gives exact dydt_out */
  &bsimp_alloc,
  &bsimp_apply,
  &bsimp_reset,
  &bsimp_order,
  &bsimp_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_bsimp = &bsimp_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified