[Bug tree-optimization/126551] New: missed vectorization due to bad reduction transformation in miniBUDE

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=126551

            Bug ID: 126551
           Summary: missed vectorization due to bad reduction
                    transformation in miniBUDE
           Product: gcc
           Version: 17.0
            Status: UNCONFIRMED
          Keywords: missed-optimization
          Severity: normal
          Priority: P3
         Component: tree-optimization
          Assignee: unassigned at gcc dot gnu.org
          Reporter: arsen at gcc dot gnu.org
                CC: pinskia at gcc dot gnu.org
  Target Milestone: ---
            Target: x86_64-linux-gnu

the following testcase:

  #include <stddef.h>
  #include <stdint.h>
  #include <cmath>
  #include <limits>

  #define ZERO 0.0f
  #define QUARTER 0.25f
  #define HALF 0.5f
  #define ONE 1.0f
  #define TWO 2.0f
  #define FOUR 4.0f
  #define CNSTNT 45.0f

  // Energy evaluation parameters
  #define HBTYPE_F 70
  #define HBTYPE_E 69
  #define HARDNESS 38.0f
  #define NPNPDIST 5.5f
  #define NPPDIST 1.0f

  static constexpr auto FloatMax = std::numeric_limits<float>::max();
  constexpr int PPWI = 1;

  struct __attribute__((__packed__)) Atom {
    float x, y, z;
    int32_t type;
  };

  struct __attribute__((__packed__)) FFParams {
    int32_t hbtype;
    float radius;
    float hphb;
    float elsc;
  };


  void fasten_main(size_t group, size_t ntypes, size_t nposes, size_t natlig,
size_t natpro,        //
                                   const Atom *protein, const Atom *ligand,    
                                    //
                                   const float *transforms_0, const float
*transforms_1, const float *transforms_2, //
                                   const float *transforms_3, const float
*transforms_4, const float *transforms_5, //
                                   const FFParams *forcefield, float *energies 
                                    //
    ) {

      float transform[3][4];
      float etot;

        int ix = group;

        // Compute transformation matrix
        const float sx = std::sin(transforms_0[ix]);
        const float cx = std::cos(transforms_0[ix]);
        const float sy = std::sin(transforms_1[ix]);
        const float cy = std::cos(transforms_1[ix]);
        const float sz = std::sin(transforms_2[ix]);
        const float cz = std::cos(transforms_2[ix]);

        transform[0][0] = cy * cz;
        transform[0][1] = sx * sy * cz - cx * sz;
        transform[0][2] = cx * sy * cz + sx * sz;
        transform[0][3] = transforms_3[ix];
        transform[1][0] = cy * sz;
        transform[1][1] = sx * sy * sz + cx * cz;
        transform[1][2] = cx * sy * sz - sx * cz;
        transform[1][3] = transforms_4[ix];
        transform[2][0] = -sy;
        transform[2][1] = sx * cy;
        transform[2][2] = cx * cy;
        transform[2][3] = transforms_5[ix];

        etot = 0.f;

      // Loop over ligand atoms
      for (int il = 0; il < natlig; il++) {
        // Load ligand atom data
        const Atom l_atom = ligand[il];
        const FFParams l_params = forcefield[l_atom.type];
        const int lhphb_ltz = l_params.hphb < 0.f;
        const int lhphb_gtz = l_params.hphb > 0.f;

        // Transform ligand atom
        float lpos_x, lpos_y, lpos_z;

          lpos_x = transform[0][3] + l_atom.x * transform[0][0] + l_atom.y *
transform[0][1] +
                      l_atom.z * transform[0][2];
          lpos_y = transform[1][3] + l_atom.x * transform[1][0] + l_atom.y *
transform[1][1] +
                      l_atom.z * transform[1][2];
          lpos_z = transform[2][3] + l_atom.x * transform[2][0] + l_atom.y *
transform[2][1] +
                      l_atom.z * transform[2][2];

        // Loop over protein atoms
        for (int ip = 0; ip < natpro; ip++) {
          // Load protein atom data
          const Atom p_atom = protein[ip];
          const FFParams p_params = forcefield[p_atom.type];

          const float radij = p_params.radius + l_params.radius;
          const float r_radij = ONE / radij;

          const float elcdst = (p_params.hbtype == HBTYPE_F && l_params.hbtype
== HBTYPE_F) ? FOUR : TWO;
          const float elcdst1 = (p_params.hbtype == HBTYPE_F && l_params.hbtype
== HBTYPE_F) ? QUARTER : HALF;
          const int type_E = ((p_params.hbtype == HBTYPE_E || l_params.hbtype
== HBTYPE_E));

          const int phphb_ltz = p_params.hphb < 0.f;
          const int phphb_gtz = p_params.hphb > 0.f;
          const int phphb_nz = p_params.hphb != 0.f;
          const float p_hphb = p_params.hphb * (phphb_ltz && lhphb_gtz ? -ONE :
ONE);
          const float l_hphb = l_params.hphb * (phphb_gtz && lhphb_ltz ? -ONE :
ONE);
          const float distdslv = (phphb_ltz ? (lhphb_ltz ? NPNPDIST : NPPDIST)
: (lhphb_ltz ? NPPDIST : -FloatMax));
          const float r_distdslv = ONE / distdslv;

          const float chrg_init = l_params.elsc * p_params.elsc;
          const float dslv_init = p_hphb + l_hphb;

            // Calculate distance between atoms
            const float x = lpos_x - p_atom.x;
            const float y = lpos_y - p_atom.y;
            const float z = lpos_z - p_atom.z;
            const float distij = std::sqrt(x * x + y * y + z * z);

            // Calculate the sum of the sphere radii
            const float distbb = distij - radij;

            const int zone1 = (distbb < ZERO);

            // Calculate steric energy
            etot += (ONE - (distij * r_radij)) * (zone1 ? TWO * HARDNESS :
0.f);

            // Calculate formal and dipole charge interactions
            float chrg_e = chrg_init * ((zone1 ? ONE : (ONE - distbb *
elcdst1)) * (distbb < elcdst ? ONE : ZERO));
            float neg_chrg_e = -std::abs(chrg_e);
            chrg_e = type_E ? neg_chrg_e : chrg_e;
            etot += chrg_e * CNSTNT;

            // Calculate the two cases for Nonpolar-Polar repulsive
interactions
            float coeff = (ONE - (distbb * r_distdslv));
            float dslv_e = dslv_init * ((distbb < distdslv && phphb_nz) ? ONE :
0.f);
            dslv_e *= (zone1 ? ONE : coeff);
            etot += dslv_e;
        }
     }

      // Write result
        energies[group] = etot * HALF;
    }

