arrayh4.c

/* Copyright (c) 1999-2023 Massachusetts Institute of Technology
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 * 
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#include <stdio.h>
#include <stdlib.h>

#include "config.h"
#include "arrayh4.h"

#ifdef VERBOSE
#  define WHEN_VERBOSE(x) x
#else
#  define WHEN_VERBOSE(x) 
#endif

int arrayh4_create(arrayh4 *b, int32 numtype, intn rank, const int32 *dims)
{
     int i;
     
     b->rank = rank;
     b->N = 1;
     for (i = 0; i < rank; ++i) {
	  b->dims[i] = dims[i];
	  b->N *= dims[i];
	  if (numtype == DFNT_FLOAT64)
	       b->scale.d[i] = NULL;
	  else
	       b->scale.f[i] = NULL;
     }
     b->numtype = numtype;
     if (numtype == DFNT_FLOAT64) {
	  b->p.d = (float64 *) malloc(b->N * sizeof(float64));
	  if (!b->p.d)
	       return 0;
     } else if (numtype == DFNT_FLOAT32) {
	  b->p.f = (float32 *) malloc(b->N * sizeof(float32));
	  if (!b->p.f)
	       return 0;
     }
     return 1;
}

int arrayh4_clone(arrayh4 *b, arrayh4 a)
{
     int dim, i;

     i = arrayh4_create(b, a.numtype, a.rank, a.dims);

     if (!i)
	  return 0;

     /* copy the axes scales: */

     if (a.numtype == DFNT_FLOAT64) {
	  for (dim = 0; dim < a.rank; ++dim)
	       if (a.scale.d[dim]) {
		    b->scale.d[dim] = (float64 *) malloc(a.dims[dim] * 
							 sizeof(float64));
		    if (!b->scale.d[dim])
			 return 0;
		    for (i = 0; i < a.dims[dim]; ++i)
			 b->scale.d[dim][i] = a.scale.d[dim][i];
	       }
     }
     else if (a.numtype == DFNT_FLOAT32) {
	  for (dim = 0; dim < a.rank; ++dim)
	       if (a.scale.f[dim]) {
		    b->scale.f[dim] = (float32 *) malloc(a.dims[dim] * 
							 sizeof(float32));
		    if (!b->scale.f[dim])
			 return 0;
		    for (i = 0; i < a.dims[dim]; ++i)
			 b->scale.f[dim][i] = a.scale.f[dim][i];
	       }
     }

     return 1;
}

void arrayh4_destroy(arrayh4 a)
{
     int i;

     if (a.numtype == DFNT_FLOAT64) {
	  free(a.p.d);
	  for (i = 0; i < a.rank; ++i)
	       free(a.scale.d[i]);
     }
     else if (a.numtype == DFNT_FLOAT32) {
	  free(a.p.f);
	  for (i = 0; i < a.rank; ++i)
	       free(a.scale.f[i]);
     }
}

int arrayh4_read(char *fname, arrayh4 *a, int require_rank)
{
     int32 numtype;
     intn rank;
     int32 dims[ARRAYH4_MAX_RANK];
     int dim;
     
     WHEN_VERBOSE(printf("Reading arrayh4 from \"%s\"...\n", fname));
     
     if (require_rank < 0 || require_rank > ARRAYH4_MAX_RANK)
	  return 0;
     
     DFSDclear();
     DFSDgetdims(fname, &rank, dims, ARRAYH4_MAX_RANK);
     
     WHEN_VERBOSE(printf("   rank = %d", rank));
     
     if (require_rank && require_rank != rank)
	  return 0;
     
     WHEN_VERBOSE(printf(", dimensions are "));
     
     for (dim = 0; dim < rank; ++dim) {
	  WHEN_VERBOSE(printf("%s%d", dim ? "x" : "", (int) dims[dim]));
     }
     
     WHEN_VERBOSE(printf("\n"));
     
     DFSDgetNT(&numtype);
     
     arrayh4_create(a, numtype, rank, dims);
     
     if (a->numtype == DFNT_FLOAT64) {
	  WHEN_VERBOSE(printf("   double precision\n"));
	  DFSDgetdata(fname, a->rank, a->dims, (VOIDP) a->p.d);
	  
	  for (dim = 0; dim < a->rank; ++dim) {
	       a->scale.d[dim] = (float64 *) malloc(a->dims[dim] * 
						    sizeof(float64));
	       if (FAIL == DFSDgetdimscale(dim+1, a->dims[dim], 
					   (VOIDP) a->scale.d[dim])) {
		    free(a->scale.d[dim]);
		    a->scale.d[dim] = NULL;
	       }
	  }
     } else if (a->numtype == DFNT_FLOAT32) {
	  WHEN_VERBOSE(printf("   single precision\n"));
	  DFSDgetdata(fname, a->rank, a->dims, (VOIDP) a->p.f);
	  
	for (dim = 0; dim < a->rank; ++dim) {
	     a->scale.f[dim] = (float32 *) malloc(a->dims[dim] * 
						  sizeof(float32));
	     if (FAIL == DFSDgetdimscale(dim+1, a->dims[dim], 
					 (VOIDP) a->scale.f[dim])) {
		  free(a->scale.f[dim]);
		  a->scale.f[dim] = NULL;
	     }
	}
     } else
	  return 0;
     
     return 1;
}

#define MIN2(a,b) ((a) < (b) ? (a) : (b))
#define MAX2(a,b) ((a) > (b) ? (a) : (b))

int arrayh4_write(char *fname, arrayh4 a)
{
     int i;
     WHEN_VERBOSE(printf("Writing arrayh4 to \"%s\"...\n", fname));
     
     remove(fname);
     
     DFSDclear();
     DFSDsetdims(a.rank, a.dims);
     DFSDsetNT(a.numtype);
     
     if (a.numtype == DFNT_FLOAT64) {
	  float64 minval = 1e40, maxval = -1e40;
	  
	  for (i = 0; i < a.N; ++i) {
	       minval = MIN2(minval, a.p.d[i]);
	       maxval = MAX2(maxval, a.p.d[i]);
	  }
	  DFSDsetrange(&maxval, &minval);
	  DFSDadddata(fname, a.rank, a.dims, (VOIDP) a.p.d);
	  
	  for (i = 0; i < a.rank; ++i)
	       if (a.scale.d[i])
		    DFSDgetdimscale(i+1, a.dims[i], (VOIDP) a.scale.d[i]);
     } else if (a.numtype == DFNT_FLOAT32) {
	  float32 minval = 1e20, maxval = -1e20;
	  
	  for (i = 0; i < a.N; ++i) {
	       minval = MIN2(minval, a.p.f[i]);
	       maxval = MAX2(maxval, a.p.f[i]);
	  }
	  DFSDsetrange(&maxval, &minval);
	  DFSDadddata(fname, a.rank, a.dims, (VOIDP) a.p.f);
	  
	  for (i = 0; i < a.rank; ++i)
	       if (a.scale.f[i])
		    DFSDgetdimscale(i+1, a.dims[i], (VOIDP) a.scale.f[i]);
     } else
	  return 0;
     
     return 1;
}

short arrayh4_conformant(arrayh4 a, arrayh4 b)
{
     int dim;
     
     if (a.numtype != b.numtype || a.N != b.N || a.rank != b.rank)
	  return 0;
     for (dim = 0; dim < a.rank; ++dim)
	  if (a.dims[dim] != b.dims[dim])
	       return 0;
     return 1;
}