dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Anyone using or interested in the ODE solvers should take a look
below.

On Thu, Jun 02, 2011 at 02:17:17PM +0200, Benjamin Kehlet wrote:
> On 2 June 2011 14:02, Anders Logg <logg@xxxxxxxxx> wrote:
> > On Thu, Jun 02, 2011 at 01:10:01PM +0200, 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.
> >> >
> >> >
> >> >
> >> >> 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.
> >
> > If it's not possible to make it work, we need to revert back.
>
> I don't know of any solution to this. This is the same problem that we
> discussed some months back (then related to Armadillo): Templated
> libraries which rely on non-templated  code (often old and implemented
> i c), so they only support the types which these underlying libraries
> can handle. I think the only solution here is a change in
> boost::math::Legendre.
>
> Of course another solution would be to split the ODE solver from
> Dolfin and let it continue as a separate project, and then import code
> from that when we are going to look at automation/generating code for
> time-dependent problems.

Yes, perhaps it's time for that. Since it is going to be removed soon
(and replaced by code generation), the best option might be to remove
it before the release of 1.0.

Are there any objections? Is anyone using the ODE solvers?

(They will make a comeback later in new form.)

--
Anders