dolfin team mailing list archive

[Benjamin Kehlet <benjamik@xxxxxxxxx>]

On 2 June 2011 14:59, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>
>
> On 02/06/11 12:10, Benjamin Kehlet wrote:
>> On 2 June 2011 11:51, Anders Logg <logg@xxxxxxxxx> wrote:
>>> On Thu, Jun 02, 2011 at 10:46:29AM +0100, Garth N. Wells wrote:
>>>>
>>>>
>>>> On 02/06/11 10:26, Anders Logg wrote:
>>>>> On Thu, Jun 02, 2011 at 10:07:59AM +0100, Garth N. Wells wrote:
>>>>>>
>>>>>>
>>>>>> On 01/06/11 23:46, Anders Logg wrote:
>>>>>>> Have you checked that there is no performance penalty?
>>>>>>
>>>>>> I just have - evaluating a Legendgre polynomial 10k times at the same
>>>>>> point is just noise with both methods (of the order 10^-5 - 10^-4 s).
>>>>>
>>>>> It may be noise for some applications, but not for others. I'm not
>>>>> sure this is a bottle-neck for the ODE code (Benjamin will know) but
>>>>> we need to evaluate Legendre polynomials of degree > 100 many times
>>>>> and then it may not be noise.
>>>>>
>>>>
>>>> For very high degree (e.g. 200) Boost is marginally faster.
>>>
>>> Sounds promising then.
>>>
>>>>>> The Boost code is slightly slower because it doesn't cache the values
>>>>>> (which is nice not to do), but may be faster if the call is inlined.
>>>>>> It's not possible to inline it at the moment because of clashes between
>>>>>> tr1:tuple and boost::tuple (Boost bug, I suspect). Old and new are the
>>>>>> same when evaluating at different points.
>>>>>
>>>>> Let's wait for Benjamin to comment.
>>>>>
>>>>
>>>> The speed is about the same (with scope to improve the speed for Boost)
>>>> for unique values. The caller should be responsible for caching, if
>>>> desired, since it can lead to memory blow out.
>>>>
>>>> Legendre does not appear in the ode code. It only appears in the
>>>> computation of quadrature schemes.
>>>
>>> True, but the quadrature schemes are used in the ode code.
>>>
>>> --
>>> Anders
>>>
>>>
>>>> Garth
>>>>
>>>>
>>>>
>>>>>> Garth
>>>>>>
>>>>>>
>>>>>>> Benjamin has
>>>>>>> worked quite hard on optimizing some of the basic math routines (in
>>>>>>> some cases by many many orders of magnitude).
>>>>>>>
>>>>>>> Benjamin, can you take a look that it still works?
>>
>> Yes, the performance seems to be about the same, but I'm unable to
>> compile it with support for GMP.
>>
>> /usr/include/boost/math/special_functions/legendre.hpp:178:
>> instantiated from ‘typename boost::math::tools::promote_args<RT,
>> float, float, float, float, float>::type boost::math::legendre_p(int,
>> int, T, const Policy&) [with T = __gmp_expr<__mpf_struct [1],
>> __mpf_struct [1]>, Policy =
>> boost::math::policies::policy<boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy,
>> boost::math::policies::default_policy>]’
>> /usr/include/boost/math/special_functions/legendre.hpp:185:
>> instantiated from ‘typename boost::math::tools::promote_args<RT,
>> float, float, float, float, float>::type boost::math::legendre_p(int,
>> int, T) [with T = __gmp_expr<__mpf_struct [1], __mpf_struct [1]>]’
>> /home/benjamik/fenics/dolfin-wells_gmp/dolfin/math/Legendre.cpp:42:
>> instantiated from here
>> /usr/include/boost/math/special_functions/legendre.hpp:167: error: no
>> matching function for call to ‘pow(__gmp_expr<__mpf_struct [1],
>> __gmp_binary_expr<long int, __gmp_expr<__mpf_struct [1],
>> __gmp_binary_expr<__gmp_expr<__mpf_struct [1], __mpf_struct [1]>,
>> __gmp_expr<__mpf_struct [1], __mpf_struct [1]>,
>> __gmp_binary_multiplies> >, __gmp_binary_minus> >,
>> __gmp_expr<__mpf_struct [1], __gmp_binary_expr<__gmp_expr<__mpf_struct
>> [1], __mpf_struct [1]>, long int, __gmp_binary_divides> >)’
>> /usr/include/bits/mathcalls.h:154: note: candidates are: double
>> pow(double, double)
>> /usr/include/c++/4.4/cmath:358: note:                 float
>> std::pow(float, float)
>> /usr/include/c++/4.4/cmath:362: note:                 long double
>> std::pow(long double, long double)
>> /usr/include/c++/4.4/cmath:369: note:                 double
>> std::pow(double, int)
>> /usr/include/c++/4.4/cmath:373: note:                 float std::pow(float, int)
>> /usr/include/c++/4.4/cmath:377: note:                 long double
>> std::pow(long double, int)
>> [...]
>>
>> boost::math::legendre seems to rely on std::pow which is not
>> templated, only implemented with the most common types.
>>
>
> Looks like some tweaks are required to work with GMP:
>
> http://www.boost.org/doc/libs/1_43_0/libs/math/doc/sf_and_dist/html/math_toolkit/using_udt/use_mpfr.html

That's not a bad solution, but it requires changing the
multi-precision type from mpf (provided by GMP) to mpfr (which is a
library that extends the floating point functionality in GMP). For
floating-point arithmetic MPFR is much better than pure GMP. I think
CGAL depends on MPFR, so it wouldn't even introduce new dependencies.
The problem is that MPFR doesn't ship with a C++-wrapper (as opposed
to GMP). Although several independent wrappers exists, none of them
are avalilable in Debian/Ubuntu through apt. The one Boost requires is
not updated since 2008 (MPFR has gone from version 2.3 to 3.0.1 since
then).

(Another option would be to take the same approach as Boost ourself:
Implement the few functions that are required (pow() plus possibly a
few more) and place it in the global namespace before including
boost::mat::legendre), but GMP does not provide pow() when the
exponent is a floating point number, so this is not straight forward
without switching to MPFR).

So I guess the question is whether we want to switch to MPFR now, to
get rid of the few lines of code in Legendre.cpp (which performs
reasonably well), when the code is likely to be thrown out pretty soon
anyway. I vote for "no", but I have no problems with moving the entire
ODE solvers to a separate project, then adding it back (without
supporting extended precision) later in the form of code generation
for time dependent problems.

Benjamin

>
> Garth
>
>
>> Benjamin
>>
>>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> Mailing list: https://launchpad.net/~dolfin
>>>>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>>>> More help   : https://help.launchpad.net/ListHelp
>>>
>>> _______________________________________________
>>> Mailing list: https://launchpad.net/~dolfin
>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help   : https://help.launchpad.net/ListHelp
>>>
>
>