From: Wolfgang Bangerth Date: Fri, 16 Feb 2001 09:38:25 +0000 (+0000) Subject: Make Legendre class thread-safe. X-Git-Tag: v8.0.0~19721 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=10fb9366150e8f9769cdcfb9cddc4c8fc04c2364;p=dealii.git Make Legendre class thread-safe. git-svn-id: https://svn.dealii.org/trunk@3949 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 20c4bf8298..8939d4eb8b 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -174,31 +174,40 @@ class LagrangeEquidistant: public Polynomial template class Legendre : public Polynomial { -public: - /** - * Constructor for polynomial of - * order @p{k}. - */ - Legendre (unsigned int k); -private: - /** - * Vector with already computed - * coefficients. - */ - static std::vector > 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& 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*> 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 & + get_coefficients (const unsigned int k); }; diff --git a/deal.II/base/source/legendre.cc b/deal.II/base/source/legendre.cc index b0cdf7176b..471c57c84d 100644 --- a/deal.II/base/source/legendre.cc +++ b/deal.II/base/source/legendre.cc @@ -13,54 +13,140 @@ #include +#include -// 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 -std::vector > -Legendre::coefficients(20,0); +std::vector*> +Legendre::coefficients(20, + static_cast*>(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 void -Legendre::compute_coefficients (unsigned int k) +Legendre::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 *c0 = new std::vector(1); + (*c0)[0] = 1.; + + std::vector *c1 = new std::vector(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 *ck = new std::vector(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 -const std::vector& -Legendre::get_coefficients (unsigned int k) +const std::vector & +Legendre::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 *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 -Legendre::Legendre (unsigned int k) - : Polynomial (get_coefficients(k)) +Legendre::Legendre (const unsigned int k) + : + Polynomial (get_coefficients(k)) {} + + + +// explicit instantiations +template class Legendre;