Skip to content

Error in coordinate transformation in neBEM (undefined behavior)

Joseph Mc Kenna requested to merge jmckenna/garfieldpp:neBEM_errors into master

The compiler has found some undefined behavior in neBEM.c

[ 66%] Building CXX object CMakeFiles/Garfield.dir/NeBem/neBEMInterface.c.o
/builds/jmckenna/garfieldpp/NeBem/neBEM.c: In function 'double neBEM::ValueKnCh(int)':
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2095:54: warning: iteration 3 invokes undefined behavior [-Waggressive-loop-optimizations]
 2095 |           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
      |                             ~~~~~~~~~~~~~~~~~~~~~~~~~^
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2094:27: note: within this loop
 2094 |         for (int j = 0; j < 4; ++j) {
      |                         ~~^~~
[ 66%] Building CXX object CMakeFiles/Garfield.dir/NeBem/nrutil.c.o
[ 67%] Building CXX object CMakeFiles/Garfield.dir/NeBem/svdcmp.c.o
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2095:26: warning: iteration 3 invokes undefined behavior [-Waggressive-loop-optimizations]
 2095 |           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
      |           ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2093:25: note: within this loop
 2093 |       for (int i = 0; i < 4; ++i) {
      |                       ~~^~~
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2095:73: warning: array subscript 3 is above array bounds of 'double [3]' [-Warray-bounds]
 2095 |           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
      |                                                          ~~~~~~~~~~~~~~~^
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2091:14: note: while referencing 'InitialVector'
 2091 |       double InitialVector[3] = {xfld - pOrigin->X, yfld - pOrigin->Y, zfld - pOrigin->Z};
      |              ^~~~~~~~~~~~~
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2095:54: warning: array subscript 3 is above array bounds of 'double [3]' [-Warray-bounds]
 2095 |           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
      |                             ~~~~~~~~~~~~~~~~~~~~~~~~~^
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2095:26: warning: array subscript 3 is above array bounds of 'double [3]' [-Warray-bounds]
 2095 |           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
      |           ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/builds/jmckenna/garfieldpp/NeBem/neBEM.c:2092:14: note: while referencing 'FinalVector'
 2092 |       double FinalVector[3] = {0., 0., 0.};
      |              ^~~~~~~~~~~

Its clear this is a valid warning, we are reading and writing from beyond the length of the array

      double FinalVector[3] = {0., 0., 0.};
      for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
           FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
         }
       }

In this pull request I fix this for loop... however looking closer at this file, this error might have propagated from elsewhere in the same file. It would be excellent for an expert to review neBEM.c lines 826 - 860.

        { // Rotate point3D from global to local system
        double InitialVector[4];
        double TransformationMatrix[4][4] = {{0.0, 0.0, 0.0, 0.0},
                                                {0.0, 0.0, 0.0, 0.0},
                                                {0.0, 0.0, 0.0, 0.0},
                                                {0.0, 0.0, 0.0, 1.0}};
        DirnCosn3D *DirCos = &(EleArr+elesrc-1)->G.DC;
        double FinalVector[4];

        InitialVector[0] = xfld - xsrc; InitialVector[1] = yfld - ysrc;
        InitialVector[2] = zfld - zsrc; InitialVector[3] = 1.0;

        TransformationMatrix[0][0] = DirCos->XUnit.X;
        TransformationMatrix[0][1] = DirCos->XUnit.Y;
        TransformationMatrix[0][2] = DirCos->XUnit.Z;
        TransformationMatrix[1][0] = DirCos->YUnit.X;
        TransformationMatrix[1][1] = DirCos->YUnit.Y;
        TransformationMatrix[1][2] = DirCos->YUnit.Z;
        TransformationMatrix[2][0] = DirCos->ZUnit.X;
        TransformationMatrix[2][1] = DirCos->ZUnit.Y;
        TransformationMatrix[2][2] = DirCos->ZUnit.Z;

        for(int i = 0; i < 4; ++i)
        {
        FinalVector[i] = 0.0;
        for(int j = 0 ; j < 4; ++j)
        {
        FinalVector[i] += TransformationMatrix[i][j] * InitialVector[j];
        }
        }

        localP.X = FinalVector[0];
        localP.Y = FinalVector[1];
        localP.Z = FinalVector[2];
        } // Point3D rotated

There is a transformation in a 4-vector that then gets copied into a 3-vector. Do we want to resize FinalVector to [3] here, then reduce the for loops size to 3? Save CPU time by skipping unneeded float operations?

I also see this code basically repeated in line 1022 - 1058. Perhaps there should be a function for this transformation to avoid repetition?

Merge request reports