Skip to content

Move unnecessary mat_mul outside loops and implemented fixed-sized mat_mul #608

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dseyler wants to merge 8 commits into
base: main
Choose a base branch
from
Open

Move unnecessary mat_mul outside loops and implemented fixed-sized mat_mul #608

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e4fcd04
Moved D*B mat_mul out of inner loop in sv_struct.cpp. Results in ~10%…
dseyler Aug 18, 2026
4342ce0
Moved D*B mat_mul outside inner loop and cleaned up commenting
dseyler Aug 18, 2026
b709907
compute D*B in place in ustruct_3d_m
dseyler Aug 19, 2026
adc276b
Moved mat_mul outside loop in comput_visc_stress_newtownian. Reduced …
dseyler Aug 19, 2026
2eda477
Merge branch 'SimVascular:main' into perf/hoist-invariant-matmul-stru…
dseyler Aug 19, 2026
6e41b0b
Optimized mat_mul to use fixed-size eigen products when shapes are kn…
dseyler Aug 19, 2026
88fa9ca
Added delegation so more mat_mul calls reach fixed-shape path
dseyler Aug 19, 2026
f139eeb
Update Code/Source/solver/mat_fun.cpp
dseyler Aug 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 58 additions & 39 deletions Code/Source/solver/mat_fun.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -507,15 +507,10 @@ mat_mul(const Array<double>& A, const Array<double>& B)

Array<double> 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
// Fall back onto the overload of this function that takes the result matrix
// as an output argument.
mat_mul(A, B, result);

return result;
}
Expand All @@ -524,8 +519,62 @@ mat_mul(const Array<double>& A, const Array<double>& 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.
Comment on lines +524 to +526

@michelebucelli michelebucelli Aug 20, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

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;
}

/// @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.
Comment on lines +537 to +541

@michelebucelli michelebucelli Aug 20, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above about documenting the meaning of template parameters.

template <int M, int K>
inline void mat_mul_fixed_rows(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, Eigen::Dynamic>> b(B.data(), K, B.ncols());
Eigen::Map<Eigen::Matrix<double, M, Eigen::Dynamic>> c(C.data(), M, C.ncols());

c.noalias() = a * b;
}

} // namespace

void mat_mul(const Array<double>& A, const Array<double>& B, Array<double>& result)

@michelebucelli michelebucelli Aug 20, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 michelebucelli Aug 20, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

{
// 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) {
Comment on lines +556 to +564
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();
Expand Down Expand Up @@ -900,36 +949,6 @@ transpose(const Array<double>& A)
return result;
}

void mat_mul6x3(const Array<double>& A, const Array<double>& B, Array<double>& 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
}

};

Expand Down
1 change: 0 additions & 1 deletion Code/Source/solver/mat_fun.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,6 @@ namespace mat_fun {
Vector<double> mat_mul(const Array<double>& A, const Vector<double>& v);
Array<double> mat_mul(const Array<double>& A, const Array<double>& B);
void mat_mul(const Array<double>& A, const Array<double>& B, Array<double>& result);
void mat_mul6x3(const Array<double>& A, const Array<double>& B, Array<double>& C);

Array<double> mat_symm(const Array<double>& A, const int nd);
Array<double> mat_symm_prod(const Vector<double>& u, const Vector<double>& v, const int nd);
Expand Down
5 changes: 3 additions & 2 deletions Code/Source/solver/mat_models.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
6 changes: 3 additions & 3 deletions Code/Source/solver/sv_struct.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -735,6 +735,9 @@ 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)
mat_mul(Dm, Bm.rslice(b), DBm);

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

// Geometric stiffness
Expand All @@ -746,9 +749,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) +
Expand Down
7 changes: 5 additions & 2 deletions Code/Source/solver/ustruct.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1403,16 +1403,19 @@ 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<double> DBm(6,3);

for (int b = 0; b < eNoNw; 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)
+ 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));

// 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) +
Expand Down
Loading