Move material stiffness mat_mul outside innner loop

[dseyler]

Problem

Currently, the mat_mul operation to compute material stiffness for both struct and ustruct is contained inside an inner loop over each element's nodes. The mat_mul operation can be safely hoisted ouside the inner loop, resulting in a ~10% speedup timestep for a 500k element simulation.

This occurs around line 750 in sv_struct.cpp:

  for (int b = 0; b < eNoN; b++) {

    for (int a = 0; a < eNoN; a++) {

      // Geometric stiffness
      NxSNx = Nx(0,a)*S(0,0)*Nx(0,b) + Nx(1,a)*S(1,0)*Nx(0,b) +
              Nx(2,a)*S(2,0)*Nx(0,b) + Nx(0,a)*S(0,1)*Nx(1,b) +
              Nx(1,a)*S(1,1)*Nx(1,b) + Nx(2,a)*S(2,1)*Nx(1,b) +
              Nx(0,a)*S(0,2)*Nx(2,b) + Nx(1,a)*S(1,2)*Nx(2,b) +
              Nx(2,a)*S(2,2)*Nx(2,b);

      T1 = amd*N(a)*N(b) + afu*NxSNx;

      // Material Stiffness (Bt*D*B)
      mat_mul(Dm, Bm.rslice(b), DBm);

A similar pattern occurs around line 1410 in ustruct.cpp:

for (int b = 0; b < eNoNw; b++) {
    for (int a = 0; a < eNoNw; a++) {
      NxSNx = Nwx(0,a)*Siso(0,0)*Nwx(0,b)
       + Nwx(0,a)*Siso(0,1)*Nwx(1,b) + Nwx(0,a)*Siso(0,2)*Nwx(2,b)
       + Nwx(1,a)*Siso(1,0)*Nwx(0,b) + Nwx(1,a)*Siso(1,1)*Nwx(1,b)
       + Nwx(1,a)*Siso(1,2)*Nwx(2,b) + Nwx(2,a)*Siso(2,0)*Nwx(0,b)
       + Nwx(2,a)*Siso(2,1)*Nwx(1,b) + Nwx(2,a)*Siso(2,2)*Nwx(2,b);

      auto DBm = mat_mul(Dm, Bm.rslice(b));

Solution

The two mat_mul operations should occur once in the outer loop rather than being repeated for each node. I implemented this change in this branch and confirmed byte-identical results on the relevant test cases.

I profiled a 500k-element struct simulation before and after the change and found a 10% speedup per time step. Similar improvements should be seen in ustruct simulations.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[ktbolt]

@dseyler Very nice ! Did you find this using profiling tools or was this identified by an AI bot ?

The original svFSI code used the Fortran MATMUL intrinsic function for this matrix multiplication so the compiler probably optimized this out.

Note that

mat_mul(Dm, Bm.rslice(b), DBm);

is probably more efficient because

auto DBm = mat_mul(Dm, Bm.rslice(b));

allocates a new DBm; not sure if the compiler can optimize this.

[dseyler]

@ktbolt I used the macOS sampling profiler, Xcode, and a branch of svMultiPhysics with timing across different sections of the code (facilitated by AI). More issues to follow, saving us an additional 10- 20%!

On that note, this made me realize it might be a good idea to routinely profile our code to look for these kinds of optimizations. Curious what your thoughts are, but I can think of a few different ways to do this (in order of increasing complexity):

• A profiling script that developers can run and analyze before and after implementing a feature.
• A skill developers can use to instruct an AI agent to analyze changes made to a branch, identify what physics is affected, select a relevant test case (likely needs to be a larger mesh than our existing test cases), and profile and analyze the runtime.
• An automatic agent in GitHub that does the above checks after every pull request. I haven't looked into this much, but I think there's a way to do this.

[dseyler]

@ktbolt Is this ready for a PR after I implement the change you recommended?

[ktbolt]

@dseyler It would be good to create a profiling document describing what you've done and provide any scripts that you might have used; maybe a group meeting to discuss this and other optimization tips.

We do need to start looking for changes in the histor.dat files. I'm not sure if this should be done for the CI tests but definitely for the large simulations from the soon to be creaeted svMultiPhysics/Examples that would be run every couple of weeks or so.

[dseyler]

@ktbolt Sounds good! I can write something up and put together some scripts.

I noticed that there's a mat_mul6x3 function in mat_fun that matches the shape of the problem at hand. Is that supposed to be more efficient than generic mat_mul? If so, I can swap that in and check if it does the same math.

[ktbolt]

@dseyler mat_mul6x3 looks dangerous; maybe better not to use it.

[dseyler]

@ktbolt Confirmed that mat_mul6x3 does give slightly different results (I think the compiler fuses things differently), but if I write each component as a sequential sum, it gives byte identical results to the generic mat_mul function and runs 3.5 x faster per call, shaving another 3% off the total runtime.

An example component from the old mat_mul6x3 that doesn't match mat_mul:

c[0] = a[0]*b[0] + a[6]*b[1] + a[12]*b[2] + a[18]*b[3] + a[24]*b[4] + a[30]*b[5];

An example component from the new mat_mul6x3 that matches mat_mul exactly and passes all test cases:

s = a[0]*b[0];
  s += a[6]*b[1];
  s += a[12]*b[2];
  s += a[18]*b[3];
  s += a[24]*b[4];
  s += a[30]*b[5];
  c[0] = s;

Do you think it's worth swapping this in since it's called so many times throughout struct assembly? If so, would it be better to have it be a separate function from mat_mul, or to check the shapes of inputs within the generic mat_mul function and handle these cases differently?

I'm leaning toward the latter since it would be less confusing to read sv_struct.cpp, where the function would be exposed, and there may be other matrix multiplications in the codebase that could be similarly optimized and dispatched within mat_mul (they would take the generic path if the size doesn't match).

Update: I looked into this more and it looks like fixed-sized eigen maps seem to be the best way to handle these instances where the dimensions of mat_mul are known at compile time. Still gives about a 3% speedup. This increases to 10% if we repeat the same for other fixed-size operations like 3x3 * 3x3

Thinking something like this within mat_mul:

{ 
  // Fixed-shape fast path for 6x6 * 6x3. Runs 6x faster than the generic path
  if (A.nrows() == 6 && A.ncols() == 6 && B.nrows() == 6 && B.ncols() == 3 &&
      result.nrows() == 6 && result.ncols() == 3) {
    mat_mul6x6_6x3(A, B, result);
    return;
  }

where

static void mat_mul6x6_6x3(const Array<double>& A, const Array<double>& B, Array<double>& C)
{
  Eigen::Map<const Eigen::Matrix<double,6,6>> a(A.data());
  Eigen::Map<const Eigen::Matrix<double,6,3>> b(B.data());
  Eigen::Map<Eigen::Matrix<double,6,3>>       c(C.data());

  c.noalias() = a * b;
}```

[ktbolt]

@dseyler I would use Eigen since it is a cleaner implementation and we are using it in other places in the code.

[dseyler]

Sounds good! Let me know if you'd like me to open these as two separate PRs (one for moving mat_mul outside the loop and another for eigen) or if one PR is fine for both

[ktbolt]

One PR will be OK; not too many changes and they are recorded here.