Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted
##

Otto_L_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-25-2016
11:10 PM

77 Views

LambertW function

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

Thanks in advance,

OL

8 Replies

Highlighted
##

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-26-2016
04:40 AM

77 Views

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); }

Highlighted
##

Otto_L_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-26-2016
06:00 AM

77 Views

Thanks.

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

Best,

OL.

Highlighted
##

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

05-26-2016
07:32 AM

77 Views

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.

Highlighted
##

Otto_L_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-13-2016
07:23 AM

77 Views

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.

Highlighted
##

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-13-2016
07:36 AM

77 Views

>> 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

Highlighted
##

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-13-2016
10:02 AM

77 Views

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.

Highlighted
##

TimP

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-13-2016
11:01 AM

77 Views

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.

Highlighted
##

Greg_T_

New Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
06:49 AM

77 Views

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

For more complete information about compiler optimizations, see our Optimization Notice.