sundials generic math

View: New views
1 Messages — Rating Filter:   Alert me  

sundials generic math

by David Lorenzetti :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hello,

Having just started with CVODE (literally ran config-make-make-install
this morning), I ran into something that looks a little suspect in the
use of macro SUNDIALS_USE_GENERIC_MATH.  Can somebody comment on this?

Header file "sundials_config.h" says this:

 > /* Use generic math functions
 >  * If it was decided that generic math functions can be used, then
 >  *     #define SUNDIALS_USE_GENERIC_MATH 1
 >  * otherwise
 >  *     #define SUNDIALS_USE_GENERIC_MATH 0
 >  */
 > #define SUNDIALS_USE_GENERIC_MATH 1

This seems to indicate that preprocessor tests of
SUNDIALS_USE_GENERIC_MATH ought to concern its value (0 or 1), not
whether it's defined.

With that in mind, look at source file "sundials_math.c", which
contains a few blocks of code like this:

 > #if defined(SUNDIALS_USE_GENERIC_MATH)
 >   return((realtype) pow((double) base, (double) exponent));
 > #elif defined(SUNDIALS_DOUBLE_PRECISION)
 >   return(pow(base, exponent));
 > #elif defined(SUNDIALS_SINGLE_PRECISION)
 >   return(powf(base, exponent));
 > #elif defined(SUNDIALS_EXTENDED_PRECISION)
 >   return(powl(base, exponent));
 > #endif

Here, SUNDIALS_USE_GENERIC_MATH gets tested for existence, not for
value.  Thus the generic math code always will be used.

In other words, suppose the user #defines SUNDIALS_USE_GENERIC_MATH as
0.  Then that macro is defined, and the source code uses the "generic
math" version of the code-- contrary to the user's wishes.

Therefore it seems to me those blocks of code ought to look something
like this:

 > #if defined(SUNDIALS_DOUBLE_PRECISION)
 >   return(pow(base, exponent));
 > #elif(SUNDIALS_USE_GENERIC_MATH == 1)
 >   return((realtype) pow((double) base, (double) exponent));
 > #elif defined(SUNDIALS_SINGLE_PRECISION)
 >   return(powf(base, exponent));
 > #elif defined(SUNDIALS_EXTENDED_PRECISION)
 >   return(powl(base, exponent));
 > #endif

This form has two advantages:

**  It tests SUNDIALS_USE_GENERIC_MATH against 1, rather than for
whether it's defined.  This avoids the problem of always using the
generic math code, regardless of the user's intent.

**  It places the SUNDIALS_DOUBLE_PRECISION piece first, thus avoiding
the unnecessary coercion of a double to double in cases where
SUNDIALS_USE_GENERIC_MATH is set to 1.  While an optimizing compiler
probably removes the unneeded coercion, there's no reason to require
that it does (and there's no reason to introduce another possible
cause of differences between debug and release versions of the code).

So-- comments?  Have I missed something here?

-Dave

David Lorenzetti
Lawrence Berkeley National Laboratory

----------------------------------------------------------------------------------------------------------------
Center for Applied Scientific Computing
Lawrence Livermore National Laboratory
P.O. Box 808, L-561
Livermore, CA  94551

(925) 424-6013  (Office)            cswoodward@...
(925) 423-2993  (Fax)               http://people.llnl.gov/woodward6