LUT.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#include "LUT.h"
#include "basics.h"
#include "MULTI_ION_IF.h"
 
#include <string.h>
#include <assert.h>
#include <cmath>
 
namespace limpet {
 
using ::opencarp::dupstr;
using ::opencarp::f_open;
using ::opencarp::FILE_SPEC;
using ::opencarp::is_big_endian;
using ::opencarp::log_msg;
 
/*  create and initialize LUT
 *
 * \param p pointer to lookup table structure
 * \param rows  number of lookup variables
 * \param mn  lower boundary of LUT
 * \param mx  upper boundary of LUT
 * \param res resolution of LUT
 * \param name  table name
 *
 */
void LUT_alloc(LUT *plut,int cols,float mn,float mx,float res,const char *name, Target target)
{
  plut->name = dupstr(name);
 
  // table bounds
  plut->cols   = cols;
  plut->mn     = mn;
  plut->mx     = mx;
  plut->res    = res;
  plut->step   = 1/res;
 
  // in terms of indices
  plut->mn_ind = (int)(mn * plut->step);
  plut->mx_ind = (int)(mx * plut->step);
  plut->rows   = plut->mx_ind - plut->mn_ind + 1;
 
  // build a table and offset it
  plut->tab  = build_matrix_ns<LUT_data_t>( plut->rows, plut->cols, sizeof(LUT_data_t), target);
  plut->tab -= plut->mn_ind;
}
 
 
/*  dump LUT to file
 *
 * \param p pointer to lookup table structure
 * \param fname filename to store LUT dump
 *
 * \return 0 if successull, -1 otherwise
 */
int LUT_dump( LUT *plut, const char *fname )
{
  FILE_SPEC fdump = f_open(fname, "wb");
 
  if (!fdump)
    return -1;
 
  // write endianness and table dimensions first
  int endian_flg = is_big_endian();
  fwrite(&endian_flg, sizeof(int), 1, fdump->fd);
  fwrite(&plut->rows, sizeof(plut->rows), 1, fdump->fd);
  fwrite(&plut->cols, sizeof(plut->cols), 1, fdump->fd);
 
  for (int i=plut->mn_ind; i<=plut->mx_ind; i++)
    for (int j=0; j<plut->cols; j++)
      fwrite(&plut->tab[i][j], sizeof(LUT_data_t), 1, fdump->fd);
 
  f_close(fdump);
 
  return 0;
}
 
 
void destroy_lut( LUT *plut, Target target )
{
  if ( plut != NULL && plut->tab != NULL )  {
    if (plut->tab != NULL) {
      plut->tab += plut->mn_ind;
      deallocate_on_target<LUT_data_t>(target, plut->tab[0]);
      deallocate_on_target(target, plut->tab);
    }
    free(plut->name);
  }
  if( plut )
    memset( plut, 0, sizeof(LUT) );
}
 
 
 
 
#define MAX_LUT_NUMINF 10       
 
int check_LUT( LUT* lut )
{
  int retval = 0;
  int numinf = 0;
 
  for ( int i=lut->mn_ind; i<=lut->mx_ind; i++ )
    for ( int j=0; j<lut->cols; j++ )
      if ( ! std::isfinite(lut->tab[i][j])) {
        log_msg(0, 2, 0, "LUT WARNING: %s=%g produces %g in entry number %d!\n",lut->name,i*lut->res,lut->tab[i][j],j);
        retval = 1;
    if( numinf++ > MAX_LUT_NUMINF ) {
          log_msg(0, 3, 0, "suppressing further errors!!!\n" );
      return retval;
    }
      }
 
  return retval;
}
 
void IIF_warn(const int wv, const char error[]) {
  // log_msg(_nc_logf, 2, 0, "danger [tm = %f]: %s : local node = %d, global node = %d\n",
  //         current_global_time(), error, wv, current_global_node(wv) ); //FIXME
}
 
void LUT_problem( LUT *lt, double val, int wv, const char *tabname )
{
  // struct exception_type except =  { .exit=0 };
  char error[5000]; //FIXME, replace me with a C++ string.
 
  if ( std::isfinite(val) ) {
#ifdef CHATTY_LUT
    snprintf(error, sizeof error, "bounds exceeded for %s-table = %g (limits: %g - %g)",
            tabname, val, lt->mn, lt->mx);
    IIF_warn(wv, error);
#endif
    return;
  } else if (std::isinf(val)) {
    snprintf(error, sizeof error, "inf passed to LUT_index() for %s-table", tabname);
    IIF_warn(wv, error);
  } else if ( std::isnan(val) ) {
    snprintf(error, sizeof error, "NaN passed to LUT_index() for %s-table", tabname);
    IIF_warn(wv, error);
  }
}
 
#ifndef IMP_FAST
int LUT_index( LUT *tab, GlobalData_t val, int locind )
{
  int indx =  (int)(tab->step*val);
 
  if (indx < tab->mn_ind) {
    indx = tab->mn_ind;
  } else if (indx > tab->mx_ind) {
    indx = tab->mx_ind;
  } else {
    return indx;
  }
 
  LUT_problem(tab, val, locind, tab->name );
  return indx;
}
#endif
 
int LUT_out_of_bounds(LUT *tab, GlobalData_t val) {
  return (val < tab->mn || val > tab->mx) ;
}
 
inline LUT_data_t LUT_interp( LUT *t, int i, int j, GlobalData_t x )
{
  return (1.-x)*t->tab[i][j] + x*t->tab[i+1][j];
}
 
 
inline LUT_data_t LUT_derror( LUT *t, int idx, GlobalData_t x )
{
  return  (x-idx*t->res)/t->res;
}
 
LUT_data_t LUT_interpRow(LUT *const tab,GlobalData_t val,int i,LUT_data_t* row)
{
  int    idx  = LUT_index(tab, val, i);
  GlobalData_t derr = LUT_derror(tab, idx, val);
  if (LUT_out_of_bounds(tab, val)) {
    for (int j=0;j<tab->cols;j++)
      row[j] = tab->tab[idx][j];
  } else {
    for (int j=0;j<tab->cols;j++)
      row[j] = LUT_interp(tab, idx, j, derr );
  }
  return derr;
}
 
 
LUT_data_t*
LUT_row( LUT *lut, GlobalData_t val, int locind )
{
  return lut->tab[LUT_index( lut, val, locind )];
}
 
#ifdef __cplusplus
extern "C"
{
#endif  // ifdef __cplusplus
 
 
void LUT_interpRow_n_elements(char *table, char *val_ptr, int offset,
                              int distance, int index, int n, char* row_ptr,
                              int lut_numelements) {
  if (n <= 0)
    return;
  LUT *tab = (LUT *)table;
  GlobalData_t *val = (GlobalData_t *)(val_ptr + offset);
  LUT_data_t* row = (LUT_data_t*)row_ptr;
  if (distance == 1) {
    for (unsigned i = 0; i < n; ++i)
      LUT_interpRow(tab, *(val + i), i + index*n,
                    row + (i * (lut_numelements)));
  } else {
    for (unsigned i = 0; i < n; ++i)
      LUT_interpRow(tab, *(val + i*distance + index*n), i + index*n,
                    row + (i * (lut_numelements)));
  }
}
 
#ifdef HAS_MLIR_CPU_MODEL
void compute_LUT_interpRow_mlir_8xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);
void compute_LUT_interpRow_mlir_4xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);
void compute_LUT_interpRow_mlir_2xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);
 
void LUT_interpRow_mlir(char *table, int i, char* row_ptr,
                        int lut_num_elements, int vector_size)
{
  LUT *const tab = (LUT *const)table;
    LUT_data_t* row = (LUT_data_t*)row_ptr;
    
    if (lut_num_elements > 0)
    {
    #ifndef IMP_FAST
      if (vector_size == 8)
      {
        compute_LUT_interpRow_mlir_8xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 4)
      {
        compute_LUT_interpRow_mlir_4xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 2)
      {
        compute_LUT_interpRow_mlir_2xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
       
    #else
      if (vector_size == 8)
      {
        compute_LUT_interpRow_mlir_8xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 4)
      {
        compute_LUT_interpRow_mlir_4xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 2)
      {
        compute_LUT_interpRow_mlir_2xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
    #endif      
    }
    
}
 
 
void LUT_problem_mlir( char *table, GlobalData_t val, int locind)
{
  LUT *const tab = (LUT *const)table;
  LUT_problem(tab, val, locind, tab->name );
}
#endif
 
#ifdef __cplusplus
}
#endif  // ifdef __cplusplus
 
}  // namespace limpet