From 2716ff8dfe219411c27cf672665aa6ea2c4b6b30 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 24 Oct 2023 11:04:00 -0600 Subject: [PATCH] Use a shared mutex instead of a unique one. --- source/base/polynomial.cc | 267 +++++++++++++++++++------------------- 1 file changed, 133 insertions(+), 134 deletions(-) diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index dcac24e43c..4dba286918 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -25,26 +25,12 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN -// have a lock that guarantees that at most one thread is changing and -// accessing the @p{coefficients} arrays of classes implementing -// polynomials with tables. make this lock local to this file. -// -// having only one lock for all of these classes is probably not going -// to be a problem since we only need it on very rare occasions. if -// someone finds this is a bottleneck, feel free to replace it by a -// more fine-grained solution -namespace -{ - std::mutex coefficients_lock; -} - - - namespace Polynomials { // -------------------- class Polynomial ---------------- // @@ -861,6 +847,12 @@ namespace Polynomials // ------------------ class Hierarchical --------------- // + // have a lock that guarantees that at most one thread is changing and + // accessing the `recursive_coefficients array. + namespace + { + std::shared_mutex coefficients_lock; + } // Reserve space for polynomials up to degree 19. Should be sufficient // for the start. @@ -879,131 +871,139 @@ namespace Polynomials Hierarchical::compute_coefficients(const unsigned int k_) { unsigned int k = k_; - - // first make sure that no other - // thread intercepts the operation - // of this function - // for this, acquire the lock - // until we quit this function - std::lock_guard lock(coefficients_lock); - // The first 2 coefficients // are hard-coded if (k == 0) k = 1; - // check: does the information - // already exist? - if ((recursive_coefficients.size() < k + 1) || - (recursive_coefficients[k].get() == nullptr)) - // no, then generate the - // respective coefficients + + // First see whether the coefficients we need have already been + // computed. This is a read operation, and so we can do that + // with a shared lock. + // + // (We could have gotten away without any lock at all if the + // inner pointers were std::atomic>, but + // first, there is no such specialization of std::atomic that + // is mutex-free, and then there is also the issue that the + // outer vector may be resized and that can definitely not + // be guarded against without a mutex of some sort.) + { + std::shared_lock lock(coefficients_lock); + + if ((recursive_coefficients.size() >= k + 1) && + (recursive_coefficients[k].get() != nullptr)) + return; + } + + // Having gotten here, we know that we need to compute a new set + // of coefficients. This has to happen under a unique lock because + // we're not only reading, but writing into the data structures: + std::unique_lock lock(coefficients_lock); + + // First make sure that there is enough + // space in the array for the + // coefficients, so we have to resize + // it to size k+1 + + // but it's more complicated than + // that: we call this function + // recursively, so if we simply + // resize it to k+1 here, then + // compute the coefficients for + // degree k-1 by calling this + // function recursively, then it will + // reset the size to k -- not enough + // for what we want to do below. the + // solution therefore is to only + // resize the size if we are going to + // *increase* it + if (recursive_coefficients.size() < k + 1) + recursive_coefficients.resize(k + 1); + + if (k <= 1) { - // make sure that there is enough - // space in the array for the - // coefficients, so we have to resize - // it to size k+1 - - // but it's more complicated than - // that: we call this function - // recursively, so if we simply - // resize it to k+1 here, then - // compute the coefficients for - // degree k-1 by calling this - // function recursively, then it will - // reset the size to k -- not enough - // for what we want to do below. the - // solution therefore is to only - // resize the size if we are going to - // *increase* it - if (recursive_coefficients.size() < k + 1) - recursive_coefficients.resize(k + 1); - - if (k <= 1) - { - // create coefficients - // vectors for k=0 and k=1 - // - // allocate the respective - // amount of memory and - // later assign it to the - // coefficients array to - // make it const - std::vector c0(2); - c0[0] = 1.; - c0[1] = -1.; - - std::vector c1(2); - c1[0] = 0.; - c1[1] = 1.; - - // now make these arrays - // const - recursive_coefficients[0] = - std::make_unique>(std::move(c0)); - recursive_coefficients[1] = - std::make_unique>(std::move(c1)); - } - else if (k == 2) - { - coefficients_lock.unlock(); - compute_coefficients(1); - coefficients_lock.lock(); + // create coefficients + // vectors for k=0 and k=1 + // + // allocate the respective + // amount of memory and + // later assign it to the + // coefficients array to + // make it const + std::vector c0(2); + c0[0] = 1.; + c0[1] = -1.; + + std::vector c1(2); + c1[0] = 0.; + c1[1] = 1.; + + // now make these arrays + // const + recursive_coefficients[0] = + std::make_unique>(std::move(c0)); + recursive_coefficients[1] = + std::make_unique>(std::move(c1)); + } + else if (k == 2) + { + coefficients_lock.unlock(); + compute_coefficients(1); + coefficients_lock.lock(); - std::vector c2(3); + std::vector c2(3); - const double a = 1.; // 1./8.; + const double a = 1.; // 1./8.; - c2[0] = 0. * a; - c2[1] = -4. * a; - c2[2] = 4. * a; + c2[0] = 0. * a; + c2[1] = -4. * a; + c2[2] = 4. * a; - recursive_coefficients[2] = - std::make_unique>(std::move(c2)); - } - else + recursive_coefficients[2] = + std::make_unique>(std::move(c2)); + } + else + { + // for larger numbers, + // compute the coefficients + // recursively. to do so, + // we have to release the + // lock temporarily to + // allow the called + // function to acquire it + // itself + coefficients_lock.unlock(); + compute_coefficients(k - 1); + coefficients_lock.lock(); + + std::vector ck(k + 1); + + const double a = 1.; // 1./(2.*k); + + ck[0] = -a * (*recursive_coefficients[k - 1])[0]; + + for (unsigned int i = 1; i <= k - 1; ++i) + ck[i] = a * (2. * (*recursive_coefficients[k - 1])[i - 1] - + (*recursive_coefficients[k - 1])[i]); + + ck[k] = a * 2. * (*recursive_coefficients[k - 1])[k - 1]; + // for even degrees, we need + // to add a multiple of + // basis fcn phi_2 + if ((k % 2) == 0) { - // for larger numbers, - // compute the coefficients - // recursively. to do so, - // we have to release the - // lock temporarily to - // allow the called - // function to acquire it - // itself - coefficients_lock.unlock(); - compute_coefficients(k - 1); - coefficients_lock.lock(); - - std::vector ck(k + 1); - - const double a = 1.; // 1./(2.*k); - - ck[0] = -a * (*recursive_coefficients[k - 1])[0]; - - for (unsigned int i = 1; i <= k - 1; ++i) - ck[i] = a * (2. * (*recursive_coefficients[k - 1])[i - 1] - - (*recursive_coefficients[k - 1])[i]); - - ck[k] = a * 2. * (*recursive_coefficients[k - 1])[k - 1]; - // for even degrees, we need - // to add a multiple of - // basis fcn phi_2 - if ((k % 2) == 0) - { - double b = 1.; // 8.; - // for (unsigned int i=1; i<=k; ++i) - // b /= 2.*i; + double b = 1.; // 8.; + // for (unsigned int i=1; i<=k; ++i) + // b /= 2.*i; - ck[1] += b * (*recursive_coefficients[2])[1]; - ck[2] += b * (*recursive_coefficients[2])[2]; - } - // finally assign the newly - // created vector to the - // const pointer in the - // coefficients array - recursive_coefficients[k] = - std::make_unique>(std::move(ck)); + ck[1] += b * (*recursive_coefficients[2])[1]; + ck[2] += b * (*recursive_coefficients[2])[2]; } + // finally assign the newly + // created vector to the + // const pointer in the + // coefficients array + recursive_coefficients[k] = + std::make_unique>(std::move(ck)); } } @@ -1012,14 +1012,13 @@ namespace Polynomials const std::vector & Hierarchical::get_coefficients(const unsigned int k) { - // first make sure the coefficients - // get computed if so necessary + // First make sure the coefficients get computed if so necessary compute_coefficients(k); - // then get a pointer to the array - // of coefficients. do that in a MT - // safe way - std::lock_guard lock(coefficients_lock); + // Then get a pointer to the array of coefficients. Do that in a MT + // safe way, but since we're only reading information we can do + // that with a shared lock + std::shared_lock lock(coefficients_lock); return *recursive_coefficients[k]; } -- 2.39.5