LambertW function

[Otto_L_]

Does any body know of a library that implements the LambertW function (also known as ProductLog) in FORTRAN?

Thanks in advance,

OL

[mecej4]

The NAG FL library contains two such routines: C05BAF and C05BBF. For example, for real arguments, see http://www.nag.com/numeric/fl/nagdoc_fl25/html/c05/c05baf.html .

The GNU Scientific Library (GSL) also contains the Lambert W function. You many have to write a wrapper function to call it from Fortran. This worked for the Cygwin-64 version of GSL with GFortran on Windows:

#include <gsl/gsl_sf_lambert.h>
double lambertw_(double *x){
return gsl_sf_lambert_W0(*x);
}

[Otto_L_]

Thanks.

Unfortunately, this is for proprietary work. The GNU will not do. I have to see if my organization has the NAG.

Best,

OL.

[mecej4]

If you do not have NAG, you may use the TOMS 743 source code, available at Netlib/TOMS. An F90 version is provided by John Burkardt at http://people.sc.fsu.edu/~jburkardt/f_src/toms743/toms743.html, subject to the LGPL license terms.

It would be fairly easy to write your own code to solve the defining equation using Newton-Raphson. For starting values, you could use the truncated series expansion or a Pade' approximation to the series.

[Otto_L_]

Thanks mecej4 for the reply. Sorry for not replying earlier, I have been out on vacation :) Now it's over :(

I think I could try what you propose. However, I am trying to improve performance on an existing solution that used something like Newton's method, but that did not recast the original problem into one of computing the LambertW function. I am concerned that I may end up coding the same thing I already have in just a slightly different way. Basically, in order for this to make a difference, I would have to come up with (or find) a really good (fast and accurate) implementation of this function. There are series expansions around which are probably faster than a search by iterations (every time one does an iteration one needs to evaluate the Exp function, which itself uses a series expansion). Fortunately, I think I only need the 'easy' branch (W_0).

Thanks,

OL.

[jimdempseyatthecove]

>> I am trying to improve performance on an existing solution that used something like Newton's method

Newton's method is typically implemented as an iterative convergence process whereby a range is split, the new split point is evaluated, and one of the two sub ranges is selected for the next split. The process repeats until the sub range is the size of the desired precision.

This works effectively as a scalar method. This method may be modified slightly to use SIMD vectorization to make the split into three (2 interior points), five (4 interior points), nine (8 interior points), seventeen (16 interior points).

Jim Dempsey

[mecej4]

Jim's description matches bisection, not Newton's method. In Newton's method, once the vicinity of the root is reached the convergence is quadratic.

There are hybrid methods, such as Brent's, in which the initial steps reduce the interval, followed by a switch-over to Newton-Raphson. There is no tracking of a root-bracket in Newton-Raphson. If a 4-digit approximation (series, Pade', etc.) is used to obtain the first estimate, two iterations of Newton's method will yield full double precision in the result. Many of the special function routines in TOMS (e.g., TOMS-443) use a fixed number of steps of Newton's method.

[TimP]

GSL lambert refers to https://cs.uwaterloo.ca/research/tr/1993/03/W.pdf

Apparently, it uses a Maple Pade curve fit.  The Horner polynomial evaluation might be speeded up by using a few separate terms to take advantage of the pipeline of most CPUs introduced since it was written.  In case a benefit is seen, it shouldn't be difficult to copy the C in Fortran.

[Greg_T_]

Hi Otto,

The NIST Guide to Available Mathematical Software (GAMS) http://gams.nist.gov/ has a large library of functions.  I've used a few routines from GAMS successfully.  Their search tool may help narrow the list to the function you're looking for.

Regards,
Greg