LLNL / Silo

Mesh and Field I/O Library and Scientific Database
https://silo.llnl.gov
Other
30 stars 24 forks source link

Adjust extents calc logic for void aliasing #376

Closed markcmiller86 closed 2 months ago

markcmiller86 commented 6 months ago

A user running with icpx compiler magic wrapper encountered cases where UCD mesh extents were calculated incorrectly. I creaated a simple reproducer (below). In the code below, if the intermediate casts are not used, icpx with optimization at -O3 will get it wrong. The pointers will be bad. I was told adding -fno-strict-aliasing is the correct solution. But, I wanted a solution in the code itself.

#include <stdio.h>

#define DB_DOUBLE 20
#define DB_FLOAT 19

#define MAX(X,Y)        ((X)>(Y)?(X):(Y))
#define MIN(X,Y)        ((X)<(Y)?(X):(Y))

typedef void const *DBVCP2_t;

int
UM_CalcExtents(DBVCP2_t coord_arrays, int datatype, int ndims, int nnodes,
               void *min_extents, void *max_extents)
{
    int            i, j;
    double       **dcoord_arrays = NULL;
    double        *dmin_extents = NULL, *dmax_extents = NULL;
    float         *fmin_extents = NULL, *fmax_extents = NULL;
    float        **fcoord_arrays = NULL;

    if (nnodes <= 0) return 0;

    if (datatype == DB_DOUBLE) {

        dmin_extents = (double *)min_extents;
        dmax_extents = (double *)max_extents;
        double const * const* _dcoord_arrays = (double const* const*) coord_arrays;
        double const * dcoord_arrays[3] = {_dcoord_arrays[0], _dcoord_arrays[1], _dcoord_arrays[2]};

        /* Initialize extent arrays */
        for (i = 0; i < ndims; i++) {
            dmin_extents[i] = dcoord_arrays[i][0];
            dmax_extents[i] = dcoord_arrays[i][0];
        }

        for (i = 0; i < ndims; i++) {
            for (j = 0; j < nnodes; j++) {
                dmin_extents[i] = MIN(dmin_extents[i], dcoord_arrays[i][j]);
                dmax_extents[i] = MAX(dmax_extents[i], dcoord_arrays[i][j]);
            }
        }

    }
    else {
        fmin_extents = (float *)min_extents;
        fmax_extents = (float *)max_extents;
        float const * const* _fcoord_arrays = (float const* const*) coord_arrays;
        float const * fcoord_arrays[3] = {_fcoord_arrays[0], _fcoord_arrays[1], _fcoord_arrays[2]};

        /* Initialize extent arrays */
        for (i = 0; i < ndims; i++) {
            fmin_extents[i] = fcoord_arrays[i][0];
            fmax_extents[i] = fcoord_arrays[i][0];
        }

        for (i = 0; i < ndims; i++) {
            for (j = 0; j < nnodes; j++) {
                fmin_extents[i] = MIN(fmin_extents[i], fcoord_arrays[i][j]);
                fmax_extents[i] = MAX(fmax_extents[i], fcoord_arrays[i][j]);
            }
        }

    }

    return 0;
}

int main(int argc, char **argv)
{
    double coords[3][10] = {
         {1, 1, 1, 1, 1, 2, 2, 2, 2, 2},
         {3, 4, 3, 4, 3, 4, 4, 5, 4, 5},
         {0, 1, 0, 2, 0, 3, 0, 4, 0, 5}};
    float fcoords[3][10] = {
         {1, 1, 1, 1, 1, 2, 2, 2, 2, 2},
         {3, 4, 3, 4, 3, 4, 4, 5, 4, 5},
         {0, 1, 0, 2, 0, 3, 0, 4, 0, 5}};
    /*void *pcoords[3] = {&coords[0][0], &coords[1][0], &coords[2][0]};*/
    void *pcoords[3];
    double mins[3] = {-1, -1, -1}, maxs[3] = {99,99,99};
    float fmins[3] = {-1, -1, -1}, fmaxs[3] = {99,99,99};

    pcoords[0] = &coords[0][0];
    pcoords[1] = &coords[1][0];
    pcoords[2] = &coords[2][0];
    UM_CalcExtents(pcoords, DB_DOUBLE, 3, 10, mins, maxs);
    printf("%s: mins = %16.16g, %16.16g, %16.16g\n", argv[0], mins[0], mins[1], mins[2]);
    printf("%s: maxs = %16.16g, %16.16g, %16.16g\n", argv[0], maxs[0], maxs[1], maxs[2]);
    pcoords[0] = &fcoords[0][0];
    pcoords[1] = &fcoords[1][0];
    pcoords[2] = &fcoords[2][0];
    UM_CalcExtents(pcoords, DB_FLOAT, 3, 10, fmins, fmaxs);
    printf("%s: mins = %16.16g, %16.16g, %16.16g\n", argv[0], (double) fmins[0], (double) fmins[1], (double) fmins[2]);
    printf("%s: maxs = %16.16g, %16.16g, %16.16g\n", argv[0], (double) fmaxs[0], (double) fmaxs[1], (double) fmaxs[2]);

    return 0;
}