/******************************************************************************
* MEMTOOLS.C                                                                  *
*                                                                             *
* Wrapper function for DAA to allocate arbitrary length arrays of memory.     *
* Thanks and acknowledgements to Richard Hogaboom for original source code    *
* and for releasing the concept as freeware.                                  *
*                                                                             *
* James Holliday                                                              *
* University of California - Davis                                            *
*                                                                             *
******************************************************************************/

#include <stddef.h>
#include <stdlib.h>

#define MAX_DIM 10

/* DAA function prototypes ------------------------------------------------- */
int   das(int,int,int*);
void* daa(int,int,int*,int*,char*,char*);

/* ALLOCATE (wrapper for DAA) ---------------------------------------------- */
void* allocate(int size, int N, int* dims, char* free)
{
  /* Declare pointer to memory space */
  void *tensor = NULL;
  char *p_free = NULL;

  /* DAA variables */
  int buffer = 0;
  int* st = NULL;

  /* Initialize array indices to begin at 0 (by using calloc) */
  st = calloc(N , sizeof(int));
  if ( st == NULL ) return NULL;

  /* Calculate the size of the packed memory block */
  buffer = das(size, N, dims);

  /* Use calloc to allocate and initialize elements to zero */
  p_free = (char *)calloc(1 , buffer);
  if ( p_free == NULL ) return NULL;

  /* Repack and reorder the elements */
  tensor = daa(size, N, dims, st, p_free, NULL);

  /* Return pointer to our tensor */
  return tensor;
}

/******************************************************************************
* Internal DAA source code (simplified)                                       *
******************************************************************************/
static int doff(int level,int *dim_ind,int dim_prod,int *dp)
{
  if ( level == 0 ) return 0;
  else return ( doff(level-1, dim_ind, dim_prod, dp) + (dim_prod/dp[level-1]) *
	 dim_ind[level-1] );
}

static int off(int level,int *dp)
{
  if ( level == 0 )  return 0;
  else return ( dp[level-1] + off(level-1, dp) );
}

static char* ptr_init(int level,int *dim_ind,int data_size,int num_dim,
		      char *data_ptr,char *ptr_ptr,int *dim,int *st,int *dp)
{
  char **ptrs;
  unsigned long i;

  if ( level < (num_dim - 1) )
  {
    ptrs = (char **) (ptr_ptr + off(level, dp) * sizeof(char *) +
	    ((level == 0)?0:doff(level, dim_ind, dp[level], dp) *
	     sizeof(char *)));

    for ( i=0 ; i<dim[level] ; i++)
    {
      dim_ind[level] = i;

      ptrs[i] = ptr_init(level+1, dim_ind, data_size, num_dim, data_ptr,
			 ptr_ptr, dim, st, dp) -
	st[level+1] * ((level+1 == num_dim-1)?data_size:sizeof(char *));
    }

      if ( level == 0 )  ptrs -= st[0];
  }
  else
  {
    ptrs = (char **) (data_ptr + doff(level, dim_ind, dp[level], dp) *
		      data_size - ((level == 0)?(st[0] * data_size):0));
  }

  return (char *) ptrs;
}

int das(int data_size,int num_dim,int *dim)
{
  int i;
  int dp[MAX_DIM];

  if ( num_dim < 1 || num_dim > MAX_DIM )  return -1;
  if ( data_size < 1 )  return -1;

  dp[0] = dim[0];
  for ( i=0 ; i<num_dim ; i++)
  {
    if ( dim[i] <= 0 )  return -1;
    if ( i > 0 )  dp[i] = dp[i-1] * dim[i];
  }

  return dp[num_dim-1] * data_size + off(num_dim-1, dp) * sizeof(char *) +
    sizeof(char *);
}

void* daa(int data_size,int num_dim,int *dim,int *st,char *base_ptr,
	  char *init_ptr)
{
  int i, j;
  int dim_ind[MAX_DIM];
  char *p_data,*p;
  int dp[MAX_DIM];
  char *ptr_ptr;

  if ( num_dim < 1 || num_dim > MAX_DIM )  return NULL;
  if ( data_size < 1 )  return NULL;

  dp[0] = dim[0];
  for ( i=0 ; i<num_dim ; i++)
  {
    dim_ind[i] = 0;
    if ( dim[i] <= 0 )  return NULL;
    if ( i > 0 )  dp[i] = dp[i-1] * dim[i];
  }

  ptr_ptr = base_ptr + dp[num_dim-1] * data_size;
  for ( i=0 ; i<sizeof(char *) ; i++, ++ptr_ptr )
    if ( ((int)ptr_ptr)%sizeof(char *) == 0 )  break;

  if ( init_ptr != NULL )
  {
    p_data = base_ptr;

    for ( i=0 ; i<dp[num_dim-1] ; i++ )
    {
      p = init_ptr;
      for ( j=0 ; j<data_size ; j++)  *p_data++ = *p++;
    }
  }

   return ptr_init(0,dim_ind,data_size,num_dim,base_ptr,ptr_ptr,dim,st,dp);
}

/* Clean up shop ----------------------------------------------------------- */
#undef MAX_DIM