From e4fcd04fe81e6f544617db86c24001723f2a6d58 Mon Sep 17 00:00:00 2001 From: dseyler Date: Tue, 18 Aug 2026 09:33:23 -0700 Subject: [PATCH 1/7] Moved D*B mat_mul out of inner loop in sv_struct.cpp. Results in ~10% speedup --- Code/Source/solver/sv_struct.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index 9ab5feca6..b2fa04ce0 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -735,6 +735,12 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, for (int b = 0; b < eNoN; b++) { + // Material stiffness (D*B). This depends only on b, so compute it once per + // b rather than once per (a,b) pair: Dm and Bm are both filled before these + // loops and are only read inside them, and mat_mul overwrites DBm rather + // than accumulating into it. Same arithmetic, eNoN times fewer calls. + mat_mul(Dm, Bm.rslice(b), DBm); + for (int a = 0; a < eNoN; a++) { // Geometric stiffness @@ -746,9 +752,6 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, T1 = amd*N(a)*N(b) + afu*NxSNx; - // Material Stiffness (Bt*D*B) - mat_mul(Dm, Bm.rslice(b), DBm); - // dM1/du1 // Material stiffness: Bt*D*B BmDBm = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + From 4342ce0af5417b4e56eb60fff819a0161ce5ccdf Mon Sep 17 00:00:00 2001 From: dseyler Date: Tue, 18 Aug 2026 10:06:31 -0700 Subject: [PATCH 2/7] Moved D*B mat_mul outside inner loop and cleaned up commenting --- Code/Source/solver/sv_struct.cpp | 5 +---- Code/Source/solver/ustruct.cpp | 5 +++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index b2fa04ce0..41848e485 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -735,10 +735,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, for (int b = 0; b < eNoN; b++) { - // Material stiffness (D*B). This depends only on b, so compute it once per - // b rather than once per (a,b) pair: Dm and Bm are both filled before these - // loops and are only read inside them, and mat_mul overwrites DBm rather - // than accumulating into it. Same arithmetic, eNoN times fewer calls. + // Material stiffness (D*B) mat_mul(Dm, Bm.rslice(b), DBm); for (int a = 0; a < eNoN; a++) { diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 69be045c7..52f8c9426 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -1404,6 +1404,9 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, double Tv{0.0}, Ku{0.0}; for (int b = 0; b < eNoNw; b++) { + + auto DBm = mat_mul(Dm, Bm.rslice(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) @@ -1411,8 +1414,6 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + 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)); - // dM1_dV1 + af/am *dM_1/dU_1 BtDB = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + Bm(2,0,a)*DBm(2,0) + Bm(3,0,a)*DBm(3,0) + From b70990702946c2af36269442e53c365ec97a417f Mon Sep 17 00:00:00 2001 From: dseyler Date: Wed, 19 Aug 2026 09:34:26 -0700 Subject: [PATCH 3/7] compute D*B in place in ustruct_3d_m --- Code/Source/solver/ustruct.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 52f8c9426..a698050e7 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -1403,9 +1403,11 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, double NxSNx{0.0}, BtDB{0.0}; double Tv{0.0}, Ku{0.0}; + Array DBm(6,3); + for (int b = 0; b < eNoNw; b++) { - auto DBm = mat_mul(Dm, Bm.rslice(b)); + mat_mul(Dm, Bm.rslice(b), DBm); for (int a = 0; a < eNoNw; a++) { NxSNx = Nwx(0,a)*Siso(0,0)*Nwx(0,b) From adc276b56e1db7d58607405be9f1e6b251bb93a9 Mon Sep 17 00:00:00 2001 From: dseyler Date: Wed, 19 Aug 2026 13:12:39 -0700 Subject: [PATCH 4/7] Moved mat_mul outside loop in comput_visc_stress_newtownian. Reduced runtime by 6% --- Code/Source/solver/mat_models.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index cda281a68..caef69cea 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -1692,10 +1692,11 @@ void compute_visc_stress_newtonian(const double mu, const int eNoN, const Array< Nx_Fi(i,a) += Nx(j,a) * Fi(j,i); } } - ddev_Nx_Fi = mat_mul(ddev, Nx_Fi); - vx_Fi_Nx_Fi = mat_mul(vx_Fi, Nx_Fi); } + mat_mul(ddev, Nx_Fi, ddev_Nx_Fi); + mat_mul(vx_Fi, Nx_Fi, vx_Fi_Nx_Fi); + // 2nd Piola-Kirchhoff stress due to viscosity // Svis = 2 * mu * J * F^-1 * d_dev * F^-T auto Fit = transpose(Fi); From 6e41b0b1e3d667b1b86b9ab169243c0cf9e7b6e3 Mon Sep 17 00:00:00 2001 From: dseyler Date: Wed, 19 Aug 2026 09:40:53 -0700 Subject: [PATCH 5/7] Optimized mat_mul to use fixed-size eigen products when shapes are known at compile time --- Code/Source/solver/mat_fun.cpp | 84 ++++++++++++++++++++++------------ Code/Source/solver/mat_fun.h | 1 - 2 files changed, 54 insertions(+), 31 deletions(-) diff --git a/Code/Source/solver/mat_fun.cpp b/Code/Source/solver/mat_fun.cpp index 582db47da..cbc241d18 100644 --- a/Code/Source/solver/mat_fun.cpp +++ b/Code/Source/solver/mat_fun.cpp @@ -524,8 +524,62 @@ mat_mul(const Array& A, const Array& B) /// /// Compute result directly into the passed argument. // +namespace { + +/// @brief Fixed-shape matrix product, C = A*B, mapped onto the existing buffers. +/// +/// Used when sizes are known at compile time. +template +inline void mat_mul_fixed(const Array& A, const Array& B, Array& C) +{ + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data()); + Eigen::Map> c(C.data()); + + c.noalias() = a * b; +} + +/// @brief As mat_mul_fixed, but with the column count known only at run time. +/// +/// Used where the right operand has one column per element node, so its width +/// depends on the element type. The row counts are still compile-time, which is +/// where most of the benefit comes from. +template +inline void mat_mul_fixed_rows(const Array& A, const Array& B, Array& C) +{ + Eigen::Map> a(A.data()); + Eigen::Map> b(B.data(), K, B.ncols()); + Eigen::Map> c(C.data(), M, C.ncols()); + + c.noalias() = a * b; +} + +} // namespace + void mat_mul(const Array& A, const Array& B, Array& result) { + // 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; + } + + if (A.nrows() == 6 && A.ncols() == 6 && B.nrows() == 6 && B.ncols() == 3 && + result.nrows() == 6 && result.ncols() == 3) { + mat_mul_fixed<6, 6, 3>(A, B, result); + return; + } + int A_num_rows = A.nrows(); int A_num_cols = A.ncols(); int B_num_rows = B.nrows(); @@ -900,36 +954,6 @@ transpose(const Array& A) return result; } -void mat_mul6x3(const Array& A, const Array& B, Array& C) -{ - #define mat_mul6x3_unroll - #ifdef mat_mul6x3_unroll - auto a = A.data(); - auto b = B.data(); - auto c = C.data(); - - 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]; - c[1] = a[1]*b[0] + a[7]*b[1] + a[13]*b[2] + a[19]*b[3] + a[25]*b[4] + a[31]*b[5]; - c[2] = a[2]*b[0] + a[8]*b[1] + a[14]*b[2] + a[20]*b[3] + a[26]*b[4] + a[32]*b[5]; - c[3] = a[3]*b[0] + a[9]*b[1] + a[15]*b[2] + a[21]*b[3] + a[27]*b[4] + a[33]*b[5]; - c[4] = a[4]*b[0] + a[10]*b[1] + a[16]*b[2] + a[22]*b[3] + a[28]*b[4] + a[34]*b[5]; - c[5] = a[5]*b[0] + a[11]*b[1] + a[17]*b[2] + a[23]*b[3] + a[29]*b[4] + a[35]*b[5]; - - c[6] = a[0]*b[6] + a[6]*b[7] + a[12]*b[8] + a[18]*b[9] + a[24]*b[10] + a[30]*b[11]; - c[7] = a[1]*b[6] + a[7]*b[7] + a[13]*b[8] + a[19]*b[9] + a[25]*b[10] + a[31]*b[11]; - c[8] = a[2]*b[6] + a[8]*b[7] + a[14]*b[8] + a[20]*b[9] + a[26]*b[10] + a[32]*b[11]; - c[9] = a[3]*b[6] + a[9]*b[7] + a[15]*b[8] + a[21]*b[9] + a[27]*b[10] + a[33]*b[11]; - c[10] = a[4]*b[6] + a[10]*b[7] + a[16]*b[8] + a[22]*b[9] + a[28]*b[10] + a[34]*b[11]; - c[11] = a[5]*b[6] + a[11]*b[7] + a[17]*b[8] + a[23]*b[9] + a[29]*b[10] + a[35]*b[11]; - - c[12] = a[0]*b[12] + a[6]*b[13] + a[12]*b[14] + a[18]*b[15] + a[24]*b[16] + a[30]*b[17]; - c[13] = a[1]*b[12] + a[7]*b[13] + a[13]*b[14] + a[19]*b[15] + a[25]*b[16] + a[31]*b[17]; - c[14] = a[2]*b[12] + a[8]*b[13] + a[14]*b[14] + a[20]*b[15] + a[26]*b[16] + a[32]*b[17]; - c[15] = a[3]*b[12] + a[9]*b[13] + a[15]*b[14] + a[21]*b[15] + a[27]*b[16] + a[33]*b[17]; - c[16] = a[4]*b[12] + a[10]*b[13] + a[16]*b[14] + a[22]*b[15] + a[28]*b[16] + a[34]*b[17]; - c[17] = a[5]*b[12] + a[11]*b[13] + a[17]*b[14] + a[23]*b[15] + a[29]*b[16] + a[35]*b[17]; - #endif -} }; diff --git a/Code/Source/solver/mat_fun.h b/Code/Source/solver/mat_fun.h index 5428479c7..caf83b628 100644 --- a/Code/Source/solver/mat_fun.h +++ b/Code/Source/solver/mat_fun.h @@ -97,7 +97,6 @@ namespace mat_fun { Vector mat_mul(const Array& A, const Vector& v); Array mat_mul(const Array& A, const Array& B); void mat_mul(const Array& A, const Array& B, Array& result); - void mat_mul6x3(const Array& A, const Array& B, Array& C); Array mat_symm(const Array& A, const int nd); Array mat_symm_prod(const Vector& u, const Vector& v, const int nd); From 88fa9ca509f4e58bf48e5869f55199a07cd15770 Mon Sep 17 00:00:00 2001 From: dseyler Date: Wed, 19 Aug 2026 11:51:14 -0700 Subject: [PATCH 6/7] Added delegation so more mat_mul calls reach fixed-shape path --- Code/Source/solver/mat_fun.cpp | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/Code/Source/solver/mat_fun.cpp b/Code/Source/solver/mat_fun.cpp index cbc241d18..d2e4c6592 100644 --- a/Code/Source/solver/mat_fun.cpp +++ b/Code/Source/solver/mat_fun.cpp @@ -507,15 +507,9 @@ mat_mul(const Array& A, const Array& B) Array result(A_num_rows, B_num_cols); - for (int i = 0; i < A_num_rows; i++) { - for (int j = 0; j < B_num_cols; j++) { - double sum = 0.0; - for (int k = 0; k < A_num_cols; k++) { - sum += A(i,k) * B(k,j); - } - result(i,j) = sum; - } - } + // Delegate rather than repeat the loop, so callers of this form reach the + // fixed-shape fast paths too. Most call sites use this overload. + mat_mul(A, B, result); return result; } From f139eeba08612ab65aee488aa38db7f4c01da2f5 Mon Sep 17 00:00:00 2001 From: dseyler <43078283+dseyler@users.noreply.github.com> Date: Thu, 20 Aug 2026 12:48:04 -0700 Subject: [PATCH 7/7] Update Code/Source/solver/mat_fun.cpp Co-authored-by: Michele Bucelli --- Code/Source/solver/mat_fun.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/mat_fun.cpp b/Code/Source/solver/mat_fun.cpp index d2e4c6592..bb0a4c081 100644 --- a/Code/Source/solver/mat_fun.cpp +++ b/Code/Source/solver/mat_fun.cpp @@ -508,7 +508,8 @@ mat_mul(const Array& A, const Array& B) Array result(A_num_rows, B_num_cols); // Delegate rather than repeat the loop, so callers of this form reach the - // fixed-shape fast paths too. Most call sites use this overload. + // Fall back onto the overload of this function that takes the result matrix + // as an output argument. mat_mul(A, B, result); return result;