]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make Legendre class thread-safe.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 Feb 2001 09:38:25 +0000 (09:38 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 Feb 2001 09:38:25 +0000 (09:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@3949 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h
deal.II/base/source/legendre.cc

index 20c4bf82982fdfec1445d406acf4c26e1693a74b..8939d4eb8b66f717280215c250e2e93d7306d88b 100644 (file)
@@ -174,31 +174,40 @@ class LagrangeEquidistant: public Polynomial<double>
 template <typename number>
 class Legendre : public Polynomial<number>
 {
-public:
-                                  /**
-                                   * Constructor for polynomial of
-                                   * order @p{k}.
-                                   */
-  Legendre (unsigned int k);
-private:
-                                  /**
-                                   * Vector with already computed
-                                   * coefficients.
-                                   */
-  static std::vector<std::vector<number> > coefficients;
-
-                                  /**
-                                   * Compute coefficients recursively.
-                                   */
-  static void compute_coefficients (unsigned int k);
-
-                                  /**
-                                   * Get coefficients for
-                                   * constructor.  This way, it can
-                                   * use the non-standard constructor
-                                   * of @ref{Polynomial}.
-                                   */
-  static const std::vector<number>& get_coefficients (unsigned int k);
+  public:
+                                    /**
+                                     * Constructor for polynomial of
+                                     * order @p{k}.
+                                     */
+    Legendre (const unsigned int k);
+
+  private:
+                                    /**
+                                     * Vector with already computed
+                                     * coefficients. For each degree
+                                     * of the polynomial, we keep one
+                                     * pointer to the list of
+                                     * coefficients; we do so rather
+                                     * than keeping a vector of
+                                     * vectors in order to simplify
+                                     * programming multithread-safe.
+                                     */
+    static std::vector<const std::vector<number>*> coefficients;
+    
+                                    /**
+                                     * Compute coefficients recursively.
+                                     */
+    static void compute_coefficients (unsigned int k);
+    
+                                    /**
+                                     * Get coefficients for
+                                     * constructor.  This way, it can
+                                     * use the non-standard
+                                     * constructor of
+                                     * @ref{Polynomial}.
+                                     */
+    static const std::vector<number> &
+    get_coefficients (const unsigned int k);
 };
 
 
index b0cdf7176bbdbae9e57e56889f7b6da51d11a6df..471c57c84d924c1fd6beeb9ff5c4f2159dbd46cb 100644 (file)
 
 
 #include <base/polynomial.h>
+#include <base/thread_management.h>
 
-// Reserve space for polynomials up to degree 19. Should be sufficient
-// in most cases.
 
+// Reserve space for polynomials up to degree 19. Should be sufficient
+// for the start.
 template <typename number>
-std::vector<std::vector<number> >
-Legendre<number>::coefficients(20,0);
+std::vector<const std::vector<number>*>
+Legendre<number>::coefficients(20,
+                              static_cast<const std::vector<number>*>(0));
+
+
+// have a lock that guarantees that at most one thread is changing and
+// accessing the @p{coefficients} array. make this lock local to this
+// file
+namespace 
+{
+  Threads::ThreadMutex coefficients_lock;
+};
+
 
 
 template <typename number>
 void
-Legendre<number>::compute_coefficients (unsigned int k)
+Legendre<number>::compute_coefficients (const unsigned int k)
 {
-  if (k<=1)
+                                  // first make sure that no other
+                                  // thread intercepts the operation
+                                  // of this function
+  coefficients_lock.acquire ();
+
+                                  // check: does the information
+                                  // already exist?
+  if ((coefficients.size() < k+1) ||
+      ((coefficients.size() >= k+1) &&
+       (coefficients[k] == 0)))
+                                    // no, then generate the
+                                    // respective coefficients
     {
-      coefficients[0].resize(1);
-      coefficients[0][0] = 1.;
-      coefficients[1].resize(2);
-      coefficients[1][0] = 0.;
-      coefficients[1][1] = 1.;
-    } else {
-      compute_coefficients(k-1);
-      coefficients[k].resize(k+1);
-      const double a = 1./k+1;
-      const double b = a*(2*k+1);
+      coefficients.resize (k+1, 0);
       
-      coefficients[k][k]   = b*coefficients[k-1][k-1];
-      coefficients[k][k-1] = b*coefficients[k-1][k-2];
-      for (unsigned int i=1 ; i<= k-2 ; ++i)
-       coefficients[k][i] = b*coefficients[k-1][i-1]
-         - k*a*coefficients[k-2][i];
-      coefficients[k][0]   = -k*a*coefficients[k-2][0];
-    }
+      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<number> *c0 = new std::vector<number>(1);
+         (*c0)[0] = 1.;
+
+         std::vector<number> *c1 = new std::vector<number>(2);
+         (*c1)[0] = 0.;
+         (*c1)[1] = 1.;
+
+                                          // now make these arrays
+                                          // const
+         coefficients[0] = c0;
+         coefficients[1] = c1;
+       }
+      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.release ();
+         compute_coefficients(k-1);
+         coefficients_lock.acquire ();
+
+         std::vector<number> *ck = new std::vector<number>(k+1);
+         
+         const number a = 1./k+1;
+         const number b = a*(2*k+1);
+         
+         (*ck)[k]   = b*(*coefficients[k-1])[k-1];
+         (*ck)[k-1] = b*(*coefficients[k-1])[k-2];
+         for (unsigned int i=1 ; i<= k-2 ; ++i)
+           (*ck)[i] = b*(*coefficients[k-1])[i-1]
+                      - k*a*(*coefficients[k-2])[i];
+         (*ck)[0]   = -k*a*(*coefficients[k-2])[0];
+
+                                          // finally assign the newly
+                                          // created vector to the
+                                          // const pointer in the
+                                          // coefficients array
+         coefficients[k] = ck;
+       };
+    };
+
+                                  // now, everything is done, so
+                                  // release the lock again
+  coefficients_lock.release ();
 }
 
 
 
 template <typename number>
-const std::vector<number>&
-Legendre<number>::get_coefficients (unsigned int k)
+const std::vector<number> &
+Legendre<number>::get_coefficients (const unsigned int k)
 {
+                                  // first make sure the coefficients
+                                  // get computed if so necessary
   compute_coefficients (k);
-  return coefficients[k];
+
+                                  // then get a pointer to the array
+                                  // of coefficients. do that in a MT
+                                  // safe way
+  coefficients_lock.acquire ();
+  const vector<number> *p = coefficients[k];
+  coefficients_lock.release ();
+
+                                  // return the object pointed
+                                  // to. since this object does not
+                                  // change any more once computed,
+                                  // this is MT safe
+  return *p;
 }
 
 
 
 template <typename number>
-Legendre<number>::Legendre (unsigned int k)
-  : Polynomial<number> (get_coefficients(k))
+Legendre<number>::Legendre (const unsigned int k)
+               :
+               Polynomial<number> (get_coefficients(k))
 {}
+
+
+
+// explicit instantiations
+template class Legendre<double>;

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.