gsl_sf_hyperg_2F1: gsl: gamma.c:769: ERROR: error

[Bruce Rout]

I am using Ubuntu 20.04 on HP Z600 with 24G of RAM.

I installed gsl library through sudo apt-get install libgsl-dev yesterday.
I am invoking the function gsl_sf_hyperg_2F many times in a loop. I receive
error as follows:

gsl: gamma.c:769: ERROR: error
Default GSL error handler invoked.
Aborted (core dumped)

The function was called with the following parameters:

Hypergeom parameters a: 1.500000 b: 1.000000 c: 0.500000 d: 0.999013
factor1: 0.001395

The factor1 parameter at the end is not part of the call but is used in
determining a characteristic function that uses the Hypergeometric
function. The d parameter is the same as z or x in most documentation of
the function.

I have tried few loops in calling the function and the function works
without a problem. If I increase the number of calls, it crashes. I will
include the code for your examination. Near the top of the code, in the
Declaration segment, if npoints is defined to be 200, the code is robust.
If npoints is set to 5000, it crashes. I have read the documentation of
previous problems with this functiion and I hope this may add some light.

The compile is done through a script called compile and is as follows

gcc -Wall -I/usr/include/gsl -c $1.c
gcc -L/usr/local/lib $1.o -lgsl -lgslcblas -lm -o $1
echo "compile finished"

This is the code:

//**********************
// hyper.c                    **
// hypergeometric function    **
//**********************

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl_sf_hyperg.h>
//**********************************************
//       Declaration Segment
//**********************************************


    // Number of rows and columns in a matrix
    // you can interactively ask for these two numbers

#define ncoef 81  // must be an exact square
#define npoints 5000 // must follow the formula of nphi+ntheta where nphi =
sqrt(npoints/2) and ntheta = nphi * 2
// npoints must be greater than ncoef

// Matrix holders are created

double
XT[ncoef][npoints],X[npoints][ncoef],XTX[ncoef][ncoef],Y[npoints],YT[npoints];
double
XTXINV[ncoef][ncoef],XTXINVXT[ncoef][npoints],I[ncoef][ncoef],A[ncoef];

int ntheta = 2 * sqrt(npoints/2);
int nphi   = sqrt(npoints/2);

// define pi

const double pi = 4.0 * atan(1.0);

//***********************************************
//       Functions
//***********************************************


int main(){

int i,j,m,l;

    double Phi,Theta,a,b,c,d,factor1,r=1.1;
double dtheta,theta,dphi,phi;
    int nm=sqrt(ncoef);
    int nl=sqrt(ncoef);
int ipoint,jcol;

dtheta = 2*pi/(ntheta);  //  repeat theta= 0  and theta = 2*pi before
initial conditions
dphi=pi/(nphi); // stay away from the poles after initial conditions




//  set a number of graph data for different values of phi, then graph
hypergeom vs phi

theta=-dtheta/2.0;  // cannot repeat theta= 0  and theta = 2*pi during
initial conditions
ipoint=0;
for (i=0;i<ntheta;i++){
     theta=theta+dtheta;
     phi=-dphi/2.0; // stay away from the poles after initial conditions
     for (j=0;j<nphi;j++){
          phi=phi+dphi;
          jcol=0;
          for (m = 0;m<nm;m++){
               for (l = 0; l < nl;l ++){

                    Theta = cos(m*theta);

                    factor1 = sqrt(2.0)*pow(.5-.5*cos(2.*phi),m*.5);
                    double surd = 4*(l*l*r*r)+1;
                    a=m/2. + sqrt(surd)/4.+1./4.;
                    b=m/2. - sqrt(surd)/4.+1./4.;
                    c=1.0/2.0;
                    d=cos(2*phi)/2. + 1./2.;
                    printf("Hypergeom parameters a: %f b: %f c: %f d: %f
factor1: %f\n",a,b,c,d,factor1);
                    printf("Calling hypergeom i = %d j= %d m = %d l = %d
\n",i,j,m,l);
                    double H = gsl_sf_hyperg_2F1(a, b, c, d);
                    printf("Back from hypergeom\n");
                    Phi = factor1*H;

                    printf("hypergeom(a,b,c,d) = %f,\tPhi = %f\n",H,Phi);
                    if (Phi == 0.0){
                         printf("Phi is zero\n");
                         printf("Press any key\n");
                        getchar();
                    }
                    X[ipoint][jcol]=Theta*Phi;
                    XT[jcol][ipoint]=X[ipoint][jcol];
                   jcol++;
               }
          }
          ipoint++;
     }
}
}

Here is the output:

.
.
.

Hypergeom parameters a: 4.608108 b: -3.108108 c: 0.500000 d: 0.999013
factor1: 0.044422
Calling hypergeom i = 0 j= 0 m = 1 l = 7
Back from hypergeom
hypergeom(a,b,c,d) = 3190547883224283.000000, Phi = 141728991368209.000000
written
Hypergeom parameters a: 5.157097 b: -3.657097 c: 0.500000 d: 0.999013
factor1: 0.044422
Calling hypergeom i = 0 j= 0 m = 1 l = 8
Back from hypergeom
hypergeom(a,b,c,d) = 214.861179, Phi = 9.544460
written
Hypergeom parameters a: 1.500000 b: 1.000000 c: 0.500000 d: 0.999013
factor1: 0.001395
Calling hypergeom i = 0 j= 0 m = 2 l = 0
gsl: gamma.c:769: ERROR: error
Default GSL error handler invoked.
Aborted (core dumped)


I cannot find gamma.c line # 769 in the library and I don't have any error
trapping written into the code.

Thank you.

Yours,

Bruce