Real precision

[jacobwilliams]

trafficstars

Currently, we are using double precision. Consider having the option of:

1. exporting multiple precisions
2. specify the precision via preprocessor flag

For 2, I usually would do something like this:

#ifdef REAL32
    integer,parameter,public :: minpack_rk = real32   !! real kind used by this module [4 bytes]
#elif REAL64
    integer,parameter,public :: minpack_rk = real64   !! real kind used by this module [8 bytes]
#elif REAL128
    integer,parameter,public :: minpack_rk = real128  !! real kind used by this module [16 bytes]
#else
    integer,parameter,public :: minpack_rk = real64   !! real kind used by this module [8 bytes]
#endif

    integer,parameter,private :: wp = minpack_rk  !! local copy of `minpack_rk` with a shorter name

So, if you do nothing, you get double precision. The module also doesn't export wp, since that is too short a variable name to be exporting (if you want it, you will get minpack_rk).

For 1: Look at what I did for quadpack. Not sure we need that for minpack though. Are there use cases for having multiple versions in the same project with different precisions?

[awvwgk]

What would be the strategy for the C API if we support multiple precisions? We can overload multiple procedure names with _Generic and define type-generic macros to dispatch to the correct symbols, but having one symbol with different precisions will make it more difficult to reliably define a C API.

[certik]

I see. So the idea is to keep things wp and later figure out if we want things to work with multiple precisions.

[jacobwilliams]

Look at this file in quadpack. This is the highest level module, and you see I am renaming all the routines, so there is no name collisions. So the C API would use this high-level module and provide wrappers for it I guess?

[awvwgk]

For the C bindings we might want to export both single and double precision versions, if those are available in the library:

typedef void (*minpack_sfunc)(
    int /* n */,
    const float* /* x */,
    float* /* fvec */,
    int* /* iflag */,
    void* /* udata */);

typedef void (*minpack_dfunc)(
    int /* n */,
    const double* /* x */,
    double* /* fvec */,
    int* /* iflag */,
    void* /* udata */);

MINPACK_EXTERN void MINPACK_CALL
minpack_shybrd1(
    minpack_sfunc /* fcn */,
    int /* n */,
    float* /* x */,
    float* /* fvec */,
    float /* tol */,
    int* /* info */,
    float* /* wa */,
    int /* lwa */,
    void* /* udata */);

MINPACK_EXTERN void MINPACK_CALL
minpack_dhybrd1(
    minpack_dfunc /* fcn */,
    int /* n */,
    double* /* x */,
    double* /* fvec */,
    double /* tol */,
    int* /* info */,
    double* /* wa */,
    int /* lwa */,
    void* /* udata */);

/* Dispatch based on signature of callback */
#define minpack_hybrd1(fcn, ...) \
    _Generic((fcn), \
             minpack_sfunc: minpack_shybrd1, \
             minpack_dfunc: minpack_dhybrd1) \
            (fcn, VAR_ARGS)

Consuming C code could use the type-generic macro to access the correct function

minpack_hybrd(fcn, n, x, fvec, tol, &info, wa, lwa, NULL);

[ivan-pi]

I recall reading (perhaps in the MINPACK User Guide) that double precision is practically a must for difficult optimization problems. How is it with other popular CSE languages like MATLAB, Python, R, and Julia, that would be wrapping MINPACK? I believe the first three use double precision by default.

On the other hand, I can imagine applications in image processing/computer vision where it would be desirable to solve thousands of optimization problems on a GPU, and single precision might suffice. For this scenario I like @awvwgk's proposition of exporting both single and double versions.

[ivan-pi]

Does generic selection (C11) work in C++ mode too? I couldn't find an answer so I assume it's no.

In C++ we would need to use function overloading or template function specialization instead:

/* minpack.hpp */

void minpack_hybrd1(minpack_sfunc, float, ...) 
{ 
  // call minpack_shybrd1 
}
void minpack_hybrd1(minpack_dfunc, double, ...) 
{ 
  // call minpack_dhybrd1
}