Do not unroll loops for jacobian computation

[IgorBaratta]

• Compilers are able to unroll small loops
• The loop form is more suited for sum-factorization

Generates:

  double J_c4 = 0.0;
  double J_c8 = 0.0;
  double J_c5 = 0.0;
  double J_c7 = 0.0;
  double J_c0 = 0.0;
  double J_c1 = 0.0;
  double J_c6 = 0.0;
  double J_c3 = 0.0;
  double J_c2 = 0.0;
  for (int ic = 0; ic < 4; ++ic)
  {
    J_c4 += coordinate_dofs[ic * 3 + 1] * FE9_C1_D010_Q421[0][0][0][ic];
    J_c8 += coordinate_dofs[ic * 3 + 2] * FE9_C2_D001_Q421[0][0][0][ic];
    J_c5 += coordinate_dofs[ic * 3 + 1] * FE9_C2_D001_Q421[0][0][0][ic];
    J_c7 += coordinate_dofs[ic * 3 + 2] * FE9_C1_D010_Q421[0][0][0][ic];
    J_c0 += coordinate_dofs[ic * 3] * FE8_C0_D100_Q421[0][0][0][ic];
    J_c1 += coordinate_dofs[ic * 3] * FE9_C1_D010_Q421[0][0][0][ic];
    J_c6 += coordinate_dofs[ic * 3 + 2] * FE8_C0_D100_Q421[0][0][0][ic];
    J_c3 += coordinate_dofs[ic * 3 + 1] * FE8_C0_D100_Q421[0][0][0][ic];
    J_c2 += coordinate_dofs[ic * 3] * FE9_C2_D001_Q421[0][0][0][ic];
  }

Instead of:

const double J_c4 = coordinate_dofs[1] * FE9_C1_D010_Q421[0][0][0][0] + coordinate_dofs[4] * FE9_C1_D010_Q421[0][0][0][1] + coordinate_dofs[7] * FE9_C1_D010_Q421[0][0][0][2] + coordinate_dofs[10] * FE9_C1_D010_Q421[0][0][0][3];
const double J_c8 = coordinate_dofs[2] * FE9_C2_D001_Q421[0][0][0][0] + coordinate_dofs[5] * FE9_C2_D001_Q421[0][0][0][1] + coordinate_dofs[8] * FE9_C2_D001_Q421[0][0][0][2] + coordinate_dofs[11] * FE9_C2_D001_Q421[0][0][0][3];
const double J_c5 = coordinate_dofs[1] * FE9_C2_D001_Q421[0][0][0][0] + coordinate_dofs[4] * FE9_C2_D001_Q421[0][0][0][1] + coordinate_dofs[7] * FE9_C2_D001_Q421[0][0][0][2] + coordinate_dofs[10] * FE9_C2_D001_Q421[0][0][0][3];
const double J_c7 = coordinate_dofs[2] * FE9_C1_D010_Q421[0][0][0][0] + coordinate_dofs[5] * FE9_C1_D010_Q421[0][0][0][1] + coordinate_dofs[8] * FE9_C1_D010_Q421[0][0][0][2] + coordinate_dofs[11] * FE9_C1_D010_Q421[0][0][0][3];
const double J_c0 = coordinate_dofs[0] * FE8_C0_D100_Q421[0][0][0][0] + coordinate_dofs[3] * FE8_C0_D100_Q421[0][0][0][1] + coordinate_dofs[6] * FE8_C0_D100_Q421[0][0][0][2] + coordinate_dofs[9] * FE8_C0_D100_Q421[0][0][0][3];
const double J_c1 = coordinate_dofs[0] * FE9_C1_D010_Q421[0][0][0][0] + coordinate_dofs[3] * FE9_C1_D010_Q421[0][0][0][1] + coordinate_dofs[6] * FE9_C1_D010_Q421[0][0][0][2] + coordinate_dofs[9] * FE9_C1_D010_Q421[0][0][0][3];
const double J_c6 = coordinate_dofs[2] * FE8_C0_D100_Q421[0][0][0][0] + coordinate_dofs[5] * FE8_C0_D100_Q421[0][0][0][1] + coordinate_dofs[8] * FE8_C0_D100_Q421[0][0][0][2] + coordinate_dofs[11] * FE8_C0_D100_Q421[0][0][0][3];
const double J_c3 = coordinate_dofs[1] * FE8_C0_D100_Q421[0][0][0][0] + coordinate_dofs[4] * FE8_C0_D100_Q421[0][0][0][1] + coordinate_dofs[7] * FE8_C0_D100_Q421[0][0][0][2] + coordinate_dofs[10] * FE8_C0_D100_Q421[0][0][0][3];
const double J_c2 = coordinate_dofs[0] * FE9_C2_D001_Q421[0][0][0][0] + coordinate_dofs[3] * FE9_C2_D001_Q421[0][0][0][1] + coordinate_dofs[6] * FE9_C2_D001_Q421[0][0][0][2] + coordinate_dofs[9] * FE9_C2_D001_Q421[0][0][0][3];

[chrisrichardson]

I guess this is OK, but I wonder if it makes much difference.
Which kernels will benefit from it? Surely all the time is spent in the main quadrature loop.

[IgorBaratta]

I guess this is OK, but I wonder if it makes much difference. Which kernels will benefit from it? Surely all the time is spent in the main quadrature loop.

Not much difference for simplices, but for sum factorization we need this loop structure.