]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use a shared mutex instead of a unique one.
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 24 Oct 2023 17:04:00 +0000 (11:04 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 24 Oct 2023 17:04:00 +0000 (11:04 -0600)
source/base/polynomial.cc

index dcac24e43cf1d9bade21f77be74b8c1c137dbf88..4dba2869184d00bfc4b341de0e2be660340c6441 100644 (file)
 #include <cmath>
 #include <limits>
 #include <memory>
+#include <shared_mutex>
 
 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<std::mutex> 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<std::unique_ptr<...>>, 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<std::shared_mutex> 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<std::shared_mutex> 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<double> c0(2);
-            c0[0] = 1.;
-            c0[1] = -1.;
-
-            std::vector<double> c1(2);
-            c1[0] = 0.;
-            c1[1] = 1.;
-
-            // now make these arrays
-            // const
-            recursive_coefficients[0] =
-              std::make_unique<const std::vector<double>>(std::move(c0));
-            recursive_coefficients[1] =
-              std::make_unique<const std::vector<double>>(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<double> c0(2);
+        c0[0] = 1.;
+        c0[1] = -1.;
+
+        std::vector<double> c1(2);
+        c1[0] = 0.;
+        c1[1] = 1.;
+
+        // now make these arrays
+        // const
+        recursive_coefficients[0] =
+          std::make_unique<const std::vector<double>>(std::move(c0));
+        recursive_coefficients[1] =
+          std::make_unique<const std::vector<double>>(std::move(c1));
+      }
+    else if (k == 2)
+      {
+        coefficients_lock.unlock();
+        compute_coefficients(1);
+        coefficients_lock.lock();
 
-            std::vector<double> c2(3);
+        std::vector<double> 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<const std::vector<double>>(std::move(c2));
-          }
-        else
+        recursive_coefficients[2] =
+          std::make_unique<const std::vector<double>>(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<double> 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<double> 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<const std::vector<double>>(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<const std::vector<double>>(std::move(ck));
       }
   }
 
@@ -1012,14 +1012,13 @@ namespace Polynomials
   const std::vector<double> &
   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<std::mutex> 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<std::shared_mutex> lock(coefficients_lock);
     return *recursive_coefficients[k];
   }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.