bugs

math library from gnu

------------------------------------------------------------------------
BUG-ID:   21819   
STATUS:   Open/Confirmed  
CATEGORY: Runtime error
SUMMARY:  gsl_sf_bessel_J_CF1() crash for large arguments

Reply-To: Koichi Takahashi <[e-mail]>
From: Koichi Takahashi <[e-mail]>
To: [e-mail]
CC: Brian Gough <[e-mail]>, Jonathan Taylor <[e-mail]>
Subject: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 02 Dec 2007 00:11:02 -0800

Hi,

Unfortunately, I found some additional sets of numbers that still
crash gsl_sf_bessel_J_CF1() even with the cvs version.   Symptom
is exactly the same as what I reported before.  I tested on
x86_64.

main()
{
//    this used to crash, but now fixed in the current cvs.
//    double a = gsl_sf_bessel_jl( 30,  3875.6138424501978 );

//    at least the following three sets of values still crashes.
     double a = gsl_sf_bessel_jl( 49, 9912.6308956132361 );
//    double a = gsl_sf_bessel_jl( 49, 9950.3478935215589 );
//    double a = gsl_sf_bessel_jl( 52, 9930.5181281798232 );

     printf("%g\n",a);
}

Let me know if there is anything I could do to help you fixing
this issue.

thanks,
Koichi

