Problems with GSL using Legendre polynomials

Question

I'm trying to update an old code that was using a version of GSL with deprecated functions, but I'm having troubles to find how to use the new version of the normalized Legendre polynomials functions. Here's a snippet that summarizes the problem:

#include <iostream>
#include <gsl/gsl_sf_legendre.h>
#include <cmath>

#define GSL_NEW

using namespace std;

int main() {

  int order = 17;
  int ntheta = 36;
  double theta_step = M_PI / ntheta;
  double c, theta;
  double legendre[ntheta][order+1];

  for( int m = 0; m <= order; m += 2) {
    for(int l = m; l <= ntheta; l += 2 ) {
      for( int t = 0; t < ntheta; t++ ) {
        theta = ( ntheta + 0.5 ) * theta_step;
        c = cos(theta);

        if( l == m ) {
#ifdef GSL_NEW
          gsl_sf_legendre_array( GSL_SF_LEGENDRE_SPHARM, order, c, &legendre[t][l] );
          cout << legendre[t][l] << endl;
#else
          gsl_sf_legendre_sphPlm_array(order, m, c, &legendre[t][l] );
          cout << legendre[t][l] << endl;
#endif
        }
      }
    }
  }
}

When I compile using GSL 1.9 I use the deprecated function gsl_sf_legendre_sphPlm_array, while when I compute with GSL 2.5 I use the new function gsl_sf_legendre_array which explicitely requires a keyword for the normalization (GSL_SF_LEGENDRE_SPHARM) and doesn't ask for the parameter m. The old version gives me consistent results, while the new one returns me a segmentation fault after 25 t loops. Can anybody of you help me correct the code and explain me what I'm doing wrong?

Answer 1

I tried to compile with the debugging symbols (the -g flag in the command below which assumes your program is called "main.cpp")

usr@cmptr $ g++ -g main.cpp -lgsl -lgslcblas -lm -o main

And ran the program using a debugger, e.g. gdb (GNU debugger):

usr@cmptr $ gdb main
(gdb) run main
Starting program: /home/usr/Desktop/main 
0, 0, 0.282095
1, 0, 0.282095
2, 0, 0.282095
3, 0, 0.282095
4, 0, 0.282095
5, 0, 0.282095
6, 0, 0.282095
7, 0, 0.282095
8, 0, 0.282095
9, 0, 0.282095
10, 0, 0.282095
11, 0, 0.282095
12, 0, 0.282095
13, 0, 0.282095
14, 0, 0.282095
15, 0, 0.282095
16, 0, 0.282095
17, 0, 0.282095
18, 0, 0.282095
19, 0, 0.282095
20, 0, 0.282095
21, 0, 0.282095
22, 0, 0.282095
23, 0, 0.282095
24, 0, 0.282095

Program received signal SIGSEGV, Segmentation fault.
0x0000555555554ded in main () at main.cpp:26
26            cout << t << ", " << l << ", " << legendre[t][l] << endl;
(gdb) quit

The error is in line 26 where you write to screen. The reason is probably that you are trying to access the array out of bounds. Maybe having a look at how to allocate the memory "properly" in the manual for gsl Legendre polynomials. I am specifically thinking of the functions gsl_sf_legendre_array_n and gsl_sf_legendre_array_index.

Note: I have not used this part of GSL myself, but I usually find that when using GSL it is useful to use temporary variables and arrays that you pack/unpack before/after calling functions. Not pretty but less error prone, and since you are using the "plusspluss" version of C you can always wrap the implementation in a class. Hopefully it is helpful.

---

Edit:

I tried increasing the size of the legendre array and the program completes when the size is set to:

double legendre[ntheta + 10][order + 4];

Maybe someone who knows how the functions work can answer why...