303 lines
8.2 KiB
C
303 lines
8.2 KiB
C
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <assert.h>
|
|
#include "myc_fortran.h"
|
|
|
|
void F_FUN(dat_axis_range)(float *xin, int *nin, float *xstart, float *xend, float *pMargin, int *i1, int *i2) {
|
|
int i;
|
|
float x1, x2, margin;
|
|
|
|
margin=*pMargin;
|
|
if (margin<0) margin=-margin;
|
|
if (*xstart>*xend) {
|
|
x1=*xend-margin;
|
|
x2=*xstart+margin;
|
|
} else {
|
|
x1=*xstart-margin;
|
|
x2=*xend+margin;
|
|
}
|
|
*i1=*nin;
|
|
*i2=0;
|
|
for (i=0; i<*nin; i++) {
|
|
if (x1 < xin[i] && xin[i] < x2) {
|
|
if (i>*i2) *i2=i;
|
|
if (i<*i1) *i1=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
void dat_fillin(float *src, float *dst, int size, float *wsum, int j, float weight) {
|
|
int k;
|
|
|
|
dst += j*size;
|
|
if (wsum!=NULL) wsum[j] += weight;
|
|
for (k=0; k<size; k++) {
|
|
*dst += weight * *src;
|
|
dst++; src++;
|
|
}
|
|
}
|
|
|
|
void dat_compress(float *xin, float *src, int nin, int grpsize, int nblocks
|
|
, float xstart, float xs, int nout, float *dst, float *wout) {
|
|
/*
|
|
** compress (sum up) nin input channels to nout output channels
|
|
** source and destination may be multidimensional arrays, but the compression
|
|
** is performed only in one dimension (d). The shape of the array is summarized
|
|
** by the following values:
|
|
** grpsize product of the size of all dimensions lower (d)
|
|
** nblocks product of the size of all dimensions higher (d)
|
|
** nin length of compressed dimension
|
|
** Examples:
|
|
** number of dimensions 1 2 3
|
|
** x-range 0..9 0..9 0..4
|
|
** y-range 0..2 0..6
|
|
** z-range 0..8
|
|
** index formula x x+y*10 x+y*(5+z*7)
|
|
** dimension to compress x x y
|
|
** grpsize 1 1 5
|
|
** nblocks 1 3 9
|
|
** nin 10 10 7
|
|
**
|
|
** geometric principe of compression:
|
|
** x[0] x[1] x[2] x[3] x[4] x[5]
|
|
** | | | | | |
|
|
** |-----0-----|-----1-----|-----2---|---3--|--4--|--5--|
|
|
** |----0----|----1----|----2----|----3----|
|
|
** | | | |
|
|
** x-start |<--xs--->|
|
|
**
|
|
** The values from the source are distributed to the destination according to the overlap
|
|
** of the intervals. Output channel 0 will get a small part from input channel 0 and
|
|
** a large part from channel 1, output channel 5 will get all of channel 4 plus
|
|
** some part of channel 3 and channel 5.
|
|
**
|
|
** xin x-values (size nin) (may be NULL, if x-values are [0.0,1.0,2.0,3.0,4.0,....])
|
|
** src source (size grpsize*nin*nblocks)
|
|
** xstart start value
|
|
** xs x-step for output
|
|
** nout number of output channels
|
|
** dst destination (size grpsize*nout*nblocks)
|
|
** wout output weights (size nout, may be NULL if not needed)
|
|
*/
|
|
int i, j, m, ja, jb, bsize;
|
|
float w, x, xl, xr, xa, xb, dx0, dx1, dx;
|
|
|
|
assert(nin>=1 && src !=NULL && grpsize>0 && nblocks>0 && xs!=0.0 && nout>0 && dst!=NULL && wout!=NULL);
|
|
for (j=0; j<nout; j++) {
|
|
wout[j]=0;
|
|
}
|
|
bsize=grpsize*nout;
|
|
for (j=0; j<bsize*nblocks; j++) {
|
|
dst[j]=0;
|
|
}
|
|
for (m=0; m<nblocks; m++) {
|
|
if (xin==NULL) {
|
|
dx=1.0/xs;
|
|
dx0=dx/2;
|
|
} else {
|
|
dx0=(xin[1]-xin[0])/xs*0.5;
|
|
}
|
|
if (nout==1) {
|
|
for (i=0; i<nin; i++) {
|
|
if (xin==NULL) {
|
|
xa=(i-xstart)/xs-dx0;
|
|
} else {
|
|
if (i<nin-1) dx1=(xin[i+1]-xin[i])/xs*0.5;
|
|
dx=dx0+dx1;
|
|
xa=(xin[i]-xstart)/xs-dx0;
|
|
}
|
|
xb=xa+dx;
|
|
if (xa>-0.5) {
|
|
if (xb<0.5) {
|
|
w=1.0;
|
|
} else {
|
|
w=(0.5-xa)/dx;
|
|
}
|
|
} else {
|
|
if (xb<0.5) {
|
|
w=(xb+0.5)/dx;
|
|
} else {
|
|
w=1.0;
|
|
}
|
|
}
|
|
if (w>0.0) {
|
|
dat_fillin(src, dst, grpsize, wout, 0, w);
|
|
}
|
|
src+=grpsize;
|
|
}
|
|
} else {
|
|
for (i=0; i<nin; i++) {
|
|
if (xin==NULL) {
|
|
xa=(i-xstart)/xs-dx0;
|
|
} else {
|
|
if (i<nin-1) dx1=(xin[i+1]-xin[i])/xs*0.5;
|
|
dx=dx0+dx1;
|
|
xa=(xin[i]-xstart)/xs-dx0;
|
|
}
|
|
ja=xa; if (ja>xa) ja--;
|
|
xb=xa+dx;
|
|
jb=xb; if (jb<xb) jb++;
|
|
if (ja<0) {
|
|
ja=-1;
|
|
if (xa<-1.0) xa=-1.0;
|
|
}
|
|
if (jb>=nout) {
|
|
jb=nout;
|
|
if (xb>nout) xb=nout;
|
|
}
|
|
if (ja<jb) {
|
|
if (ja+1==jb) {
|
|
x=(xa+xb)*0.5-ja;
|
|
w=(xb-xa)/dx;
|
|
if (ja>=0) dat_fillin(src, dst, grpsize, wout, ja, (1-x)*w);
|
|
if (jb<nout) dat_fillin(src, dst, grpsize, wout, jb, x*w);
|
|
} else {
|
|
xl=ja+1-xa;
|
|
xr=xb+1-jb;
|
|
w=1.0/dx;
|
|
if (ja>=0) dat_fillin(src, dst, grpsize, wout, ja, xl*xl*w/2);
|
|
if (ja+2==jb) {
|
|
dat_fillin(src, dst, grpsize, wout, ja+1, ((1-xl/2)*xl+(1-xr/2)*xr)*w);
|
|
} else {
|
|
dat_fillin(src, dst, grpsize, wout, ja+1, ((1-xl/2)*xl+0.5)*w);
|
|
for (j=ja+1; j<jb-1; j++) {
|
|
dat_fillin(src, dst, grpsize, wout, j, w);
|
|
}
|
|
dat_fillin(src, dst, grpsize, wout, jb-1, ((1-xr/2)*xr+0.5)*w);
|
|
}
|
|
if (jb<nout) dat_fillin(src, dst, grpsize, wout, jb, xr*xr*w/2);
|
|
}
|
|
}
|
|
src+=grpsize;
|
|
}
|
|
}
|
|
wout=NULL; /* wout is stored only in the first block */
|
|
dst+=bsize;
|
|
}
|
|
}
|
|
|
|
int F_FUN(dat_comp)(int *flgs, float *xin, float **fin, int *nin, int *grpsize, int *nblocks
|
|
, float *xstart, float *xs, int *nout, float **weights) {
|
|
/* flgs 0,2: xin unused 1,3: xin used
|
|
** flgs 0,1: weights unused 2,3: weights used
|
|
*/
|
|
float *fout, *wout;
|
|
int l;
|
|
|
|
wout=calloc(*nout, sizeof(float));
|
|
if (wout==NULL) return 0;
|
|
l = *nout * *grpsize * *nblocks;
|
|
fout=calloc(l, sizeof(float));
|
|
if (fout==NULL) { free(wout); return 0; }
|
|
if (0==(*flgs & 1)) xin=NULL;
|
|
dat_compress(xin, *fin, *nin, *grpsize, *nblocks, *xstart, *xs, *nout, fout, wout);
|
|
free(*fin);
|
|
*fin=fout;
|
|
if (0==(*flgs &2) || weights==NULL) {
|
|
free(wout);
|
|
} else {
|
|
*weights=wout;
|
|
}
|
|
return l;
|
|
}
|
|
|
|
void F_FUN(dat_copy2f)(float **fin, float *fout, int *size, int *freeIt) {
|
|
float *p;
|
|
int i, s;
|
|
|
|
p=*fin;
|
|
s=*size;
|
|
for (i=0; i<s; i++) {
|
|
*fout++ = *p++;
|
|
}
|
|
if (*freeIt) {
|
|
free(*fin);
|
|
*fin=NULL;
|
|
}
|
|
}
|
|
|
|
void F_FUN(dat_copy2p)(float **fout, float *fin) {
|
|
float *p;
|
|
int i, s;
|
|
|
|
p=calloc(1,sizeof(float));
|
|
if (p==NULL) return;
|
|
*p=*fin;
|
|
*fout = p;
|
|
return;
|
|
}
|
|
|
|
void F_FUN(dat_extract)(float **fin, int *offset, int *step, float *fout, int *size, int *freeIt) {
|
|
float *p;
|
|
int i, s;
|
|
|
|
p=*fin + *offset;
|
|
s=*size;
|
|
for (i=0; i<s; i++) {
|
|
*fout++ = *p;
|
|
p += *step;
|
|
}
|
|
if (*freeIt) {
|
|
free(*fin);
|
|
*fin=NULL;
|
|
}
|
|
}
|
|
|
|
#include "napi.h"
|
|
|
|
int F_FUN(dat_nexus_getslab)(NXhandle handle, int *pType, int *rank,
|
|
int *start, int *size, float **fout) {
|
|
|
|
int status, i, type;
|
|
int cellsize, ncells;
|
|
void *data;
|
|
float *flts;
|
|
double *dptr;
|
|
int *ip;
|
|
int cstart[32], csize[32];
|
|
|
|
type=*pType;
|
|
switch (type) {
|
|
case NX_FLOAT64: cellsize=8; break;
|
|
case NX_UINT32: type=NX_INT32; /* fall through */
|
|
case NX_INT32: /* fall through */
|
|
case NX_FLOAT32: cellsize=4; break;
|
|
case NX_UINT16: type=NX_INT16; /* fall through */
|
|
case NX_INT16: cellsize=2; break;
|
|
case NX_UINT8: type=NX_INT8; /* fall through */
|
|
case NX_INT8: /* fall through */
|
|
case NX_CHAR: cellsize=1; break;
|
|
}
|
|
if (cellsize<4 || *rank<1 || *rank>32) return NX_ERROR; /* for now, support only f32, i32 and f64 */
|
|
ncells=1;
|
|
for (i=0; i<*rank; i++) {
|
|
cstart[i]=start[*rank-i-1];
|
|
csize[i]=size[*rank-i-1];
|
|
ncells=ncells*size[i];
|
|
}
|
|
for (;i<32;i++) {
|
|
cstart[i]=0;
|
|
csize[i]=1;
|
|
}
|
|
data=calloc(ncells, cellsize);
|
|
if (data==NULL) return NX_ERROR;
|
|
flts=(float *)data;
|
|
|
|
status=NXgetslab(handle, data, cstart, csize);
|
|
if (status!=NX_OK) { free(data); return NX_ERROR; }
|
|
|
|
if (type==NX_INT32) { /* convert int32 to float 32 */
|
|
ip=(int *)data;
|
|
for (i=0; i<ncells; i++) {
|
|
flts[i]=*ip; ip++;
|
|
}
|
|
} else if (type==NX_FLOAT64) { /* convert float 64 to float 32 */
|
|
dptr=(double *)data;
|
|
for (i=0; i<ncells; i++) {
|
|
flts[i]=*dptr; dptr++;
|
|
}
|
|
}
|
|
*fout=flts;
|
|
return NX_OK;
|
|
}
|