src/parsec/disk-image/parsec/parsec-benchmark/pkgs/libs/gsl/src/interpolation/akima.c - public/gem5-resources - Git at Google

 /* interpolation/akima.c
  *
  * Copyright (C) 1996, 1997, 1998, 1999, 2000 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.
  */

 /* Author:  G. Jungman
  */
 #include <config.h>
 #include <stdlib.h>
 #include <math.h>
 #include <gsl/gsl_errno.h>
 #include "integ_eval.h"
 #include <gsl/gsl_interp.h>

 typedef struct
 {
   double * b;
   double * c;
   double * d;
   double * _m;
 } akima_state_t;


 /* common creation */
 static void *
 akima_alloc (size_t size)
 {
   akima_state_t *state = (akima_state_t *) malloc (sizeof (akima_state_t));

   if (state == NULL)
     {
       GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
     }

   state->b = (double *) malloc (size * sizeof (double));

   if (state->b == NULL)
     {
       free (state);
       GSL_ERROR_NULL("failed to allocate space for b", GSL_ENOMEM);
     }

   state->c = (double *) malloc (size * sizeof (double));

   if (state->c == NULL)
     {
       free (state->b);
       free (state);
       GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
     }

   state->d = (double *) malloc (size * sizeof (double));

   if (state->d == NULL)
     {
       free (state->c);
       free (state->b);
       free (state);
       GSL_ERROR_NULL("failed to allocate space for d", GSL_ENOMEM);
     }

   state->_m = (double *) malloc ((size + 4) * sizeof (double));

   if (state->_m == NULL)
     {
       free (state->d);
       free (state->c);
       free (state->b);
       free (state);
       GSL_ERROR_NULL("failed to allocate space for _m", GSL_ENOMEM);
     }

   return state;
 }


 /* common calculation */
 static void
 akima_calc (const double x_array[], double b[],  double c[],  double d[], size_t size, double m[])
 {
   size_t i;

   for (i = 0; i < (size - 1); i++)
     {
       const double NE = fabs (m[i + 1] - m[i]) + fabs (m[i - 1] - m[i - 2]);
       if (NE == 0.0)
         {
           b[i] = m[i];
           c[i] = 0.0;
           d[i] = 0.0;
         }
       else
         {
           const double h_i = x_array[i + 1] - x_array[i];
           const double NE_next = fabs (m[i + 2] - m[i + 1]) + fabs (m[i] - m[i - 1]);
           const double alpha_i = fabs (m[i - 1] - m[i - 2]) / NE;
           double alpha_ip1;
           double tL_ip1;
           if (NE_next == 0.0)
             {
               tL_ip1 = m[i];
             }
           else
             {
               alpha_ip1 = fabs (m[i] - m[i - 1]) / NE_next;
               tL_ip1 = (1.0 - alpha_ip1) * m[i] + alpha_ip1 * m[i + 1];
             }
           b[i] = (1.0 - alpha_i) * m[i - 1] + alpha_i * m[i];
           c[i] = (3.0 * m[i] - 2.0 * b[i] - tL_ip1) / h_i;
           d[i] = (b[i] + tL_ip1 - 2.0 * m[i]) / (h_i * h_i);
         }
     }
 }


 static int
 akima_init (void * vstate, const double x_array[], const double y_array[],
             size_t size)
 {
   akima_state_t *state = (akima_state_t *) vstate;

   double * m = state->_m + 2; /* offset so we can address the -1,-2
                                  components */

   size_t i;
   for (i = 0; i <= size - 2; i++)
     {
       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
     }

   /* non-periodic boundary conditions */
   m[-2] = 3.0 * m[0] - 2.0 * m[1];
   m[-1] = 2.0 * m[0] - m[1];
   m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
   m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];

   akima_calc (x_array, state->b, state->c, state->d, size, m);

   return GSL_SUCCESS;
 }


 static int
 akima_init_periodic (void * vstate,
                      const double x_array[],
                      const double y_array[],
                      size_t size)
 {
   akima_state_t *state = (akima_state_t *) vstate;

   double * m = state->_m + 2; /* offset so we can address the -1,-2
                                  components */

   size_t i;
   for (i = 0; i <= size - 2; i++)
     {
       m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
     }

   /* periodic boundary conditions */
   m[-2] = m[size - 1 - 2];
   m[-1] = m[size - 1 - 1];
   m[size - 1] = m[0];
   m[size] = m[1];

   akima_calc (x_array, state->b, state->c, state->d, size, m);

   return GSL_SUCCESS;
 }

 static void
 akima_free (void * vstate)
 {
   akima_state_t *state = (akima_state_t *) vstate;

   free (state->b);
   free (state->c);
   free (state->d);
   free (state->_m);
   free (state);
 }


 static
 int
 akima_eval (const void * vstate,
             const double x_array[], const double y_array[], size_t size,
             double x,
             gsl_interp_accel * a,
             double *y)
 {
   const akima_state_t *state = (const akima_state_t *) vstate;

   size_t index;

   if (a != 0)
     {
       index = gsl_interp_accel_find (a, x_array, size, x);
     }
   else
     {
       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
     }

   /* evaluate */
   {
     const double x_lo = x_array[index];
     const double delx = x - x_lo;
     const double b = state->b[index];
     const double c = state->c[index];
     const double d = state->d[index];
     *y = y_array[index] + delx * (b + delx * (c + d * delx));
     return GSL_SUCCESS;
   }
 }


 static int
 akima_eval_deriv (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   double x,
                   gsl_interp_accel * a,
                   double *dydx)
 {
   const akima_state_t *state = (const akima_state_t *) vstate;

   size_t index;

   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

   if (a != 0)
     {
       index = gsl_interp_accel_find (a, x_array, size, x);
     }
   else
     {
       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
     }

   /* evaluate */
   {
     double x_lo = x_array[index];
     double delx = x - x_lo;
     double b = state->b[index];
     double c = state->c[index];
     double d = state->d[index];
     *dydx = b + delx * (2.0 * c + 3.0 * d * delx);
     return GSL_SUCCESS;
   }
 }


 static
 int
 akima_eval_deriv2 (const void * vstate,
                    const double x_array[], const double y_array[], size_t size,
                    double x,
                    gsl_interp_accel * a,
                    double *y_pp)
 {
   const akima_state_t *state = (const akima_state_t *) vstate;

   size_t index;

   DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

   if (a != 0)
     {
       index = gsl_interp_accel_find (a, x_array, size, x);
     }
   else
     {
       index = gsl_interp_bsearch (x_array, x, 0, size - 1);
     }

   /* evaluate */
   {
     const double x_lo = x_array[index];
     const double delx = x - x_lo;
     const double c = state->c[index];
     const double d = state->d[index];
     *y_pp = 2.0 * c + 6.0 * d * delx;
     return GSL_SUCCESS;
   }
 }


 static
 int
 akima_eval_integ (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   gsl_interp_accel * acc,
                   double a, double b,
                   double * result)
 {
   const akima_state_t *state = (const akima_state_t *) vstate;

   size_t i, index_a, index_b;

   if (acc != 0)
     {
       index_a = gsl_interp_accel_find (acc, x_array, size, a);
       index_b = gsl_interp_accel_find (acc, x_array, size, b);
     }
   else
     {
       index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
       index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
     }

   *result = 0.0;

   /* interior intervals */

   for(i=index_a; i<=index_b; i++) {
     const double x_hi = x_array[i + 1];
     const double x_lo = x_array[i];
     const double y_lo = y_array[i];
     const double dx = x_hi - x_lo;
     if(dx != 0.0) {

       if (i == index_a || i == index_b)
         {
           double x1 = (i == index_a) ? a : x_lo;
           double x2 = (i == index_b) ? b : x_hi;
           *result += integ_eval (y_lo, state->b[i], state->c[i], state->d[i],
                                  x_lo, x1, x2);
         }
       else
         {
           *result += dx * (y_lo
                            + dx*(0.5*state->b[i]
                                  + dx*(state->c[i]/3.0
                                        + 0.25*state->d[i]*dx)));
         }
     }
     else {
       *result = 0.0;
       return GSL_FAILURE;
     }
   }

   return GSL_SUCCESS;
 }


 static const gsl_interp_type akima_type =
 {
   "akima",
   5,
   &akima_alloc,
   &akima_init,
   &akima_eval,
   &akima_eval_deriv,
   &akima_eval_deriv2,
   &akima_eval_integ,
   &akima_free
 };

 const gsl_interp_type * gsl_interp_akima = &akima_type;

 static const gsl_interp_type akima_periodic_type =
 {
   "akima-periodic",
   5,
   &akima_alloc,
   &akima_init_periodic,
   &akima_eval,
   &akima_eval_deriv,
   &akima_eval_deriv2,
   &akima_eval_integ,
   &akima_free
 };

 const gsl_interp_type * gsl_interp_akima_periodic = &akima_periodic_type;
	/* interpolation/akima.c
	*
	* Copyright (C) 1996, 1997, 1998, 1999, 2000 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.
	*/

	/* Author: G. Jungman
	*/
	#include <config.h>
	#include <stdlib.h>
	#include <math.h>
	#include <gsl/gsl_errno.h>
	#include "integ_eval.h"
	#include <gsl/gsl_interp.h>

	typedef struct
	{
	double * b;
	double * c;
	double * d;
	double * _m;
	} akima_state_t;


	/* common creation */
	static void *
	akima_alloc (size_t size)
	{
	akima_state_t state = (akima_state_t ) malloc (sizeof (akima_state_t));

	if (state == NULL)
	{
	GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
	}

	state->b = (double ) malloc (size sizeof (double));

	if (state->b == NULL)
	{
	free (state);
	GSL_ERROR_NULL("failed to allocate space for b", GSL_ENOMEM);
	}

	state->c = (double ) malloc (size sizeof (double));

	if (state->c == NULL)
	{
	free (state->b);
	free (state);
	GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
	}

	state->d = (double ) malloc (size sizeof (double));

	if (state->d == NULL)
	{
	free (state->c);
	free (state->b);
	free (state);
	GSL_ERROR_NULL("failed to allocate space for d", GSL_ENOMEM);
	}

	state->_m = (double ) malloc ((size + 4) sizeof (double));

	if (state->_m == NULL)
	{
	free (state->d);
	free (state->c);
	free (state->b);
	free (state);
	GSL_ERROR_NULL("failed to allocate space for _m", GSL_ENOMEM);
	}

	return state;
	}


	/* common calculation */
	static void
	akima_calc (const double x_array[], double b[], double c[], double d[], size_t size, double m[])
	{
	size_t i;

	for (i = 0; i < (size - 1); i++)
	{
	const double NE = fabs (m[i + 1] - m[i]) + fabs (m[i - 1] - m[i - 2]);
	if (NE == 0.0)
	{
	b[i] = m[i];
	c[i] = 0.0;
	d[i] = 0.0;
	}
	else
	{
	const double h_i = x_array[i + 1] - x_array[i];
	const double NE_next = fabs (m[i + 2] - m[i + 1]) + fabs (m[i] - m[i - 1]);
	const double alpha_i = fabs (m[i - 1] - m[i - 2]) / NE;
	double alpha_ip1;
	double tL_ip1;
	if (NE_next == 0.0)
	{
	tL_ip1 = m[i];
	}
	else
	{
	alpha_ip1 = fabs (m[i] - m[i - 1]) / NE_next;
	tL_ip1 = (1.0 - alpha_ip1) * m[i] + alpha_ip1 * m[i + 1];
	}
	b[i] = (1.0 - alpha_i) * m[i - 1] + alpha_i * m[i];
	c[i] = (3.0 * m[i] - 2.0 * b[i] - tL_ip1) / h_i;
	d[i] = (b[i] + tL_ip1 - 2.0 * m[i]) / (h_i * h_i);
	}
	}
	}


	static int
	akima_init (void * vstate, const double x_array[], const double y_array[],
	size_t size)
	{
	akima_state_t state = (akima_state_t ) vstate;

	double * m = state->_m + 2; /* offset so we can address the -1,-2
	components */

	size_t i;
	for (i = 0; i <= size - 2; i++)
	{
	m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
	}

	/* non-periodic boundary conditions */
	m[-2] = 3.0 * m[0] - 2.0 * m[1];
	m[-1] = 2.0 * m[0] - m[1];
	m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
	m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];

	akima_calc (x_array, state->b, state->c, state->d, size, m);

	return GSL_SUCCESS;
	}


	static int
	akima_init_periodic (void * vstate,
	const double x_array[],
	const double y_array[],
	size_t size)
	{
	akima_state_t state = (akima_state_t ) vstate;

	double * m = state->_m + 2; /* offset so we can address the -1,-2
	components */

	size_t i;
	for (i = 0; i <= size - 2; i++)
	{
	m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
	}

	/* periodic boundary conditions */
	m[-2] = m[size - 1 - 2];
	m[-1] = m[size - 1 - 1];
	m[size - 1] = m[0];
	m[size] = m[1];

	akima_calc (x_array, state->b, state->c, state->d, size, m);

	return GSL_SUCCESS;
	}

	static void
	akima_free (void * vstate)
	{
	akima_state_t state = (akima_state_t ) vstate;

	free (state->b);
	free (state->c);
	free (state->d);
	free (state->_m);
	free (state);
	}


	static
	int
	akima_eval (const void * vstate,
	const double x_array[], const double y_array[], size_t size,
	double x,
	gsl_interp_accel * a,
	double *y)
	{
	const akima_state_t state = (const akima_state_t ) vstate;

	size_t index;

	if (a != 0)
	{
	index = gsl_interp_accel_find (a, x_array, size, x);
	}
	else
	{
	index = gsl_interp_bsearch (x_array, x, 0, size - 1);
	}

	/* evaluate */
	{
	const double x_lo = x_array[index];
	const double delx = x - x_lo;
	const double b = state->b[index];
	const double c = state->c[index];
	const double d = state->d[index];
	y = y_array[index] + delx (b + delx * (c + d * delx));
	return GSL_SUCCESS;
	}
	}


	static int
	akima_eval_deriv (const void * vstate,
	const double x_array[], const double y_array[], size_t size,
	double x,
	gsl_interp_accel * a,
	double *dydx)
	{
	const akima_state_t state = (const akima_state_t ) vstate;

	size_t index;

	DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

	if (a != 0)
	{
	index = gsl_interp_accel_find (a, x_array, size, x);
	}
	else
	{
	index = gsl_interp_bsearch (x_array, x, 0, size - 1);
	}

	/* evaluate */
	{
	double x_lo = x_array[index];
	double delx = x - x_lo;
	double b = state->b[index];
	double c = state->c[index];
	double d = state->d[index];
	dydx = b + delx (2.0 * c + 3.0 * d * delx);
	return GSL_SUCCESS;
	}
	}


	static
	int
	akima_eval_deriv2 (const void * vstate,
	const double x_array[], const double y_array[], size_t size,
	double x,
	gsl_interp_accel * a,
	double *y_pp)
	{
	const akima_state_t state = (const akima_state_t ) vstate;

	size_t index;

	DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

	if (a != 0)
	{
	index = gsl_interp_accel_find (a, x_array, size, x);
	}
	else
	{
	index = gsl_interp_bsearch (x_array, x, 0, size - 1);
	}

	/* evaluate */
	{
	const double x_lo = x_array[index];
	const double delx = x - x_lo;
	const double c = state->c[index];
	const double d = state->d[index];
	y_pp = 2.0 c + 6.0 * d * delx;
	return GSL_SUCCESS;
	}
	}


	static
	int
	akima_eval_integ (const void * vstate,
	const double x_array[], const double y_array[], size_t size,
	gsl_interp_accel * acc,
	double a, double b,
	double * result)
	{
	const akima_state_t state = (const akima_state_t ) vstate;

	size_t i, index_a, index_b;

	if (acc != 0)
	{
	index_a = gsl_interp_accel_find (acc, x_array, size, a);
	index_b = gsl_interp_accel_find (acc, x_array, size, b);
	}
	else
	{
	index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
	index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
	}

	*result = 0.0;

	/* interior intervals */

	for(i=index_a; i<=index_b; i++) {
	const double x_hi = x_array[i + 1];
	const double x_lo = x_array[i];
	const double y_lo = y_array[i];
	const double dx = x_hi - x_lo;
	if(dx != 0.0) {

	if (i == index_a \|\| i == index_b)
	{
	double x1 = (i == index_a) ? a : x_lo;
	double x2 = (i == index_b) ? b : x_hi;
	*result += integ_eval (y_lo, state->b[i], state->c[i], state->d[i],
	x_lo, x1, x2);
	}
	else
	{
	result += dx (y_lo
	+ dx(0.5state->b[i]
	+ dx*(state->c[i]/3.0
	+ 0.25state->d[i]dx)));
	}
	}
	else {
	*result = 0.0;
	return GSL_FAILURE;
	}
	}

	return GSL_SUCCESS;
	}


	static const gsl_interp_type akima_type =
	{
	"akima",
	5,
	&akima_alloc,
	&akima_init,
	&akima_eval,
	&akima_eval_deriv,
	&akima_eval_deriv2,
	&akima_eval_integ,
	&akima_free
	};

	const gsl_interp_type * gsl_interp_akima = &akima_type;

	static const gsl_interp_type akima_periodic_type =
	{
	"akima-periodic",
	5,
	&akima_alloc,
	&akima_init_periodic,
	&akima_eval,
	&akima_eval_deriv,
	&akima_eval_deriv2,
	&akima_eval_integ,
	&akima_free
	};

	const gsl_interp_type * gsl_interp_akima_periodic = &akima_periodic_type;