Move unnecessary mat_mul outside loops and implemented fixed-sized mat_mul

[dseyler]

Current situation

After profiling a 500k element structural mechanics simulation, I noticed several mat_mul function calls that were unnecessarily repeated within loops, as discussed in issue 602.

These occurred at:

• mat_models.cpp lines 1695-1696
• sv_struct.cpp line 750
• ustruct.cpp line 144

and have now been moved outside their loops.

Additionally, while investigating, it was noticed that many of these mat_mul operations were of fixed-size arrays with sizes known at compilation. These can be optimized with fixed-size Eigen maps, which run 3-4x faster than the generic mat_mul operation.

New fixed-size functions were defined in mat_mul.cpp:

• mat_mul_fixed()
• mat_mul_fixed_rows()

These functions are then dispatched within mat_mul() so developers do not need to distinguish between these differences.

///
/// Used when sizes are known at compile time.
template <int M, int K, int N>
inline void mat_mul_fixed(const Array<double>& A, const Array<double>& B, Array<double>& C)
{
  Eigen::Map<const Eigen::Matrix<double, M, K>> a(A.data());
  Eigen::Map<const Eigen::Matrix<double, K, N>> b(B.data());
  Eigen::Map<Eigen::Matrix<double, M, N>>       c(C.data());

  c.noalias() = a * b;
}

And in mat_mul:

{ 
  // Fixed-shape fast paths for the products that dominate the element loops.
  //
  //   3x3 * 3x3     F*S in struct_3d; vx*Fi, ddev*Fit, Fi*ddev_Fit and the
  //                 potential-viscosity products in mat_models; F^T*F in cep
  //   3x3 * 3xeNoN  ddev*Nx_Fi and vx_Fi*Nx_Fi in the viscous tangent
  //   6x6 * 6x3     the material stiffness product D*B in struct_3d/ustruct_3d_m
  if (A.nrows() == 3 && A.ncols() == 3 && B.nrows() == 3 &&
      result.nrows() == 3 && result.ncols() == B.ncols()) {
    if (B.ncols() == 3) {
      mat_mul_fixed<3, 3, 3>(A, B, result);
    } else {
      mat_mul_fixed_rows<3, 3>(A, B, result);
    }
    return;
  }

In total, these changes decreased runtime by ~21%.

I suspect this number could increase much more if we correctly implement eigen arrays throughout sv_struct.cpp and mat_models.cpp as many arrays are of size nsd or enon, which have a discrete number of values and can be optimized similarly. I tested this on the viscosity models alone and found another 10% improvement (can open an issue about this, but the code is too messy for this PR)

Release Notes

• Moved mat_mul outside loops in several locations
• Implemented fixed-shape mat_mul functions
• mat_mul directs arrays with know sizes at compilation down these paths

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[claude]

Claude Code Review

This pull request is from a fork — automated review is disabled. A repository maintainer can comment @claude review to run a one-time review.

[Copilot]

Pull request overview

This PR targets solver hot paths by reducing redundant matrix multiplications inside element loops and accelerating common fixed-shape mat_mul cases via Eigen fixed-size/dynamic-column maps.

Changes:

• Hoists invariant mat_mul(Dm, Bm.rslice(b), DBm) out of inner node loops in sv_struct.cpp and ustruct.cpp.
• Hoists two mat_mul calls out of the eNoN loop in compute_visc_stress_newtonian() (mat_models.cpp).
• Introduces fixed-shape fast paths inside mat_fun::mat_mul(A,B,result) and routes the 2-arg mat_mul(A,B) overload through it; removes legacy mat_mul6x3 API.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 1 comment.

Show a summary per file

File	Description
Code/Source/solver/ustruct.cpp	Computes DBm = Dm*B once per b (outer loop) instead of per (a,b) pair.
Code/Source/solver/sv_struct.cpp	Computes DBm = Dm*B once per b (outer loop) to avoid redundant work in the inner loop.
Code/Source/solver/mat_models.cpp	Moves ddev*Nx_Fi and vx_Fi*Nx_Fi products outside the a loop after Nx_Fi is fully assembled.
Code/Source/solver/mat_fun.h	Removes the mat_mul6x3 declaration.
Code/Source/solver/mat_fun.cpp	Adds Eigen-based fixed-shape fast paths and delegates the 2-arg mat_mul overload to the 3-arg implementation; removes mat_mul6x3.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 72.78%. Comparing base (a3413e0) to head (f139eeb).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #608      +/-   ##
==========================================
+ Coverage   72.73%   72.78%   +0.05%     
==========================================
  Files         255      255              
  Lines       39290    39294       +4     
  Branches     6726     6727       +1     
==========================================
+ Hits        28576    28601      +25     
+ Misses      10472    10450      -22     
- Partials      242      243       +1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ktbolt]

