lock_error

void grid::compute_boundary_interactions_multipole_multipole(gsolve_type type, const std::vector<boundary_interaction_type>& ilist_n_bnd,
    const gravity_boundary_type& mpoles) {
    PROF_BEGIN;

    std::vector<space_vector> const& com0 = com[0];
//    for (integer si = 0; si != ilist_n_bnd.size(); ++si) {
    hpx::parallel::for_loop(
        hpx::parallel::par, 0, ilist_n_bnd.size(),
        [&mpoles, &com0, &ilist_n_bnd, type, this](std::size_t si) {

            taylor<4, simd_vector> m0;
            taylor<4, simd_vector> n0;
            std::array<simd_vector, NDIM> dX;
            std::array<simd_vector, NDIM> X;
            space_vector Y;

            boundary_interaction_type const& bnd = ilist_n_bnd[si];
            integer index = mpoles.is_local ? bnd.second : si;
            load_multipole(m0, Y, mpoles, index, false);
            const integer list_size = bnd.first.size();
            for (integer li = 0; li < list_size; li += simd_len) {
                for (integer i = 0; i != simd_len && li + i < list_size; ++i) {
                    const integer iii0 = bnd.first[li + i];
                    space_vector const& com0iii0 = com0[iii0];
                    for (integer d = 0; d < NDIM; ++d) {
                        X[d][i] = com0iii0[d];
                    }
                    if (type == RHO) {
                        multipole const& Miii0 = M[iii0];
                        real const tmp = m0()[i] / Miii0();
                        for (integer j = taylor_sizes[2]; j != taylor_sizes[3]; ++j) {
                            n0[j][i] = m0[j][i] - Miii0[j] * tmp;
                        }
                    }
                }
                if (type != RHO) {
                    for (integer j = taylor_sizes[2]; j != taylor_sizes[3]; ++j) {
                        n0[j] = ZERO;
                    }
                }
                for (integer d = 0; d < NDIM; ++d) {
                    dX[d] = X[d] - simd_vector(Y[d]);
                }

                taylor<5, simd_vector> D;
                taylor<4, simd_vector> A0;
                std::array<simd_vector, NDIM> B0 = {
                    simd_vector(0.0), simd_vector(0.0), simd_vector(0.0)
                };

                D.set_basis(dX);

    //             A0() = m0() * D();
    //             for (integer a = 0; a < NDIM; ++a) {
    //                 if (type != RHO) {
    //                     A0() -= m0(a) * D(a);
    //                 }
    //                 for (integer b = a; b < NDIM; ++b) {
    //                     A0() += m0(a, b) * D(a, b) * (real(1) / real(2)) * factor(a, b);
    //                     for (integer c = b; c < NDIM; ++c) {
    //                         A0() -= m0(a, b, c) * D(a, b, c) * (real(1) / real(6)) * factor(a, b, c);
    //                     }
    //                 }
    //             }
                A0[0] = m0[0] * D[0];
                if (type != RHO) {
                    for (integer i = taylor_sizes[0]; i != taylor_sizes[1]; ++i) {
                        A0[0] -= m0[i] * D[i];
                    }
                }
                for (integer i = taylor_sizes[1]; i != taylor_sizes[2]; ++i) {
                    A0[0] += m0[i] * D[i] * (factor[i] / real(2));
                }
                for (integer i = taylor_sizes[2]; i != taylor_sizes[3]; ++i) {
                    A0[0] -= m0[i] * D[i] * (factor[i] / real(6));
                }

                for (integer a = 0; a < NDIM; ++a) {
                    A0(a) = m0() * D(a);
                    for (integer b = 0; b < NDIM; ++b) {
                        if (type != RHO) {
                            A0(a) -= m0(a) * D(a, b);
                        }
                        for (integer c = b; c < NDIM; ++c) {
                            const auto tmp = D(a, b, c) * (factor(c, b) / real(2));
                            A0(a) += m0(c, b) * tmp;
                        }
                    }
                }

                if (type == RHO) {
                    for (integer a = 0; a < NDIM; ++a) {
                        for (integer b = 0; b < NDIM; ++b) {
                            for (integer c = b; c < NDIM; ++c) {
                                for (integer d = c; d < NDIM; ++d) {
                                    const auto tmp = D(a, b, c, d) * (factor(b, c, d) / real(6));
                                    B0[a] -= n0(b, c, d) * tmp;
                                }
                            }
                        }
                    }
                }

                for (integer a = 0; a < NDIM; ++a) {
                    for (integer b = a; b < NDIM; ++b) {
                        A0(a, b) = m0() * D(a, b);
                        for (integer c = 0; c < NDIM; ++c) {
                            A0(a, b) -= m0(c) * D(a, b, c);
                        }
                    }
                }

    //             for (integer a = 0; a < NDIM; ++a) {
    //                 for (integer b = a; b < NDIM; ++b) {
    //                     for (integer c = b; c < NDIM; ++c) {
    //                         A0(a, b, c) = m0() * D(a, b, c);
    //                     }
    //                 }
    //             }
                for (integer i = taylor_sizes[2]; i != taylor_sizes[3]; ++i) {
                    A0[i] =  m0[0] * D[i];
                }

                std::lock_guard<hpx::lcos::local::spinlock> lock(*L_mtx);
                for (integer i = 0; i != simd_len && i + li < list_size; ++i) {
                    const integer iii0 = bnd.first[li + i];
                    expansion& Liii0 = L[iii0];

                    // FIXME: a similar loop below goes over the range [0, 4)
                    //        which is correct?
                    // DCM: Both
                    for (integer j = 0; j != 20; ++j) {
                        Liii0[j] += A0[j][i];
                    }

                    if (type == RHO) {
                        space_vector& L_ciii0 = L_c[iii0];
                        for (integer j = 0; j != NDIM; ++j) {
                            L_ciii0[j] += B0[j][i];
                        }
                    }
                }
            }
        });
    PROF_END;
}