From: Jonny Taylor <[e-mail]>
To: Koichi Takahashi <[e-mail]>
Cc: [e-mail]
Subject: [Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 2 Dec 2007 10:49:33 +0000

While the symptom is the same, the cause is different. For those  
numbers it seems that 10,000 iterations is simply not enough.  
Interestingly, In fact, in all three cases it requires less than 50  
additional iterations to converge!?

The naive fix is simply to increase the maximum permitted number of  
iterations, but a more sophisticated fix would probably need to  
justify a maximum number of iterations, and propose an alternative  
method of generating the result in the cases where the number of  
iterations is prohibitive...

Jonny
From: Koichi Takahashi <[e-mail]>
To: Jonny Taylor <[e-mail]>
Cc: [e-mail]
Subject: [Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 02 Dec 2007 06:06:11 -0800

At this range, x is still not large enough to use the asymptotic form?
In gsl_sf_bessel_jl_e,

   else if(x > 1000.0 && x > 100.0*l*l)
   {
     //asymptotic
   }
   else
   {
     //CF1
   }

so, for example, 100 * 50 * 50 = 250,000.
For l=50, the iteration starts to exceed 10,000 around x=9900.
If we want to stick to 10,000 max iteration, we have to switch to
the asymptotic version with something like x > 3*l*l.  Maybe this is
too small?   If so I'd consider increasing the max iteration.
Maybe we should do both.  Or there can be a better method.

koichi

------------------------------------------------------------------------
BUG-ID:   21826   
STATUS:   Open            
CATEGORY: Build
SUMMARY:  libtool problem on hp-ux

From: "Benoit, Gerard" <[e-mail]>
To: [e-mail]
Subject: [Bug-gsl] gsl-1.9 : error in compilation scripts
Date: Fri, 13 Jul 2007 12:28:35 +0200

Hello,

I wish to report a bug in the compilation scripts:
GSL version 1.9  downloaded from ftp.gnu.org
computer HP 9000/785 under HP-UX  B.11.00 A 
config.guess -> hppa2.0w-hp-hpux11.00
compiler  : cc -> /opt/ansic/bin/cc   HP C  HP92453-01 A.11.01.20
compilation options : -Ae
linker  : ld -> /bin/ld

fatal error during the build of the shared libraries( first for cblas/libgslcblas.sl )
command executed : $archive_cmds 
gsl-1.9/libtool ( line 214 ) : archive_cmds="\$CC -b \${wl}+h \${wl}\$soname \${wl}+b \${wl}\$install_libdir -o \$lib \$libobjs \$deplibs \$compiler_flags"
coming from gsl-1.9/aclocal.m4 ( line 6541 ) : 
   _LT_AC_TAGVAR(archive_cmds, $1)='$CC -b ${wl}+h ${wl}$soname ${wl}+b ${wl}$install_libdir -o $lib $libobjs $deplibs $compiler_flags'

the compiler cc don't recognize the "-b" option ( it is a linker option ) and throw it : the linker looks for a main module and don't find it
with the correct form ${wl}-b to pass the option to the linker, cc add a start-up module ( like /usr/ccs/lib/crt0.o or /opt/langtools/lib/crt0.o ) which is useless for a shared library and has not be compiled with a PIC option ( +z or +Z ) and the linker stops.

In the preceding version 1.6 ( I have never downloaded version 1.7 and 1.8 ) the command was :
gsl-1.6/libtool ( line 187 ) : archive_cmds="\$LD -b +h \$soname +b \$install_libdir -o \$lib \$libobjs \$deplibs \$linker_flags"
from
gsl-1.6/aclocal.m4 ( line 2441 ) :    *) archive_cmds='$LD -b +h $soname +b $install_libdir -o $lib $libobjs $deplibs $linker_flags' ;;

This command works fine

so, look for a way to return to the LD command in gsl-1.9/aclocal.m4 ( line 6541 ) :

Best regards

Gérard BENOIT

Ce courriel (incluant ses éventuelles pièces jointes) peut contenir des informations confidentielles et/ou protégées ou dont la diffusion est restreinte. Si vous avez reçu ce courriel par erreur, vous ne devez ni le copier, ni l'utiliser, ni en divulguer le contenu à quiconque. Merci d’en avertir immédiatement l’expéditeur et d’effacer ce courriel de votre système. Astrium décline toute responsabilité en cas de corruption par virus, d’altération ou de falsification de ce courriel lors de sa transmission par voie électronique.

This email (including any attachments) may contain confidential and/or privileged information or information otherwise protected from disclosure. If you are not the intended recipient, please notify the sender immediately, do not copy this message or any attachments, do not use it for any purpose or disclose its content to any person, but delete this message and any attachments from your system. Astrium disclaims any and all liability if this email transmission was virus corrupted, altered or falsified.
----------------------------------------------------------------
Astrium SAS (393 341 516 RCS Paris) - Siège social: 6 rue Laurent Pichat, 75016 Paris, France
_______________________________________________
Bug-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/bug-gsl

------------------------------------------------------------------------
BUG-ID:   21828   
STATUS:   Open/Confirmed  
CATEGORY: Performance
SUMMARY:  suboptimal performance of gsl_fdfsolver_lmsder

From: "Alexander Usov" <[e-mail]>
To: [e-mail]
Subject: [Help-gsl] Strange performance of gsl_fdfsolver_lmsder
Date: Wed, 24 Oct 2007 20:45:01 +0200

Hi all,

I am currently working on the problem involving source extraction from
astronomical images, which essentially boils down to fitting a number of
2d gaussians to the image.

One of the traditionally used fitters in this field is a Levenberg-Marquardt,
which gsl_fdfsolver_lmsder is and implementation of.

At some moment I have notices that for the bigger images (about 550
pixels, 20-30 parameters) gsl's lmsder algorithm spends a large fraction
of the run-time (about 50%) doing household transform.

While looking around for are different minimization algorithms I have made
a surprising finding that original netlib/minpack/lmder is almost twice faster
that that of gsl.

Could anyone explain such a big difference in performace?

-- 
Best regards,
  Alexander.

_______________________________________________
Help-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/help-gsl

Reply-To: [e-mail]
From: Brian Gough <[e-mail]>
To: "Alexander Usov" <[e-mail]>
Cc: [e-mail]
Subject: Re: [Help-gsl] Strange performance of gsl_fdfsolver_lmsder
Date: Thu, 25 Oct 2007 21:57:08 +0100

At Wed, 24 Oct 2007 20:45:01 +0200,
Alexander Usov wrote:
> At some moment I have notices that for the bigger images (about 550
> pixels, 20-30 parameters) gsl's lmsder algorithm spends a large fraction
> of the run-time (about 50%) doing household transform.
> 
> While looking around for are different minimization algorithms I have made
> a surprising finding that original netlib/minpack/lmder is almost twice faster
> that that of gsl.
> 
> Could anyone explain such a big difference in performace?

I have a vague memory that there was some quantity (Jacobian?) that
MINPACK only computes fully at the end, but in GSL it is accessible to
the user at each step so I felt I had to update it on each iteration
in the absence of some alternate scheme.  Sorry this is not a great
answer but I am not able to look at it in detail now.

-- 
Brian Gough

_______________________________________________
Help-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/help-gsl

------------------------------------------------------------------------
BUG-ID:   21830   
STATUS:   Open/Confirmed  
CATEGORY: None
SUMMARY:  provide function to retrieve size of chebyshev data

Reply-To: [e-mail]
From: Brian Gough <[e-mail]>
To: [e-mail]
Subject: Re: [Help-gsl] Writing Chebychev coefficients to disk
Date: Mon, 12 Mar 2007 21:40:21 +0000

At Mon, 12 Mar 2007 10:28:10 -0600,
Patrick Alken wrote:
> 
> It seems to me the simplest way to do this is just:
> 
> fwrite(gsl_cheb_ptr, sizeof(gsl_cheb_series), 1, fp);
> fwrite(gsl_cheb_ptr->c, sizeof(double), n, fp);
> 
> and later:
> 
> fread(gsl_cheb_ptr, sizeof(gsl_cheb_series), 1, fp);
> fread(gsl_cheb_ptr->c, sizeof(double), n, fp);
> 

Yes, the library should help in getting the correct value of "n"
(which is c->order + 1 in this case) to avoid subtle off-by-one
errors.

It's certainly more useful to expose the memory than provide the read
and write functions directly though.

-- 
Brian Gough

Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/

_______________________________________________
Help-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/help-gsl

None

------------------------------------------------------------------------
BUG-ID:   21831   
STATUS:   Open            
CATEGORY: Accuracy problem
SUMMARY:  Levý random number generator for alpha < 1

From: [e-mail]
To: [e-mail]
Subject: [Bug-gsl] Levý random number generator
Date: Mon, 26 Mar 2007 19:48:01 -0300

The Levý skew random number generator (gsl_ran_levy_skew) does not  
procuce a Levý random number when beta=0 (symmetric case), and the  
gsl_ran_levy function does not work as stated in the docs. I made some  
histograms from 10^6 samples to check the accuracy of the algorithms,  
by comparison agaisnt the numerical integration of the equation of  
Levý's PDF. For the gsl_ran_levy function there is a good precison for  
alpha [1,2], for alpha (0.3,1) you must sum a series of random numbers  
to get the same precision (tipicaly 100 or more gsl_ran_levy numbers).  
For alpha<=0.3 the algorithm does not work properly, even worse, the  
error increases as you add more random numbers. This contradicts the  
manual that says "the algoritm only works for alpha (0,2]". The  
function gsl_ran_levy_skew does not produce levy random numbers when  
beta=0, instead the pdf of the random numbers is a linear (?!?!) one.

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

_______________________________________________
Bug-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/bug-gsl

From: [e-mail]
To: Brian Gough <[e-mail]>
Cc: 
Subject: Re: [Bug-gsl] Lev� random number generator
Date: Tue, 27 Mar 2007 09:35:15 -0300

Thanks for your quick answer, and sorry about my poor english, it is  
not my natural language.

The code below generates 10^6 random numbers, and makes a normalized  
histogram wich is compared to the levy pdf. To get the levy pdf, it  
numericaly integrates the characteristic function for levy process  
(the function f in the code). The n parameter just adds a series of  
levy numbers to get better precision. The code saves 2 files:  
lhist-$alpha-$n (the normalized histogram) and lpdf-$alpha (the pdf  
for the levy process). It also prints to the stdout the absolute error  
(square of the difference) between the histogram and the pdf.

The function levy skew shows problems for alpha<1.

With this code, you can also check the problems related to the  
gsl_ran_levy function (just change gsl_ran_levy_skew by gsl_ran_levy,  
cutting the last paramenter).

I am using the pre-compiled gsl that comes with debian etch (gsl  
version 1.8.2).

If you are interested, I also encoded a routine to generate levy skew  
random numbers, it is not fully tested, but it works for beta=0 and  
alpha<1 (it suffers from the same precision problem as gsl_ran_levy  
function for small alpha)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_integration.h>

double f (double x, void *params)
{
	double alpha = *(double *) params;
	double f = exp(-pow(x,alpha))/M_PI;
	return f;
}

double *levy_pdf(double alpha,int hist_size,double a,double b)
{
	double abserr,*lpdf,dx;
	gsl_function F;
	int i;
	gsl_integration_workspace *w=gsl_integration_workspace_alloc(1000);
	gsl_integration_workspace *cw=gsl_integration_workspace_alloc(1000);
	gsl_integration_qawo_table  
*wf=gsl_integration_qawo_table_alloc(0.0,1.0,GSL_INTEG_COSINE,200);
	lpdf=(double*)calloc(hist_size,sizeof(double));
	F.function=&f;
	F.params=&alpha;
	dx=(double)(b-a)/(double)hist_size;
	for (i=0;i<hist_size;i++)
	{
		gsl_integration_qawo_table_set(wf,i*dx+a,1.0,GSL_INTEG_COSINE);
		gsl_integration_qawf (&F,0.0,1e-10,1000,w,cw,wf,&lpdf[i],&abserr);
	}
	gsl_integration_qawo_table_free(wf);
	gsl_integration_workspace_free(w);
	gsl_integration_workspace_free(cw);
	return (lpdf);
}

int main (int argc,char *argv[])
{
	double *l,*lpdf,a=-20,b=20,alpha,dx,n,errabs=0.0;
	unsigned long int rnd_seed;
	int i,j,rand_numbers=1e6,hist_size=400;
	gsl_histogram *h;
	gsl_rng *r;
	struct timeval *tv;
	struct timezone *tz;
	char filename[50];
	FILE *f1,*f2;
	if(argc!=3)
	{
		printf("\nThe program must be called with 2 parameters: alpha and n\n \n");
		exit(1);
	}
	dx=(double)(b-a)/(double)hist_size;
	h=gsl_histogram_alloc(hist_size);
	alpha=atof(argv[1]);
	n=atof(argv[2]);
	strcpy(filename,"lhist-");
	strcat(filename,argv[1]);
	strcat(filename,"-");
	strcat(filename,argv[2]);
	f1=fopen(filename,"w+");