From: wolf Date: Tue, 15 Oct 2002 16:13:27 +0000 (+0000) Subject: Move legendre.cc into polynomial.cc to get back to the unique mapping between .h... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0bba129854bfa7d819cfa1e4e2141c6cbddfe9eb;p=dealii-svn.git Move legendre.cc into polynomial.cc to get back to the unique mapping between .h and .cc files. git-svn-id: https://svn.dealii.org/trunk@6653 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/source/legendre.cc b/deal.II/base/source/legendre.cc deleted file mode 100644 index 4b059b830a..0000000000 --- a/deal.II/base/source/legendre.cc +++ /dev/null @@ -1,211 +0,0 @@ -//-------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2000, 2001, 2002 by the deal.II authors -// -// This file is subject to QPL and may not be distributed -// without copyright and license information. Please refer -// to the file deal.II/doc/license.html for the text and -// further information on this license. -// -//-------------------------------------------------------------------- - - -#include -#include - -#include - -//TODO:[?] This class leaks memory, but only at the very end of a program. -// Since it expands the Legendre::coefficients array, the elements -// of this static variable are not destroyed at the end of the program -// run. While this is not a problem (since the returned memory could -// not be used anyway then), it is a little confusing when looking at -// a memory checked such as "purify". Maybe, this should be handled somehow -// to avoid this confusion in future. - -// Reserve space for polynomials up to degree 19. Should be sufficient -// for the start. -template -std::vector *> -Legendre::recursive_coefficients( - 20, static_cast*>(0)); -template -std::vector *> -Legendre::shifted_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; -}; - - -// See polynomial.cc for an expalanation of this - -#if ((__GNUC__ == 3) && (__GNUC_MINOR__ == 1)) -#define SHIFT_TYPE double -#else -#define SHIFT_TYPE long double -#endif - -template -void -Legendre::compute_coefficients (const unsigned int k_) -{ - unsigned int k = k_; - - // first make sure that no other - // thread intercepts the operation - // of this function - coefficients_lock.acquire (); - - // 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.size() >= k+1) && - (recursive_coefficients[k] == 0))) - // no, then generate the - // respective coefficients - { - recursive_coefficients.resize (k+1, 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 - recursive_coefficients[0] = c0; - recursive_coefficients[1] = c1; - // Compute polynomials - // orthogonal on [0,1] - c0 = new std::vector(*c0); - c1 = new std::vector(*c1); - - Polynomial::shift(*c0, (SHIFT_TYPE) -1.); - Polynomial::scale(*c0, 2.); - Polynomial::shift(*c1, (SHIFT_TYPE) -1.); - Polynomial::scale(*c1, 2.); - Polynomial::multiply(*c1, std::sqrt(3.)); - shifted_coefficients[0]=c0; - shifted_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); - const number b = a*(2*k-1); - const number c = a*(k-1); - - (*ck)[k] = b*(*recursive_coefficients[k-1])[k-1]; - (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2]; - for (unsigned int i=1 ; i<= k-2 ; ++i) - (*ck)[i] = b*(*recursive_coefficients[k-1])[i-1] - -c*(*recursive_coefficients[k-2])[i]; - - (*ck)[0] = -c*(*recursive_coefficients[k-2])[0]; - - // finally assign the newly - // created vector to the - // const pointer in the - // coefficients array - recursive_coefficients[k] = ck; - // and compute the - // coefficients for [0,1] - ck = new std::vector(*ck); - shift(*ck,(SHIFT_TYPE) -1.); - Polynomial::scale(*ck, 2.); - Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); - shifted_coefficients[k] = ck; - }; - }; - - // now, everything is done, so - // release the lock again - coefficients_lock.release (); -} - - - -template -const std::vector & -Legendre::get_coefficients (const unsigned int k) -{ - // 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 - coefficients_lock.acquire (); - const std::vector *p = shifted_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 (const unsigned int k) - : - Polynomial (get_coefficients(k)) -{} - - - -template -std::vector > -Legendre::generate_complete_basis (const unsigned int degree) -{ - std::vector > v; - v.reserve(degree+1); - for (unsigned int i=0; i<=degree; ++i) - v.push_back (Legendre(i)); - return v; -}; - - - -// explicit instantiations -template class Legendre; diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index bb1c2d00d9..bb12b4452c 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -14,6 +14,27 @@ #include #include +#include + +#include + + +// 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 +{ + Threads::ThreadMutex coefficients_lock; +}; + + +// -------------------- class Polynomial ---------------- // + template Polynomial::Polynomial (const std::vector &a): @@ -219,9 +240,9 @@ Polynomial::print(std::ostream& out) const } } -// ------------------------------------------------------------ // +// ------------------ class LagrangeEquidistant --------------- // LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, const unsigned int support_point): @@ -473,6 +494,191 @@ generate_complete_basis (const unsigned int degree) +// ------------------ class Legendre --------------- // + + +//TODO:[?] This class leaks memory, but only at the very end of a program. +// Since it expands the Legendre::coefficients array, the elements +// of this static variable are not destroyed at the end of the program +// run. While this is not a problem (since the returned memory could +// not be used anyway then), it is a little confusing when looking at +// a memory checker such as "purify". Maybe, this should be handled somehow +// to avoid this confusion in future. + +// Reserve space for polynomials up to degree 19. Should be sufficient +// for the start. +template +std::vector *> +Legendre::recursive_coefficients( + 20, static_cast*>(0)); +template +std::vector *> +Legendre::shifted_coefficients( + 20, static_cast*>(0)); + + + +//TODO[WB]: Treat the same way as above +#if ((__GNUC__ == 3) && (__GNUC_MINOR__ == 1)) +#define SHIFT_TYPE double +#else +#define SHIFT_TYPE long double +#endif + +template +void +Legendre::compute_coefficients (const unsigned int k_) +{ + unsigned int k = k_; + + // first make sure that no other + // thread intercepts the operation + // of this function + coefficients_lock.acquire (); + + // 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.size() >= k+1) && + (recursive_coefficients[k] == 0))) + // no, then generate the + // respective coefficients + { + recursive_coefficients.resize (k+1, 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 + recursive_coefficients[0] = c0; + recursive_coefficients[1] = c1; + // Compute polynomials + // orthogonal on [0,1] + c0 = new std::vector(*c0); + c1 = new std::vector(*c1); + + Polynomial::shift(*c0, (SHIFT_TYPE) -1.); + Polynomial::scale(*c0, 2.); + Polynomial::shift(*c1, (SHIFT_TYPE) -1.); + Polynomial::scale(*c1, 2.); + Polynomial::multiply(*c1, std::sqrt(3.)); + shifted_coefficients[0]=c0; + shifted_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); + const number b = a*(2*k-1); + const number c = a*(k-1); + + (*ck)[k] = b*(*recursive_coefficients[k-1])[k-1]; + (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2]; + for (unsigned int i=1 ; i<= k-2 ; ++i) + (*ck)[i] = b*(*recursive_coefficients[k-1])[i-1] + -c*(*recursive_coefficients[k-2])[i]; + + (*ck)[0] = -c*(*recursive_coefficients[k-2])[0]; + + // finally assign the newly + // created vector to the + // const pointer in the + // coefficients array + recursive_coefficients[k] = ck; + // and compute the + // coefficients for [0,1] + ck = new std::vector(*ck); + shift(*ck,(SHIFT_TYPE) -1.); + Polynomial::scale(*ck, 2.); + Polynomial::multiply(*ck, std::sqrt(2.*k+1.)); + shifted_coefficients[k] = ck; + }; + }; + + // now, everything is done, so + // release the lock again + coefficients_lock.release (); +} + + + +template +const std::vector & +Legendre::get_coefficients (const unsigned int k) +{ + // 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 + coefficients_lock.acquire (); + const std::vector *p = shifted_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 (const unsigned int k) + : + Polynomial (get_coefficients(k)) +{} + + + +template +std::vector > +Legendre::generate_complete_basis (const unsigned int degree) +{ + std::vector > v; + v.reserve(degree+1); + for (unsigned int i=0; i<=degree; ++i) + v.push_back (Legendre(i)); + return v; +}; + + +// ------------------ explicit instantiations --------------- // + template class Polynomial; template class Polynomial; template class Polynomial; @@ -483,3 +689,5 @@ template void Polynomial::shift(const double offset); template void Polynomial::shift(const long double offset); template void Polynomial::shift(const long double offset); template void Polynomial::shift(const long double offset); + +template class Legendre;