(extracted from https://github.com/UoB-HPC/miniBUDE OMP w/ PPWI=1)

... fails be vectorized in the inner (natpro) loop on x86-64 w/ -Ofast
-ffast-math -march=x86-64-v4.

applying the following patch:

@@ -123,19 +123,17 @@
           const int zone1 = (distbb < ZERO);

           // Calculate steric energy
-          etot += (ONE - (distij * r_radij)) * (zone1 ? TWO * HARDNESS : 0.f);

           // Calculate formal and dipole charge interactions
           float chrg_e = chrg_init * ((zone1 ? ONE : (ONE - distbb * elcdst1))
* (distbb < elcdst ? ONE : ZERO));
           float neg_chrg_e = -std::abs(chrg_e);
           chrg_e = type_E ? neg_chrg_e : chrg_e;
-          etot += chrg_e * CNSTNT;

           // Calculate the two cases for Nonpolar-Polar repulsive interactions
           float coeff = (ONE - (distbb * r_distdslv));
           float dslv_e = dslv_init * ((distbb < distdslv && phphb_nz) ? ONE :
0.f);
           dslv_e *= (zone1 ? ONE : coeff);
-          etot += dslv_e;
+          etot += ((ONE - (distij * r_radij)) * (zone1 ? TWO * HARDNESS :
0.f)) + (chrg_e * CNSTNT) + dslv_e;
       }
    }

(i.e. merging all of the writes into etot into one write)

... makes it vectorize on x86.

Pinski briefly looked at the testcase, and said:

22:29:49 <pinskia> though it looks like the reduction of etot is not being
handled correctly
22:29:57 <pinskia> etot[0] that is
22:30:44 <pinskia> missed:   reduction used in loop
22:31:38 <pinskia> _359 = etot$_93 + _360;
22:31:38 <pinskia> ...
22:31:38 <pinskia> _355 = _300 ? _359 : etot$_93;
22:33:10 <pinskia>           etot[l] += (ONE - (distij * r_radij)) * (zone1 ?
TWO * HARDNESS : 0.f);
22:34:06 <pinskia> got incorrectly changed into the above rather than temp =
_300 ? _360 : 0.0f; _355 = temp + etot$93;
22:34:16 <pinskia> this is supposed to handled in ifcvt
22:35:48 <pinskia> is_cond_scalar_reduction is supposed to detect that and
create the reduction correctly but it looks like it is not

(note that the testcase declares 'etot' as a variable and not an array, in the
original code it was an array of 1 element)

... this lead me to try merging the +=s

the vectorized version seems to produce a ~15% result on the bm2 deck benchmark
(when compiled with 16.1, I didn't run the benchmark with trunk)