Optimize double dot product and dyadic product

[dseyler]

Problem

Many of the operations in mat_fun.h could be sped up with a better use of Eigen and by reordering nested loops such that the innermost loop steps along the index in which memory is stored consecutively.

Specifically, double_dot_product(), dyadic_product, and symmetric_dyadic_product are good candidates for these improvements as they are called several times for each gauss point in struct assembly.

Solution

I propose the following changes, resulting in a ~12% decrease in struct/ustruct assembly runtime, which I've implemented in this branch.

double_dot_product()

Add a fast path for double dot products that are contracted along the last two dimensions, which occurs frequently when computing the isochoric elasticity tensor in bar_to_iso. The mechanism is as follows: nsd x nsd x nsd x nsd tensors are stored as a (nsd x nsd) x (nsd x nsd) matrix. As a result, if contracting along the last two dimensions, this can be expressed as a * b.transpose rather than using the general Eigen .contract(). Similarly, this could also be applied for contraction along the first two dimensions as a.transpose * b, but this pattern currently doesn't appear anywhere in the code, so I only implemented it for contraction along the last two dimensions.

// Fast path for dimsA = dimsB = {2,3}: C_ijmn = A_ijkl * B_mnkl.
        if (dimsA[0] == 2 && dimsA[1] == 3 && dimsB[0] == 2 && dimsB[1] == 3) {
            constexpr int N = nsd * nsd;
            Tensor<nsd> C;
            Eigen::Map<const Eigen::Matrix<double, N, N>> a(A.data());
            Eigen::Map<const Eigen::Matrix<double, N, N>> b(B.data());
            Eigen::Map<Eigen::Matrix<double, N, N>> c(C.data());
            c.noalias() = a * b.transpose();
            return C;
        }

dyadic_product()

Similarly, the nested for loops in the dyadic product can be written as an Eigen outer product:

dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
        // Initialize the result tensor
        Tensor<nsd> C;
        constexpr int N = nsd * nsd;

        Eigen::Map<const Eigen::Matrix<double, N, 1>> a(A.data());
        Eigen::Map<const Eigen::Matrix<double, N, 1>> b(B.data());
        Eigen::Map<Eigen::Matrix<double, N, N>> c(C.data());
        c.noalias() = a * b.transpose();

        return C;

symmetric_dyadic_product()

I couldn't figure out how to remove the nested for-loops entirely for this function, since the two matrices A and B each supply one of the first two indices and one of the last two, so the result doesn't reduce to a single outer product. Two of the four for loops can be eliminated, though, leaving only the last two indices, such that the indexed block is stored consecutively in memory. For context, most of the speedup in the other functions can be achieved just by correctly ordering the loops without switching to Eigen, so I think eliminating the other two loops wouldn't drastically improve performance.

template <int nsd>
    Tensor<nsd>
    symmetric_dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
        Tensor<nsd> C;

        // Compute the symmetric product: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk)
        for (int l = 0; l < nsd; ++l) {
            for (int k = 0; k < nsd; ++k) {
                Eigen::Map<Eigen::Matrix<double, nsd, nsd>> blk(C.data() + nsd * nsd * (k + nsd * l));
                blk.noalias() = 0.5 * (A.col(k) * B.col(l).transpose()
                                     + A.col(l) * B.col(k).transpose());
            }
        }

        // Return the symmetric product
        return C;

Additional context

No response

Code of Conduct

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

[ktbolt]

@dseyler What's the motivation for this Issue ? Are struct and ustruct simulations especially slow ? Or did you do a performance analysis and see that a significant amount of time was being spent in the above mentioned functions ?

How much does the total execution time decrease say for a 30 minute simulation ?

I am thinking that it might be better to spend time cleaning up things at a higher level, like bar_to_iso. If you could re-organize that code to chain operations in Eigen then the code would be cleaner and have better performance.

[dseyler]

@ktbolt Yes, struct and ustruct are particularly slow. As discussed here, our struct simulations could be running 4x faster with just a few fixes, mostly affecting assembly, which is 8x slower than it needs to be. The proposed changes in this issue amount to only a ~10% speedup, but don't require any changes to call sites.

The linked discussion also includes a cleanup of bar_to_iso, but the changes mentioned in this issue address the bulk of bar_to_iso runtime.

Do you think it would be better to rewrite bar_to_iso to eliminate references to double_dot_product and dyadic_product() entirely and just inline them with Eigen?

[ktbolt]

@dseyler I have not done much with Eigen so I don't know it as well as you but you should get better optimization by rewriting bar_to_iso to chain Eigen operations, for example, something like

  auto CC_iso_1 = PP.contract( transpose<nsd>(CC_bar.contract(PP, contractionDims)), contractionDims) +
                  (-2.0/nsd) * (dyadic_product<nsd>(Ci, S_iso) + dyadic_product<nsd>(S_iso, Ci));

The C++ compiler will probably inline dyadic_product but I don't know if an expression template will be created from this; you will need to check.

You do need to be careful using Eigen methods like contract in chained expressions.

The performance gains you already see are quite good so you will need to decide if you want to spend more time on this.