@dseyler Templates should be defined in header files and they don't need to be declared inline.

Templates will use function overloading so you can rename mat_mul_fixed and mat_mul_fixed_rows to mat_mul.

[michelebucelli]

Thanks @dseyler! I left a few comments.

[michelebucelli]

I would suggest taking the opportunity to update the documentation of the functions being modified here, including details such as the meaning of inputs and outputs (which is fairly obvious here, but still), and whether the functions perform checks on their inputs.

In this case, this might look something like this:

Suggested change

	/// @brief Multiply two matrices.
	///
	/// Given the input matrices @f$A@f$ and @f$B@f$, computes and returns
	/// their product @f$AB@f$.
	///
	/// This function checks that the matrix dimensions are compatible, and
	/// throws an std::runtime_error exception otherwise.
	///
	/// Reproduces Fortran MATMUL.

1. I don't know if it's appropriate to keep the Fortran routine reference (it is probably only meaningful to those who know the original Fortran codebase, but it might still be useful during the transition).
2. If the runtime error is replaced as I suggest below, the comment would have to be updated.

Additionally, I suggest moving this documentation comment (and maybe others for other modified functions in this file) to the header file, rather than the source file. As far as I know, that is the standard place to put documentation, the logic being that someone who wants to use the function will first look it up in the header file (so that they can see the signature), and only look at the cpp file if they really need to know the implementation details (which they shouldn't care about in most normal situations).

[ktbolt]

Yes, let's remove the remembrance of things past like the Fortran references.

Documentation should be in both the header (API documentation) and .cpp (implementation details) files.

[michelebucelli]

I suggest using one of the exceptions from Core/Exception.h or FE/Common/FEException.h. I think that InvalidArgumentException is probably appropriate here.

I suggest to look up the function throw_if and use that to bundle the check and the exception raising in a single statement. The syntax of that function is somewhat obscure, but if you look up some places where it is called it should become clear.

[michelebucelli]

I suggest giving a documentation to this function, along the lines of what I suggested above. Since this function has both input and output arguments, I suggest using the Doxygen syntax @param[in] and @param[out] to document them individually, making this explicit (as opposed to, implied by their types).

[michelebucelli]

In the same vein as above, I suggest to expand the documentation here. Here it is especially important do document what the template parameters map to (I assume they're the dimensions of the matrices, but which of the dimensions each parameter corresponds to is left implicit). I believe that Doxygen has the syntax @tparam for template parameters.

(Since this function is only defined in the cpp file, I think it's appropriate that the documentation remain here).

[michelebucelli]

Same as above about documenting the meaning of template parameters.

[michelebucelli]

I think it would be cleaner to offer only one function between this (I'll call it option A below) and the overload that returns its result (which I'll call option B). I assume that both are currently used in the code, but a search-and-replace should quickly fix that if we resolve to only leave one.

The pros and cons of the two options, as far as I can tell, are as follows:

1. option A is faster than or equal to than option B, but this can probably be at least alleviated by designing option B and the Array class so that return value optimization takes place and move construction/move copy are possible;
2. option A has more cumbersome (old style idiom, less intuitive) syntax, and requires the caller to instantiate the target object (potentially two lines instead of one); option B supports the more modern syntax C = mat_mul(A, B), which has much better readability;
3. option B supports assigning the result to const outputs, which in some cases might allow the compiler to better optimize things out of loops when possible. In other words, with option B you will be able to do const auto C = mat_mul(A, B), while option A requires Array<double> C(...); mat_mul(A, B, C);, so that C cannot be made const.

Personally, I'd lean in favor of going with option A but implementing move constructor and move assignment for Array, so that copies can be made at basically zero cost if the object being copied-from is going to be discarded.

[michelebucelli]

@ktbolt

Templates should be defined in header files and they don't need to be declared inline

I guess the bigger question is whether we want to offer the option to the other parts of the library to call the fixed-size functions directly (in which case those functions should definitely be moved to the header) or not (in which case I think it's better if they stay in the anonymous namespace inside mat_mul.cpp). I slightly lean towards the latter, so that the other modules do not need to be aware of the low-level optimizations that mat_mul does under the hood.

Templates will use function overloading so you can rename mat_mul_fixed and mat_mul_fixed_rows to mat_mul.

This is true I think, but I would actually favor the explicit names, because I usually have a hard time remembering the overload resolution rules when the same function is both a template and not a template 😅 I suppose that, without template arguments, a call like mat_mul(A, B, C) would always resolve to the function that is not a template. Is this correct?

[ktbolt]

@dseyler There are probably too many mat_mul functions. Note that many uses of mat_mul can be replaced by* for example in ustruct.cpp

auto Pdev = mat_fun::mat_mul(F, Siso);

can be replaced by

auto Pdev = F * Siso;

[michelebucelli]

@ktbolt Good point, I hadn't noticed Array::operator*, that's even more C++-ish! Wouldn't it make sense then to only keep operator* and move all the special-case, fixed-size optimizations that @dseyler did into that?

[dseyler]

@ktbolt Good point, I hadn't noticed Array::operator*, that's even more C++-ish! Wouldn't it make sense then to only keep operator* and move all the special-case, fixed-size optimizations that @dseyler did into that?

Thanks for the feedback! Is there any reason that mat_mul couldn't be completely eliminated in this way?

Maybe for later, but I think a more detailed audit and optimization of our Array class would be a good idea to do sometime soon. I wouldn't be surprised if we could cut runtime in half or more by doing so.

[dseyler]

@ktbolt

Templates should be defined in header files and they don't need to be declared inline

I guess the bigger question is whether we want to offer the option to the other parts of the library to call the fixed-size functions directly (in which case those functions should definitely be moved to the header) or not (in which case I think it's better if they stay in the anonymous namespace inside mat_mul.cpp). I slightly lean towards the latter, so that the other modules do not need to be aware of the low-level optimizations that mat_mul does under the hood.

Templates will use function overloading so you can rename mat_mul_fixed and mat_mul_fixed_rows to mat_mul.

This is true I think, but I would actually favor the explicit names, because I usually have a hard time remembering the overload resolution rules when the same function is both a template and not a template 😅 I suppose that, without template arguments, a call like mat_mul(A, B, C) would always resolve to the function that is not a template. Is this correct?

@ktbolt If I'm understanding correctly, if the templates were defined in the header, then they would also have to be explicitly included in the other modules?

I agree that developers should ideally not have to think about whether their array sizes are known at compile time, and mat_mul (and Array.cpp/h) should all handle that behind the scenes.

I'm not too opinionated but would prefer whichever option allows the modules to be unaware of the nuances in these functions, but let someone who cares about it to quickly scan the lower-level code to see what optimizations are currently in place (I can imagine other optimal paths being implemented for other cases in the future).

Unrelated, but there are a few other spots in the code where the use of these fixed-size operations could be expanded if we template other functions that take nsd or enon as arguments. That would be a bit more exposed to developers I think.