[bug #65728] Add sparse functionalities

[Patrick Alken]

URL:
  <https://savannah.gnu.org/bugs/?65728>

 Summary: Add sparse functionalities
   Group: GNU Scientific Library
   Submitter: psa
   Submitted: Sun 12 May 2024 01:09:04 PM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any


___

Follow-up Comments:


---
Date: Sun 12 May 2024 01:09:04 PM UTC By: Patrick Alken 
from Brijesh Upadhaya 

I don't know if this is necessary, but the scaling of a complex (sparse)
matrix with a real vector, complex sparse matrix-vector multiplication, and
unpacking(in place) a complex vector could be added. I had to somehow do the
above-mentioned stuff while solving a time-harmonic FE problem.

There are alloc and calloc for the blocks and vectors. Is there any particular
reason for not having anything to resize a block or a vector?

In gsl/spmatrix/getset_source.c, _set() replaces the matrix elements, could it
be possible to assemble the matrix as

*(BASE*) ptr = *(BASE*) + x;

and copy the entire code to a new name.







___

Reply to this item at:

  <https://savannah.gnu.org/bugs/?65728>

___
Message sent via Savannah
https://savannah.gnu.org/

[bug #63679] C99 compatibility fix for the configure script

[Patrick Alken]

Update of bug #63679 (group gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #2:

Patch applied in commit ac250bf26d80f723319acde73c4f42475d81fc98


___

Reply to this item at:

  

___
Message sent via Savannah
https://savannah.gnu.org/

Re: Typo in Mersenne Twister PRNG's leading comment

[Patrick Alken]

Thank you David, I fixed it on the git

On 6/26/23 10:47, David Grimm wrote:

Hello,

as of commit b6c4b6e077679a4ce9b4a6e0889d0e59672acc3a of the git 
repository git.savannah.gnu.org/gsl.git,
in the file rng/mt.c, lines 43-44, the title of the paper describing 
the Mersenne Twister algorithm is cited as "Mersenne Twister: A 
623-dimensionally equidistributerd uniform pseudorandom number 
generator", adding an additional r in "equidistributed".


yours,
David Grimm

Re: -lcblas needed with -latlas

[Patrick Alken]

If you're using the ATLAS library, you need both -latlas and -lcblas 
together. If you are using the GSL cblas, you only need -lgslcblas


On 4/25/22 15:15, Crawford, Christopher wrote:

https://www.gnu.org/software/gsl/doc/html/usage.html#an-example-program
I was looking through this documentation, and couldn't test it, but it
seems like the following command

$ gcc example.o -lgsl -lcblas -latlas -lm

shouldn't have -lcblas.  I could be wrong, though.
--thanks, Chris

Christopher Crawford, Professor, Director of Graduate Studies
Department of Physics and Astronomy, University of Kentucky
859-257-2504 (o)  859-421-9612 (c)   skype: chris2crawford
schedule: www.pa.uky.edu/~crawford   uky.zoom.us/my/c.crawford

Re: Segmentation fault in gsl_integration_fixed (GSL Integration)

[Patrick Alken]

I don't think this is a bug. For Hermite integration, b must be > 0 (see 
documentation). For your input parameters, the weight function is 1, so 
you could just use Legendre integration


On 4/19/22 18:57, Zhoulai Fu wrote:

A segmentation fault is raised in gsl_integration_fixed function in GSL
Integration. The program that triggers the segmentation fault is attached.
Below is debugging info from gdb.




*Program received signal SIGSEGV, Segmentation fault.0x7ff256302bd9 in
gsl_integration_fixed (func=func@entry=0x7ffec35ea0e0,
result=result@entry=0x7ffec35ea0d8, w=w@entry=0x0)at fixed.c:134134
  const size_t n = w->n*

BR,
Zhoulai

Re: Potential security bug: Buffer overflow in gsl_stats_quantile_from_sorted_data (of library Statistics)

[Patrick Alken]

Hi Zhoulai,

  I appreciate you following up on this, and for your patience. I have 
just committed a fix to the git repository which should address the 
issue. I added a check on the input parameter 'f' which returns an error 
if it is out of range.


Best regards,

Patrick

On 4/16/22 08:46, Zhoulai Fu wrote:

I just found out that the buffer overflow issue I reported years ago (see
below) remains in the recent version of GSL (just tested on the newest one
on GitHub). Is this issue something we plan to fix? I am asking since
now computing quantiles with GSL, or anything depending on it,  seems not
secure.

BR,
Zhoulai

On Thu, Dec 3, 2020 at 11:12 PM Zhoulai Fu@Gmail 
wrote:


Running the following code (also attached as a file) triggers a
segmentation error.













*#include #include #include
int main(void){  double upperq;  double data[5] =
{17.2, 18.1, 16.5, 18.3, 12.6};  gsl_sort (data, 1, 5);  upperq =
gsl_stats_quantile_from_sorted_data (data, 1, 5, 675);  return 0;}// gcc
statsort_bug.c -lgsl -lgslcblas; ./a.out*

The error points to statistics/quantiles_source.c:41:


*  result = (1 - delta) * sorted_data[lhs * stride] + delta *
sorted_data[(lhs + 1) * stride] ;*
The segmentation error is due to a stack buffer overflow (where
lhs=2700, strid=1 as shown in GDB). The bug could be exploited for
security attack, knowing that it occurs when the quantile "f" is
beyond the expected [0,1] range (f=675 in this case).

BR,
Zhoulai

[bug #59624] Buffer overflow in gsl_stats_quantile_from_sorted_data

[Patrick Alken]

Update of bug #59624 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #2:

fixed in commit 989a193268b963aa1047814f7f1402084fb7d859


___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: Fwd: bug #59781

[Patrick Alken]

Michael, thank you again for tracking this down, your code and 
explanation was very useful. I have added your fix to the git repository


On 3/24/22 14:26, Michael Dunlap wrote:

Hello, please forgive the briefness of my previous email, I was rather
excited.

Allow me to introduce myself.

My name is Michael Dunlap. I graduated from the University of Arizona in
2018 with a
bachelor's degree. I majored in math and minored in computer science.
Having made use of
the GNU Science Library, I thought it would be worth checking to see if
there was a way I
could make a contribution to the project. The website says that anyone
interested in
contributing should start with working on the issues mentioned in the bug
tracker. Being
acquainted with the QR algorithm I decided to work on bug #59781.

Thank you very much to Patrick Alken for including well written commentary
in the code.
Without the guidance to Matrix Computations and DLAHQR, working with the
code would have
been substantially more difficult.

Part of algorithm 7.5.2 from Matrix Computations involves looking for
subdiagonal entries
that are sufficiently close to zero. If any are found we can break the
problem up into two
smaller subproblems ("deflation").

In the book Matrix Computations the criterion by which a subdiagonal entry
is deemed small
enough is:

|h_i,i-1| <= tol * ( |h_i,i| + |h_i-1,i-1| )

where |.| refers to absolute value, and tol is a suitably small number. In
the code the
value GSL_DBL_EPSILON is used for tol.

Here is a snippet of output generated by adding print statements to the
code that shows the
problem quite clearly:

[image: fig1.png]
(fig 1)

In the above we are checking the subdiagonal for sufficiently small
entries. The program
fails to find any and proceeds to perform a francis qrstep using Wilkinson
shift.
However it is quite clear that the subdiagonal entry -2.47033e-323 is small
enough
according to the aforementioned criterion. The program should have set it
to zero.

The precise nature of the bug can be stated succinctly as:
In the function francis_search_subdiag_small_elements the variables del and
dpel are
backwards.

Here is the original code:

static inline size_t
francis_search_subdiag_small_elements(gsl_matrix * A)
{
   const size_t N = A->size1;
   size_t i;
   double dpel = gsl_matrix_get(A, N - 2, N - 2);

   for (i = N - 1; i > 0; --i)
 {
   double sel  = gsl_matrix_get(A, i, i - 1);
   double del  = gsl_matrix_get(A, i, i);

   if ((sel == 0.0) ||
   (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel
 {
   gsl_matrix_set(A, i, i - 1, 0.0);
   return (i);
 }

   dpel = del;
 }

   return (0);
} /* francis_search_subdiag_small_elements() */

At first del is set to h_N-1,N-1 and dpel is set to h_N-2,N-2. In the
second pass of the
loop when the comparison is made del = h_N-2,N-2 and dpel is set to
h_N-1,N-1.
This does not match algorithm 7.5.2 from Matrix Computations.

Going back to the snippet we see that the program is comparing
-2.47033e-323 to the two
zero entries on the diagonal. Hence "deflation" of the problem fails to
occur. The
algorithm should be comparing -2.47033e-323 to 1.21667e-16 and 0. Instead
it is comparing
-2.47033e-323 to the two zero valued diagonal entries. "Circling" the entry
for sel with
<,> and the entries for del and dpel with [,] we have the following:


[image: fig2.png]
(fig 2)

Now, the determination of "close enough to zero" is not necessarily unique.
So perhaps it could be argued that this
"bug" is not really a bug but an alternative way of making such
determination.

As has already been mentioned, gsl_eigen_nonsymm works on the transpose of
the matrix in
the bug submission. Additionally, balancing the matrix prior to handing it
off to  gsl_eigen_nonsymm also results in
getting the correct eigenvalues.

Sincerely,
Michael Dunlap

-- Forwarded message -
From: Michael Dunlap 
Date: Fri, Mar 18, 2022 at 12:23 AM
Subject: bug #59781
To: 


Hello!
I believe that I have found the bug and fixed it.

The bug was pertaining to the variable dpel in the  function:
francis_search_subdiag_small_elements.

I have attached a corrected version with comments to show where the issue
was.

I have tested this version of francis.c against the matrix submitted by
Dmitry Cheshkov and also checked it against test.c, both were successful.

Thanks,
Michael Dunlap

[bug #59781] Eigenvalue failure

[Patrick Alken]

Update of bug #59781 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #2:

fixed in commit cf8f4650b8aef8ac8172656939d25080dcbb33ef


___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: Fwd: bug #59781

[Patrick Alken]

Dear Michael,

  Thank you for your great work on this bug. I apologize for the slow 
response. It seems you have identified the cause of the problem. I am 
going to look at your patch and do some additional testing to see if 
everything looks ok. I had written a more comprehensive test program for 
the eigenvalue solver some time ago so I will run that program to make 
sure everything is ok.


Please give me a few days to look into everything and make sure the 
testing looks ok.


Best regards,

Patrick

On 3/24/22 14:26, Michael Dunlap wrote:

Hello, please forgive the briefness of my previous email, I was rather
excited.

Allow me to introduce myself.

My name is Michael Dunlap. I graduated from the University of Arizona in
2018 with a
bachelor's degree. I majored in math and minored in computer science.
Having made use of
the GNU Science Library, I thought it would be worth checking to see if
there was a way I
could make a contribution to the project. The website says that anyone
interested in
contributing should start with working on the issues mentioned in the bug
tracker. Being
acquainted with the QR algorithm I decided to work on bug #59781.

Thank you very much to Patrick Alken for including well written commentary
in the code.
Without the guidance to Matrix Computations and DLAHQR, working with the
code would have
been substantially more difficult.

Part of algorithm 7.5.2 from Matrix Computations involves looking for
subdiagonal entries
that are sufficiently close to zero. If any are found we can break the
problem up into two
smaller subproblems ("deflation").

In the book Matrix Computations the criterion by which a subdiagonal entry
is deemed small
enough is:

|h_i,i-1| <= tol * ( |h_i,i| + |h_i-1,i-1| )

where |.| refers to absolute value, and tol is a suitably small number. In
the code the
value GSL_DBL_EPSILON is used for tol.

Here is a snippet of output generated by adding print statements to the
code that shows the
problem quite clearly:

[image: fig1.png]
(fig 1)

In the above we are checking the subdiagonal for sufficiently small
entries. The program
fails to find any and proceeds to perform a francis qrstep using Wilkinson
shift.
However it is quite clear that the subdiagonal entry -2.47033e-323 is small
enough
according to the aforementioned criterion. The program should have set it
to zero.

The precise nature of the bug can be stated succinctly as:
In the function francis_search_subdiag_small_elements the variables del and
dpel are
backwards.

Here is the original code:

static inline size_t
francis_search_subdiag_small_elements(gsl_matrix * A)
{
   const size_t N = A->size1;
   size_t i;
   double dpel = gsl_matrix_get(A, N - 2, N - 2);

   for (i = N - 1; i > 0; --i)
 {
   double sel  = gsl_matrix_get(A, i, i - 1);
   double del  = gsl_matrix_get(A, i, i);

   if ((sel == 0.0) ||
   (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel
 {
   gsl_matrix_set(A, i, i - 1, 0.0);
   return (i);
 }

   dpel = del;
 }

   return (0);
} /* francis_search_subdiag_small_elements() */

At first del is set to h_N-1,N-1 and dpel is set to h_N-2,N-2. In the
second pass of the
loop when the comparison is made del = h_N-2,N-2 and dpel is set to
h_N-1,N-1.
This does not match algorithm 7.5.2 from Matrix Computations.

Going back to the snippet we see that the program is comparing
-2.47033e-323 to the two
zero entries on the diagonal. Hence "deflation" of the problem fails to
occur. The
algorithm should be comparing -2.47033e-323 to 1.21667e-16 and 0. Instead
it is comparing
-2.47033e-323 to the two zero valued diagonal entries. "Circling" the entry
for sel with
<,> and the entries for del and dpel with [,] we have the following:


[image: fig2.png]
(fig 2)

Now, the determination of "close enough to zero" is not necessarily unique.
So perhaps it could be argued that this
"bug" is not really a bug but an alternative way of making such
determination.

As has already been mentioned, gsl_eigen_nonsymm works on the transpose of
the matrix in
the bug submission. Additionally, balancing the matrix prior to handing it
off to  gsl_eigen_nonsymm also results in
getting the correct eigenvalues.

Sincerely,
Michael Dunlap

-- Forwarded message -
From: Michael Dunlap 
Date: Fri, Mar 18, 2022 at 12:23 AM
Subject: bug #59781
To: 


Hello!
I believe that I have found the bug and fixed it.

The bug was pertaining to the variable dpel in the  function:
francis_search_subdiag_small_elements.

I have attached a corrected version with comments to show where the issue
was.

I have tested this version of francis.c against the matrix submitted by
Dmitry Cheshkov and also checked it against test.c, both were successful.

Thanks,
Michael Dunlap

[bug #59781] Eigenvalue failure

[Patrick Alken]

Follow-up Comment #1, bug #59781 (project gsl):

emails from M. Dunlap mdunlap1 =at= email =dot= arizona =dot= edu

Part of algorithm 7.5.2 from Matrix Computations involves looking for
subdiagonal entries
that are sufficiently close to zero. If any are found we can break the
problem up into two
smaller subproblems ("deflation").

In the book Matrix Computations the criterion by which a subdiagonal entry
is deemed small
enough is:

|h_i,i-1| <= tol * ( |h_i,i| + |h_i-1,i-1| )

where |.| refers to absolute value, and tol is a suitably small number. In
the code the
value GSL_DBL_EPSILON is used for tol.

Here is a snippet of output generated by adding print statements to the
code that shows the
problem quite clearly:

[image: fig1.png]
(fig 1)

In the above we are checking the subdiagonal for sufficiently small
entries. The program
fails to find any and proceeds to perform a francis qrstep using Wilkinson
shift.
However it is quite clear that the subdiagonal entry -2.47033e-323 is small
enough
according to the aforementioned criterion. The program should have set it
to zero.

The precise nature of the bug can be stated succinctly as:
In the function francis_search_subdiag_small_elements the variables del and
dpel are
backwards.

Here is the original code:

static inline size_t
francis_search_subdiag_small_elements(gsl_matrix * A)
{
  const size_t N = A->size1;
  size_t i;
  double dpel = gsl_matrix_get(A, N - 2, N - 2);

  for (i = N - 1; i > 0; --i)
{
  double sel  = gsl_matrix_get(A, i, i - 1);
  double del  = gsl_matrix_get(A, i, i);

  if ((sel == 0.0) ||
  (fabs(sel) < GSL_DBL_EPSILON * (fabs(del) + fabs(dpel
{
  gsl_matrix_set(A, i, i - 1, 0.0);
  return (i);
}

  dpel = del;
}

  return (0);
} /* francis_search_subdiag_small_elements() */

At first del is set to h_N-1,N-1 and dpel is set to h_N-2,N-2. In the
second pass of the
loop when the comparison is made del = h_N-2,N-2 and dpel is set to
h_N-1,N-1.
This does not match algorithm 7.5.2 from Matrix Computations.

Going back to the snippet we see that the program is comparing
-2.47033e-323 to the two
zero entries on the diagonal. Hence "deflation" of the problem fails to
occur. The
algorithm should be comparing -2.47033e-323 to 1.21667e-16 and 0. Instead
it is comparing
-2.47033e-323 to the two zero valued diagonal entries. "Circling" the entry
for sel with
<,> and the entries for del and dpel with [,] we have the following:


[image: fig2.png]
(fig 2)

Now, the determination of "close enough to zero" is not necessarily unique.
So perhaps it could be argued that this
"bug" is not really a bug but an alternative way of making such
determination.

As has already been mentioned, gsl_eigen_nonsymm works on the transpose of
the matrix in
the bug submission. Additionally, balancing the matrix prior to handing it
off to  gsl_eigen_nonsymm also results in
getting the correct eigenvalues.

Sincerely,
Michael Dunlap

===
original email
===


Hello!
I believe that I have found the bug and fixed it.

The bug was pertaining to the variable dpel in the  function:
francis_search_subdiag_small_elements.

I have attached a corrected version with comments to show where the issue
was.

I have tested this version of francis.c against the matrix submitted by
Dmitry Cheshkov and also checked it against test.c, both were successful.

Thanks,
Michael Dunlap


fig1.png

fig2.png



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #39473] more efficient algorithm for 3j, 6j, 9j calculations (gsl_sf_coupling_{3j, 6j, 9j}_e

[Patrick Alken]

Additional Item Attachment, bug #39473 (project gsl):

File name: report.pdf Size:33 KB


File name: coupling3j.patch   Size:15 KB




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #61094] gsl_linalg_LU_decomp fails to compute LU decomposition for some singular matrix

[Patrick Alken]

Update of bug #61094 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

Thank you for the report. This issue is fixed in commit
ed23edf3f88bf6b3d7d254521b19056bee6ce1ea

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: missing declarations in Basis Splines

[Patrick Alken]

Hello,

  I mistakenly updated the doc website with some experimental docs for 
some new B-spline features. These features are not in v2.7, or on the 
git yet. I have now updated the website with the correct 2.7 
documentation. Sorry for the trouble.


Thanks,

Patrick

On 7/3/21 3:01 PM, florian wrote:

Hi,

I could not run multiple examples from Basis Splines, e.g.: 
"https://www.gnu.org/software/gsl/doc/html/bspline.html#example-3-least-squares-and-breakpoints;, 
because
multiple documented functions could not be found in gsl include path 
(or gsl_bspline.h). Not found, but documented, functions are for example:

- gsl_bspline_ncontrol
- gsl_bspline_init_uniform
- gsl_bspline_wlssolve
- gsl_bspline_calc
- ...

Tested with gsl 2.7 from 
https://mirror.ibcp.fr/pub/gnu/gsl/gsl-latest.tar.gz and 
git://git.savannah.gnu.org/gsl.git

Tested with gsl 2.6.

Could not find any possible configure enable option within INSTALL or 
README files or from documentation.


Best regards
Flo

Re: Undefined reference to `gsl_spmatrix_uchar_norm1'

[Patrick Alken]

For now you can comment out that line. I will log this issue into the 
bug tracker to make a more permanent fix for the next release.


Thanks,

Patrick

On 6/28/21 10:20 AM, Håkon Hægland wrote:
[[ I am resending this since I forgot to select "Reply all", and 
bug-gsl@gnu.org <mailto:bug-gsl@gnu.org> was by accident not included 
in the recepients ]]
Yes, I agree that seems like a good reason to exclude those functions. 
Otherwise, I have no strong opinion on whether to include it or not.

 I am trying to write Perl bindings to the GSL library using swig, see

https://github.com/leto/math--gsl/issues/231 
<https://github.com/leto/math--gsl/issues/231>


According to the gsl_spmatrix_uchar.h header file (line 175):

https://git.savannah.gnu.org/cgit/gsl.git/tree/spmatrix/gsl_spmatrix_uchar.h#n175 
<https://git.savannah.gnu.org/cgit/gsl.git/tree/spmatrix/gsl_spmatrix_uchar.h#n175>


there should be a function gsl_spmatrix_uchar_norm1(), so swig gets 
confused and creates a wrapper for a non-existing function.


Maybe that line 175 could be commented out?

Håkon

On Fri, Jun 25, 2021 at 6:46 PM Patrick Alken <mailto:al...@colorado.edu>> wrote:


Hello,

   Unsigned versions of this function were excluded (uchar, uint,
ushort), because there is the possibility of adding two large uchar
values together during the 1-norm calculation and ending up with a
smaller value due to it wrapping back to 0. This didn't make much
sense
to me when I wrote the function, so I excluded the unsigned versions.

Do you think it makes sense to have a uchar_norm1 function?

Patrick

On 6/25/21 11:45 AM, Håkon Hægland wrote:
> I am on Ubuntu 21.04, using GSL 2.7. Here is my test program
norm.c :
>
>
> #include 
> #include 
> int
> main()
> {
>      gsl_spmatrix_uchar *A = gsl_spmatrix_uchar_alloc(2, 2);
>
>      gsl_spmatrix_uchar_set(A, 0, 0, 1);
>      gsl_spmatrix_uchar_set(A, 0, 1, 2);
>      gsl_spmatrix_uchar_set(A, 1, 0, 3);
>      gsl_spmatrix_uchar_set(A, 1, 1, 4);
>
>      double norm = gsl_spmatrix_uchar_norm1(A);
>      printf("norm = %g\n", norm);
>      gsl_spmatrix_uchar_free(A);
>
>      return 0;
> }
>
> When I compile this I get an error message from the linker:
>
> $ gcc -Wall -I/opt/gsl/gsl-2.7/include -c norm.c
> $ gcc -L/opt/gsl/gsl-2.7/lib -o norm norm.o -lgsl -lgslcblas -lm
> /bin/ld: norm.o: in function `main':
> norm.c:(.text+0x93): undefined reference to
`gsl_spmatrix_uchar_norm1'
> collect2: error: ld returned 1 exit status
>
> Best regards,
> Håkon Hægland

Re: Undefined reference to `gsl_spmatrix_uchar_norm1'

[Patrick Alken]

Hello,

  Unsigned versions of this function were excluded (uchar, uint, 
ushort), because there is the possibility of adding two large uchar 
values together during the 1-norm calculation and ending up with a 
smaller value due to it wrapping back to 0. This didn't make much sense 
to me when I wrote the function, so I excluded the unsigned versions.


Do you think it makes sense to have a uchar_norm1 function?

Patrick

On 6/25/21 11:45 AM, Håkon Hægland wrote:

I am on Ubuntu 21.04, using GSL 2.7. Here is my test program norm.c :


#include 
#include 
int
main()
{
 gsl_spmatrix_uchar *A = gsl_spmatrix_uchar_alloc(2, 2);

 gsl_spmatrix_uchar_set(A, 0, 0, 1);
 gsl_spmatrix_uchar_set(A, 0, 1, 2);
 gsl_spmatrix_uchar_set(A, 1, 0, 3);
 gsl_spmatrix_uchar_set(A, 1, 1, 4);

 double norm = gsl_spmatrix_uchar_norm1(A);
 printf("norm = %g\n", norm);
 gsl_spmatrix_uchar_free(A);

 return 0;
}

When I compile this I get an error message from the linker:

$ gcc -Wall -I/opt/gsl/gsl-2.7/include -c norm.c
$ gcc -L/opt/gsl/gsl-2.7/lib -o norm norm.o -lgsl -lgslcblas -lm
/bin/ld: norm.o: in function `main':
norm.c:(.text+0x93): undefined reference to `gsl_spmatrix_uchar_norm1'
collect2: error: ld returned 1 exit status

Best regards,
Håkon Hægland

Re: Inaccurate Results for Lambert W function

[Patrick Alken]

Thank you, I have logged the issue into the bug tracker

On 6/6/21 10:56 AM, Daming Zou wrote:

Hi,

The gsl_sf_lambert_W0(x) function gives inaccurate results for inputs 
near 0. I attached a code snippet to demonstrate it at the end of this 
mail. The details are as following:


[Error Examples]
- For very small input, e.g., 1e-30, 1e-40, the returned results are 
inaccurate;

- For smaller ones, e.g., 1e-50, 1e-60, the returned results are 0.

[Accurate Result]
- The series of lambertW(x) near 0 is W(x) = x - x^2 + 3/2 x^3 - 8/3 
x^4 + 125/24 x^5 + O(x^6), for more terms, please refer the Table 4 of 
this paper[1].
- For inputs |x|<1e-20 (more precisely, the double epsilon 
GSL_DBL_EPSILON), and for double precision, we can consider W(x) = x, 
since the -x^2 term can be neglected. Thus, this applies to 1e-30, 
1e-40, etc.


[Reason for Inaccurate Results]
- Currently, for input values larger than -1/E+1e-3, 
gsl_sf_lambert_W0_e(x) uses *halley_iteration* for solving the result. 
However, for extremely small values near 0, there is severe 
cancellation at line 61 of lambert.c : `w -= t;`, so the results are 
inaccurate or even 0.


[Possible Solution]
- Just return x for |x| < 1e-20.
- More aggressively, using polynomial to approximate the function for 
|x| < 1e-4. (A 5-terms polynomial should be sufficient for double 
precision). The interval and the number of terms could be further 
investigated.


[Reference]
[1] Fukushima, Toshio. "Precise and fast computation of Lambert 
W-functions without transcendental function evaluations." 
https://www.sciencedirect.com/science/article/pii/S0377042712005213


Not sure if this is defined as a bug in GSL, since from the aspect of 
absolute error, it may be acceptable; but from the aspect of relative 
error, it is significant. However, I believe the result W(1e-50) = 
1e-50 is better than W(1e-50) = 0 in most scenarios.


Thanks,
Daming


[Attached Code Example]

#include 
#include 

const double c[5] = {
 1.0,
    -1.0,
 3.0/2.0,
    -8.0/3.0,
 125.0/24.0
};

double horner_eval(double x) {
    double y = x*(c[0] + x*(c[1] + x*(c[2] + x*(c[3] + x*c[4];
    return y;
}

double inputs[7] = {
    1e-5,
    1e-10,
    1e-20,
    1e-30,
    1e-40,
    1e-50,
    1e-60
};

int main() {
    double x, y, y_poly;
    for (int i = 0; i < 7; i++) {
    x = inputs[i];
    y = gsl_sf_lambert_W0(x);
    y_poly = horner_eval(x);
    printf("input:  %.20e\n", x);
    printf("gsl_sf_lambertW:    %.20e\n", y);
    printf("poly_output:    %.20e\n", y_poly);
    printf("\n");
    }
    return 0;
}


$ ./gsl-lambert.out
input:  1.8180e-05
gsl_sf_lambertW:    9.149997397560e-06
poly_output:    9.149997397560e-06

input:  1.3643e-10
gsl_sf_lambertW:    9.89974982e-11
poly_output:    9.89974982e-11

input:  9.99945153e-21
gsl_sf_lambertW:    9.99945153e-21
poly_output:    9.99945153e-21

input:  1.8334e-30
gsl_sf_lambertW:    1.000275223352e-30
poly_output:    1.8334e-30

input:  9.99929293e-41
gsl_sf_lambertW:    1.03314933177740113005e-40
poly_output:    9.99929293e-41

input:  1.0762e-50
gsl_sf_lambertW:    0.e+00
poly_output:    1.0762e-50

input:  9.99970433e-61
gsl_sf_lambertW:    0.e+00
poly_output:    9.99970433e-61

[bug #60741] Inaccurate Results for Lambert W function

[Patrick Alken]

URL:
  

 Summary: Inaccurate Results for Lambert W function
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Mon 07 Jun 2021 03:49:16 AM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from daming.zou =at= inf =dot= ethz =dot= ch

Hi,

The gsl_sf_lambert_W0(x) function gives inaccurate results for inputs near 0.
I attached a code snippet to demonstrate it at the end of this mail. The
details are as following:

[Error Examples]
- For very small input, e.g., 1e-30, 1e-40, the returned results are
inaccurate;
- For smaller ones, e.g., 1e-50, 1e-60, the returned results are 0.

[Accurate Result]
- The series of lambertW(x) near 0 is W(x) = x - x^2 + 3/2 x^3 - 8/3 x^4 +
125/24 x^5 + O(x^6), for more terms, please refer the Table 4 of this
paper[1].
- For inputs |x|<1e-20 (more precisely, the double epsilon GSL_DBL_EPSILON),
and for double precision, we can consider W(x) = x, since the -x^2 term can be
neglected. Thus, this applies to 1e-30, 1e-40, etc.

[Reason for Inaccurate Results]
- Currently, for input values larger than -1/E+1e-3, gsl_sf_lambert_W0_e(x)
uses *halley_iteration* for solving the result. However, for extremely small
values near 0, there is severe cancellation at line 61 of lambert.c : `w -=
t;`, so the results are inaccurate or even 0.

[Possible Solution]
- Just return x for |x| < 1e-20.
- More aggressively, using polynomial to approximate the function for |x| <
1e-4. (A 5-terms polynomial should be sufficient for double precision). The
interval and the number of terms could be further investigated.

[Reference]
[1] Fukushima, Toshio. "Precise and fast computation of Lambert W-functions
without transcendental function evaluations."
https://www.sciencedirect.com/science/article/pii/S0377042712005213

Not sure if this is defined as a bug in GSL, since from the aspect of absolute
error, it may be acceptable; but from the aspect of relative error, it is
significant. However, I believe the result W(1e-50) = 1e-50 is better than
W(1e-50) = 0 in most scenarios.

Thanks,
Daming


[Attached Code Example]

#include 
#include 

const double c[5] = {
 1.0,
-1.0,
 3.0/2.0,
-8.0/3.0,
 125.0/24.0
};

double horner_eval(double x) {
double y = x*(c[0] + x*(c[1] + x*(c[2] + x*(c[3] + x*c[4];
return y;
}

double inputs[7] = {
1e-5,
1e-10,
1e-20,
1e-30,
1e-40,
1e-50,
1e-60
};

int main() {
double x, y, y_poly;
for (int i = 0; i < 7; i++) {
x = inputs[i];
y = gsl_sf_lambert_W0(x);
y_poly = horner_eval(x);
printf("input:  %.20e\n", x);
printf("gsl_sf_lambertW:%.20e\n", y);
printf("poly_output:%.20e\n", y_poly);
printf("\n");
}
return 0;
}


$ ./gsl-lambert.out
input:  1.8180e-05
gsl_sf_lambertW:9.149997397560e-06
poly_output:9.149997397560e-06

input:  1.3643e-10
gsl_sf_lambertW:9.89974982e-11
poly_output:9.89974982e-11

input:  9.99945153e-21
gsl_sf_lambertW:9.99945153e-21
poly_output:9.99945153e-21

input:  1.8334e-30
gsl_sf_lambertW:1.000275223352e-30
poly_output:1.8334e-30

input:  9.99929293e-41
gsl_sf_lambertW:1.03314933177740113005e-40
poly_output:9.99929293e-41

input:  1.0762e-50
gsl_sf_lambertW:0.e+00
poly_output:1.0762e-50

input:  9.99970433e-61
gsl_sf_lambertW:0.e+00
poly_output:9.99970433e-61








___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: Describe the error

[Patrick Alken]

Thanks, I fixed this bug on the git

On 5/19/21 9:13 AM, Mark Galassi wrote:
>> The descriptions of function gsl_histogram_max_bin and
>> function gsl_histogram_min_bin have something in common.
>> Is that a problem?
> I think you mean to say "the description of gsl_histogram_min_bin has an 
> error: it says maximum instead of minimum in one place".  Looking at it I 
> think you might be right.  This patch should fix it:
>
> diff --git a/doc/histogram.rst b/doc/histogram.rst
> index c25e342a..f5287b60 100644
> --- a/doc/histogram.rst
> +++ b/doc/histogram.rst
> @@ -322,7 +322,7 @@ Histogram Statistics
>  .. function:: size_t gsl_histogram_min_bin (const gsl_histogram * h)
>  
> This function returns the index of the bin containing the minimum
> -   value. In the case where several bins contain the same maximum value the
> +   value. In the case where several bins contain the same minimum value the
> smallest index is returned.
>  
>  .. index::
>
>
> Your other question:
>
>> And  I am trying to translate the manual into Chinese.
>> But for me, some words are difficult to understand. Could
>> you please help me to proofread them?
> It's really wonderful that you are translating the manual into Chinese.  
> Unfortunately I cannot help you because I don't know how to read Chinese, but 
> I hope someone else can.
>

Re: Problem with Schur decomposition

[Patrick Alken]

Note that the documentation says:

"If T is desired, it is stored in the upper portion of A on output".

This means that the gsl_eigen_nonsymm / gsl_eigen_nonsymm_Z functions do
not zero out the lower portion of the output matrix. The lower portion
of A is used as temporary workspace by the eigenvalue solver, so the
values stored here should not be referenced. You can explicitly set the
lower portion of A to zero with the function:

gsl_linalg_hessenberg_set_zero()

to get the behavior you want.

Patrick

On 5/18/21 3:05 PM, david allen wrote:
> What is supposed to be a upper quasi-triangular matrix has non-zero
> elements below the first sub-diagonal. This is demonstrated in the
> attached file. Computer and OS information are also in the file.
>

[bug #60635] physical constants may need updating

[Patrick Alken]

URL:
  

 Summary: physical constants may need updating
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 19 May 2021 02:23:28 AM UTC
Category: Accuracy problem
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from M. Galassi


I happened to come across a definition of the gravitational constant, which we
have as:

usr/include/gsl/gsl_const_cgs.h:#define GSL_CONST_CGS_GRAVITATIONAL_CONSTANT
(6.673e-8) /* cm^3 / g s^2 */
/usr/include/gsl/gsl_const_cgsm.h:#define
GSL_CONST_CGSM_GRAVITATIONAL_CONSTANT (6.673e-8) /* cm^3 / g s^2 */
/usr/include/gsl/gsl_const_mksa.h:#define
GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT (6.673e-11) /* m^3 / kg s^2 */
/usr/include/gsl/gsl_const_mks.h:#define GSL_CONST_MKS_GRAVITATIONAL_CONSTANT
(6.673e-11) /* m^3 / kg s^2 */

and the one I saw was different: 6.67430 m^3/(kg*s^2).

Looking it up, it turns out that since 1998 (when our number was valid), the
measurement has been updated.  You can see the history of updates at:

https://en.wikipedia.org/wiki/Gravitational_constant#Modern_value

and indeed the most recent is 6.67430.

We might want to do an audit to make sure that we follow the NIST or LBL or
another standards body on *all* the physical constants, since you want to
update them in lockstep to make sure that some identities still apply.  Our
constants were introduced around 2000 and have not been touched since some
additions in 2006.

At the same time we might want to do a bit of soul searching on where our
cutoff is for accepting constants into GSL.  I do some work with earth
parameters, and am shocked at the lack of a clear reference set of C header
files (from NASA or whoever) with earth constants.

gsl currently has fundamental constants (boltzmann, gravitational, bohr
radius, electron charge, ...  But we also have some more solar-system bound
ones, like the astronomical unit, and we even reach pedestrian earth with
things like GSL_CONST_MKSA_GRAM_FORCE which relates to the acceleration of
gravity on the earth's surface.  In fact we also have
GSL_CONST_MKSA_GRAV_ACCEL, which is 9.80665 m/s^2.

So I'd like to propose that we add a suite of earth radius values based on
concepts from:

https://en.wikipedia.org/wiki/Earth_radius#Global_average_radii

https://gis.stackexchange.com/questions/25494/how-accurate-is-approximating-the-earth-as-a-sphere

and fact sheets like:

https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html (and the
references they give)

To start with the earth radius, 

The International Union of Geodesy and Geophysics (IUGG) defines various types
of mean earth radius, and they tend to stem from the definition of earth
equatorial radius A and earth polar radius B.  Setting R1 might be enough to
cover what most people need.

A patch for the gravity constant and addition of the earth stuff I mentioned
would be like below, and similar for mks, cgs, and cgsm.

At this point I'm just proposing it as food for thought - there is no urgency,
and it might be best to apply this patch when we have audited *all* the
physical constants we use, and possibly documented their provenance.

diff --git a/const/gsl_const_mksa.h b/const/gsl_const_mksa.h
index 5d91d1ca..4c82e272 100644
--- a/const/gsl_const_mksa.h
+++ b/const/gsl_const_mksa.h
@@ -22,12 +22,15 @@
 #define __GSL_CONST_MKSA__
 
 #define GSL_CONST_MKSA_SPEED_OF_LIGHT (2.99792458e8) /* m / s */
-#define GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT (6.673e-11) /* m^3 / kg s^2 */
+#define GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT (6.67430-11) /* m^3 / (kg s^2)
*/
 #define GSL_CONST_MKSA_PLANCKS_CONSTANT_H (6.62606896e-34) /* kg m^2 / s */
 #define GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR (1.05457162825e-34) /* kg m^2 /
s */
 #define GSL_CONST_MKSA_ASTRONOMICAL_UNIT (1.49597870691e11) /* m */
 #define GSL_CONST_MKSA_LIGHT_YEAR (9.46053620707e15) /* m */
 #define GSL_CONST_MKSA_PARSEC (3.08567758135e16) /* m */
+#define GSL_CONST_MKSA_EARTH_EQUATORIAL_RADIUS_A (6378137.0) /* m */
+#define GSL_CONST_MKSA_EARTH_POLAR_RADIUS_B (6356752.3) /* m */
+#define GSL_CONST_MKSA_EARTH_MEAN_RADIUS_R1
((2*GSL_CONST_MKS_EARTH_EQUATORIAL_RADIUS_A +
GSL_CONST_MKS_EARTH_POLAR_RADIUS_B) / 3.0) /* m */
 #define GSL_CONST_MKSA_GRAV_ACCEL (9.80665e0) /* m / s^2 */
 #define GSL_CONST_MKSA_ELECTRON_VOLT (1.602176487e-19) /* kg m^2 / s^2 */
 #define GSL_CONST_MKSA_MASS_ELECTRON (9.10938188e-31) /* kg */






___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: Bug in the complex no header file of the gsl_complex.h .

[Patrick Alken]

Try this instead:

gsl_complex z;

GSL_SET_COMPLEX(,5.7,8.0);
printf("output %f",gsl_complex_arg(z));

On 5/13/21 12:50 AM, Chetan wrote:
>    error: expected identifier or ‘(’ before ‘do’
>        95 | #define GSL_SET_COMPLEX(zp,x,y) do {(zp)->dat[0]=(x);
>    (zp)->dat[1]=(y);} while(0)
>
>    here is the piece of code where i was trying to use the library.
>
>         gsl_complex z;
>         gsl_complex GSL_SET_COMPLEX(,5.7,8.0);
>
>         printf("output %f",gsl_complex_arg(z));
>
>    Thank you
>

Re: [bug #60457] Feature request: complex tridiagonal solvers

[Patrick Alken]

Hello, certainly I can add that function. It would be best if you send
me a patch with 'git diff'

Patrick

On 4/26/21 12:23 PM, Ed wrote:
> URL:
>   
>
>  Summary: Feature request: complex tridiagonal solvers
>  Project: GNU Scientific Library
> Submitted by: mohawk
> Submitted on: Mon 26 Apr 2021 06:23:56 PM UTC
> Category: None
> Severity: 3 - Normal
> Operating System: 
>   Status: None
>  Assigned to: None
>  Open/Closed: Open
>  Release: 2.6
>  Discussion Lock: Any
>
> ___
>
> Details:
>
> I was most surprised to find that while there is the real
> gsl_linalg_solve_tridiag, there is no complex equivalent. Would you like me to
> make a pull request to add it (and maybe other complex tridiagonal functions)?
>
>
>
>
> ___
>
> Reply to this item at:
>
>   
>
> ___
>   Message sent via Savannah
>   https://savannah.gnu.org/
>
>

[bug #55412] Bessel function returns incorrect value

[Patrick Alken]

Update of bug #55412 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

Note that the documentation uses the cylindrical Bessel function
gsl_sf_bessel_J0, while your program uses the spherical Bessel function
gsl_sf_bessel_j0. Both outputs are correct, and the documentation is correct.

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60370] Interpolation inconsistencies

[Patrick Alken]

Update of bug #60370 (project gsl):

  Status:None => Duplicate  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

This is a duplicate, and has been fixed in a previous commmit. The e_extrap
function is marked as deprecated

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60371] Interpolation domain error handling

[Patrick Alken]

URL:
  

 Summary: Interpolation domain error handling
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Sun 11 Apr 2021 03:59:07 AM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

reported by Damien Lebrun-Grandie dalg24 =at= gmail =dot= com

With the current version of GSL (without the patch), the user do not get a
chance to check the return status of gsl_interp2d_eval_e() because GSL_ERROR()
gets called direclty instead of returning an error code as advertised in the
documentation.

With the changes I propose, the domain error will still be caught by macro
DISCARD_STATUS() that checks the return code of gsl_interp2d_eval_e() in
gsl_interp2d_eval() but the intended behavior for gsl_interp2d_eval_e() will
be achieved.  The only downside I see with this patch is that the user now
gets an out-of-domain error but he is not told anymore whether the error was
caused by the x or the y coordinate.





___

File Attachments:


---
Date: Sun 11 Apr 2021 03:59:07 AM UTC  Name:
bug_interp2d_domain_error_handling.c  Size: 2KiB   By: psa


---
Date: Sun 11 Apr 2021 03:59:07 AM UTC  Name:
0001-Fix-error-handling-in-2D-interpolation-functions.patch  Size: 1KiB   By:
psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60370] Interpolation inconsistencies

[Patrick Alken]

URL:
  

 Summary: Interpolation inconsistencies
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Sun 11 Apr 2021 03:54:57 AM UTC
Category: Documentation
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

reported by Damien Lebrun-Grandie dalg24 =at= gmail =dot= com

The function gsl_interp2d_eval_e_extrap is misnamed. It should be:
gsl_interp2d_eval_extrap_e. The documentation is correct but the function name
is wrong.




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60335] Test failure in spmatrix

[Patrick Alken]

Update of bug #60335 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

spmatrix test failure fixed in commit
15ea0bf6bb19c048010fed2879ff2669ed37a1b7

Note that gsl_linalg_hessenberg is deprecated and is replaced by
gsl_linalg_hessenberg_decomp

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60335] Test failure in spmatrix

[Patrick Alken]

URL:
  

 Summary: Test failure in spmatrix
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Sun 04 Apr 2021 05:13:07 AM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from J.D.Lamb =at= johndlamb =dot= net

I notice that gsl_linalg_hessenberg is still not declared in any header file.
Presumably the intention was to declare it in gsl_linalg.h

Also I get the following fails in spmatrix:
Linux 4.15.0-58-generic #64-Ubuntu SMP Tue Aug 6 11:12:41 UTC 2019 x86_64
x86_64 x86_64 GNU/Linux
gcc (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0
CFLAGS = -march=native -O2 -malign-double -Wall -pipe -pedantic
-Wno-long-long

FAIL: gsl_spmatrix_complex_scale_columns[53,107](COO) [1132]
FAIL: gsl_spmatrix_complex_scale_rows[53,107](COO) [1133]
FAIL: gsl_spmatrix_complex_scale_columns[53,107](CSC) [1135]
FAIL: gsl_spmatrix_complex_scale_rows[53,107](CSC) [1136]
FAIL: gsl_spmatrix_complex_scale_columns[53,107](CSR) [1138]
FAIL: gsl_spmatrix_complex_scale_rows[53,107](CSR) [1139]
FAIL: gsl_spmatrix_complex_float_scale[53,107](COO) [1211]
FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](COO) [1212]
FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](COO) [1213]
FAIL: gsl_spmatrix_complex_float_scale[53,107](CSC) [1214]
FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](CSC) [1215]
FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](CSC) [1216]
FAIL: gsl_spmatrix_complex_float_scale[53,107](CSR) [1217]
FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](CSR) [1218]
FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](CSR) [1219]
FAIL: gsl_spmatrix_complex_scale_columns[40,20](COO) [2450]
FAIL: gsl_spmatrix_complex_scale_rows[40,20](COO) [2451]
FAIL: gsl_spmatrix_complex_scale_columns[40,20](CSC) [2453]
FAIL: gsl_spmatrix_complex_scale_rows[40,20](CSC) [2454]
FAIL: gsl_spmatrix_complex_scale_columns[40,20](CSR) [2456]
FAIL: gsl_spmatrix_complex_scale_rows[40,20](CSR) [2457]
FAIL: gsl_spmatrix_complex_float_scale[40,20](COO) [2529]
FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](COO) [2530]
FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](COO) [2531]
FAIL: gsl_spmatrix_complex_float_scale[40,20](CSC) [2532]
FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](CSC) [2533]
FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](CSC) [2534]
FAIL: gsl_spmatrix_complex_float_scale[40,20](CSR) [2535]
FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](CSR) [2536]
FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](CSR) [2537]
FAIL: gsl_spmatrix_complex_scale_columns[30,30](COO) [3768]
FAIL: gsl_spmatrix_complex_scale_rows[30,30](COO) [3769]
FAIL: gsl_spmatrix_complex_scale_columns[30,30](CSC) [3771]
FAIL: gsl_spmatrix_complex_scale_rows[30,30](CSC) [3772]
FAIL: gsl_spmatrix_complex_scale_columns[30,30](CSR) [3774]
FAIL: gsl_spmatrix_complex_scale_rows[30,30](CSR) [3775]
FAIL: gsl_spmatrix_complex_float_scale[30,30](COO) [3847]
FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](COO) [3848]
FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](COO) [3849]
FAIL: gsl_spmatrix_complex_float_scale[30,30](CSC) [3850]
FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](CSC) [3851]
FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](CSC) [3852]
FAIL: gsl_spmatrix_complex_float_scale[30,30](CSR) [3853]
FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](CSR) [3854]
FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](CSR) [3855]







___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #36577] vector/matrix _add, _scale routines take wrong scalar type

[Patrick Alken]

Update of bug #36577 (project gsl):

  Status:   Confirmed => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

fixed in commit a0f31bd68dbeb467669d72bd5d53caeb8f310f1c

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #60026] Incorporate MIXMAX random number extension into GSL

[Patrick Alken]

URL:
  

 Summary: Incorporate MIXMAX random number extension into GSL
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Tue 09 Feb 2021 04:09:05 AM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from savvidy =at= inp =dot= demokritos =dot= gr

We have news on MIXMAX PRNG. 

1.It was fully implemented at CERN into the CLHEP/Geant4 software
a toolkit for the simulation of the passage of particles through matter as
default generator:

https://proj-clhep.web.cern.ch/proj-clhep/(see Notes)

2.Became a part of the CMS Simulation Software at CERN:


https://indico.cern.ch/event/731433/contributions/3015654/attachments/1680131/2698971/CMSsim.pdf
  (see Pages 13-17)

https://indico.cern.ch/event/587955/contributions/2937635/attachments/1679273/2706817/PosterCMS_SIM_v4.pdf


3.It was recently used by NASA  in their design of the new mission to the Sun
with the neutrino telescope  - a Solar Neutrino Spacecraft - the article
is attached (see page 13)

The authors of this project stated that: 

" Early in this project, the random number generation functions available in
standard C++ packages proved unsuitable for use;

they produce values which fall into a regular pattern of magnitude and will
generally avoid generating numbers very near zero or one.

For these reasons, the MIXMAX family of pseudo-random number generators have
been utilised here. Based off on Anosov C-systems—

a mathematical function group with useful instabilities in their structure—
this family of generators allows for the most random numerical

output in the shortest computational time [7] ".

We think that the problems with C++ packages above were connected with the
resolution power of the generators.
In a simple terms -  the MIXMAX is genuine 64 bit generator and have superior
resolution power. 

Please let me know if there were questions regarding MIXMAX.




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: gsl_ran_gaussian_ziggurat doesn't work with gsl_rng_slatec and gsl_rng_uni

[Patrick Alken]

Hello,

  I would welcome your help with fixing rng issues. We don't currently
have a GSL developer with expertise in rngs (I myself have never worked
on them). The best approach would be to clone the git repository
(instructions on the GSL webpage), make your changes, and then send me a
git diff when ready. Then I can take a look at things.

Best,
Patrick

On 2/5/21 12:52 AM, camel-cdr--- via Bug reports for the GNU Scientific
Library wrote:
> I've been developing a ziggurat implementation and used 
> https://www.seehuhn.de/pages/ziggurat to benchmark mine against the one used 
> in gsl.
> I ran[timegauss.c,](https://m.seehuhn.de/data/ziggurat/timegauss.c) but it 
> fails with gsl_rng_slatec and gsl_rng_uni.
> Consequently, I checked randist/gausszig.c and gsl_ran_gaussian_ziggurat 
> doesn't support generators with smaller ranges than 16777216.
>
> On a related note: I'd also be interested in improving the performance of 
> gsl_ran_gaussian_ziggurat outside of this bug, but I'm new to this project, 
> and I'm not sure what the procedure would be. (In my current benchmarks my 
> implementation is 1.3 times faster and uses 2/3 off the memory)

[bug #45521] Erroneous use of GSL_ERROR_NULL instead of GSL_ERROR macro

[Patrick Alken]

Update of bug #45521 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

fixed in commit adc8360d469bfa7e6e22a563baa9d011b724

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59910] Update information for Python wrapper

[Patrick Alken]

Update of bug #59910 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

added link for CythonGSL to webpage

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59494] Bug in gsl_rstat

[Patrick Alken]

Update of bug #59494 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

The issue with n=5 has been corrected in commit
42933198abec39071c1ae1bdf2d7e3adeed233a2

For n > 5, the median calculation only produces an estimate of the true
median, not the exact median, using the algorithm of Jain and Chlamtec 1985.
The only way to compute the exact median is to store all values in an array
and use gsl_stats_median. So this is not in fact a bug in the code.

I also uploaded bug_corrected.c, which fixes some syntax errors in the
original example program.

(file #50758)
___

Additional Item Attachment:

File name: bug_corrected.cSize:1 KB




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59758] Doc bug

[Patrick Alken]

Update of bug #59758 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

fixed in commit c1a38a563eec07eab392a99efc2cd5876ff62c68

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59710] Undefined behavior error in Nonlinear Least-Squares Fitting

[Patrick Alken]

Update of bug #59710 (project gsl):

  Status:None => Fixed  


___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #58202] Wrong Result Median for gsl_rstat_quantile_get with fix

[Patrick Alken]

Update of bug #58202 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

fixed in commit 42933198abec39071c1ae1bdf2d7e3adeed233a2

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59914] Native build of GSL-2.5 on windows 10

[Patrick Alken]

URL:
  

 Summary: Native build of GSL-2.5 on windows 10
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 20 Jan 2021 02:43:53 AM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from hariseldon99 =at= gmail =dot= com

I've managed to build gsl-2.5 dlls for windows 10 by cross-compiling on my
ubuntu install using MXE
.
I had to do this for my
masters' teaching lab, as well as colleges affiliated with my university,
all of whom only have windows machines. Details are given below

https://sites.google.com/phys.buruniv.ac.in/programming-workshop-cbcs/setup

I tried to use the methods linked to at the gsl website (@
https://www.gnu.org/software/gsl/extras/native_win_builds.html) but failed
for various reasons.

If my way is deemed useful or worthy, then perhaps I could create a git
page for this and have it linked from the gsl website?




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59913] gsl 2.3.0 problem in gsl_integration_cquad

[Patrick Alken]

URL:
  

 Summary: gsl 2.3.0 problem in gsl_integration_cquad
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 20 Jan 2021 02:41:50 AM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from buchholz =at= tugraz =dot= at

I found a problem in gsl (assuming the compiler/linker is not broken), which I
would consider a bug.

A minimal working example is attached.
I compiled using g++ 6.3.0 and gsl 2.3.0. The compile command is given in the
file.
The code uses an extension to cause the code to stop on floating point
exceptions (found this via stack exchange
https://linux.die.net/man/3/feenableexcept).

I would expect the code to just run through with output something like

result: 4.39016e-164
abserr: 0
nevals: 33

Instead i get the message

Gleitkomma-Ausnahme (german for floating point exception)

Running the program with gdb I get

Program received signal SIGFPE, Arithmetic exception.
0x77a1b42f in gsl_integration_cquad () from
/usr/lib/x86_64-linux-gnu/libgsl.so.19
(gdb) backtrace
#0  0x77a1b42f in gsl_integration_cquad () from
/usr/lib/x86_64-linux-gnu/libgsl.so.19
#1  0x4cee in main () at main.cpp:42


Changing the lower bound to 19.2 will cause the code to work, i.e. no floating
point exception.

If I do not trap the exception the code works as expected, but if there is no
real problem, then no exception should be raised.






___

File Attachments:


---
Date: Wed 20 Jan 2021 02:41:50 AM UTC  Name: main.cpp  Size: 1KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59912] gsl_permutation header files

[Patrick Alken]

URL:
  

 Summary: gsl_permutation header files
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 20 Jan 2021 02:38:43 AM UTC
Category: Documentation
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from brijesh =dot= upadhaya =at= gmail =dot= com

The part of gsl-2.5 documentation that describes the functions related to
permutations says
"The functions described in this chapter are defined in the header file
gsl_permutation.h"

However, I got a warning: implicit declaration of function
‘gsl_permute_matrix’, thus, I looked into the header files in
/include/gsl
and found out that gsl_permutation.h didn't include gsl_permute_matrix.h.

This probably is not a bug as I added gsl_permute_matrix.h and compiled my
program. However, it wasn't mentioned in the online documentation.





___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59911] Problem with qagui 1D integrator

[Patrick Alken]

URL:
  

 Summary: Problem with qagui 1D integrator
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 20 Jan 2021 02:36:50 AM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from takacs =dot= adam =at= wigner =dot= mta =dot= hu

Dear GSL Helpers,

We tried to use the qagiu algorithm to calculate the following integration

   \int_0^infty dx cosh(x) e^{-cosh(x)} ,

however qagiu falls and complains on "bad integrand behavior". If we're
using qags with a fixed big upper limit, the integral is fine.

Do you have any suggestion what could be the problem? We'd like to
use qagiu without upper limit.




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59910] Update information for Python wrapper

[Patrick Alken]

URL:
  

 Summary: Update information for Python wrapper
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 20 Jan 2021 02:34:41 AM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from david =dot= cortes =dot= rivera =at= gmail =dot= com

In the GNU page for GSL: https://www.gnu.org/software/gsl/
the section for wrappers suggests PyrexGsl as a Pyrex wrapper. However,
that package hasn't been updated since a long time (same for pyrex) and
doesn't cover newer additions to GSL. It would be more appropriate to
point towards its fork CythonGSL (https://github.com/twiecki/CythonGSL)
, which works with cython (successor of pyrex).




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59900] Add truncated normal distribution

[Patrick Alken]

URL:
  

 Summary: Add truncated normal distribution
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Sun 17 Jan 2021 10:26:25 PM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from msunet =at= shellblade =dot= net

Please find attached an implementation of a truncated normal distribution
based on John Burkardt's presentation linked below.

Please note that I am neither a mathematician nor an expert in numerical
computing, so any feedback is appreciated. This is also my first contribution
to GSL (and a GNU project in general), so please excuse any inefficiencies.

Thank you,

Marc



Adds a truncated normal distribution as described in John Burkardt's
"The Truncated Normal Distribution":

https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
---
 cdf/gauss.c   | 22 
 cdf/gaussinv.c| 20 ++
 cdf/gsl_cdf.h |  6 +
 cdf/test.c| 61 +++
 filter/Makefile.am|  2 +-
 movstat/Makefile.am   |  2 +-
 randist/gauss.c   | 19 ++
 randist/gsl_randist.h |  3 +++
 randist/test.c| 20 ++
 rstat/Makefile.am |  2 +-
 10 files changed, 154 insertions(+), 3 deletions(-)

diff --git a/cdf/gauss.c b/cdf/gauss.c
index cdc8b0a3..110e58f1 100644
--- a/cdf/gauss.c
+++ b/cdf/gauss.c
@@ -349,3 +349,25 @@ gsl_cdf_gaussian_Q (const double x, const double sigma)
 {
   return gsl_cdf_ugaussian_Q (x / sigma);
 }
+
+double
+gsl_cdf_tgaussian_P (const double x, const double a, const double b, const
double sigma)
+{
+  if (x <= a)
+{
+  return 0.0;
+}
+  if (x >= b)
+{
+  return 1.0;
+}
+  const double n = gsl_cdf_gaussian_P (x, sigma) - gsl_cdf_gaussian_P (a,
sigma);
+  const double d = gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a,
sigma);
+  return n / d;
+}
+
+double
+gsl_cdf_tgaussian_Q (const double x, const double a, const double b, const
double sigma)
+{
+  return 1.0 - gsl_cdf_tgaussian_P (x, a, b, sigma);
+}
diff --git a/cdf/gaussinv.c b/cdf/gaussinv.c
index dcbdca1d..28c9866c 100644
--- a/cdf/gaussinv.c
+++ b/cdf/gaussinv.c
@@ -200,3 +200,23 @@ gsl_cdf_gaussian_Qinv (const double Q, const double
sigma)
 {
   return sigma * gsl_cdf_ugaussian_Qinv (Q);
 }
+
+double
+gsl_cdf_tgaussian_Pinv (const double P, const double a, const double b, const
double sigma)
+{
+  const double k
+= gsl_cdf_gaussian_P (a, sigma)
++ P * (gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma));
+
+  return gsl_cdf_gaussian_Pinv (k, sigma);
+}
+
+double
+gsl_cdf_tgaussian_Qinv (const double Q, const double a, const double b, const
double sigma)
+{
+  const double k
+= gsl_cdf_gaussian_P (b, sigma)
+- Q * (gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma));
+
+  return gsl_cdf_gaussian_Pinv (k, sigma);
+}
diff --git a/cdf/gsl_cdf.h b/cdf/gsl_cdf.h
index 2bc3fed5..61379a82 100644
--- a/cdf/gsl_cdf.h
+++ b/cdf/gsl_cdf.h
@@ -46,6 +46,12 @@ double gsl_cdf_gaussian_Q (const double x, const double
sigma);
 double gsl_cdf_gaussian_Pinv (const double P, const double sigma);
 double gsl_cdf_gaussian_Qinv (const double Q, const double sigma);
 
+double gsl_cdf_tgaussian_P (const double x, const double a, const double b,
const double sigma);
+double gsl_cdf_tgaussian_Q (const double x, const double a, const double b,
const double sigma);
+
+double gsl_cdf_tgaussian_Pinv (const double P, const double a, const double
b, const double sigma);
+double gsl_cdf_tgaussian_Qinv (const double Q, const double a, const double
b, const double sigma);
+
 double gsl_cdf_gamma_P (const double x, const double a, const double b);
 double gsl_cdf_gamma_Q (const double x, const double a, const double b);
 
diff --git a/cdf/test.c b/cdf/test.c
index 1b4cf777..65eeaa6e 100644
--- a/cdf/test.c
+++ b/cdf/test.c
@@ -39,6 +39,8 @@
 
 void test_ugaussian (void);
 void test_ugaussianinv (void);
+void test_tgaussian (void);
+void test_tgaussianinv (void);
 void test_exponential (void);
 void test_exponentialinv (void);
 void test_exppow (void);
@@ -398,6 +400,7 @@ main (void)
  Function values computed with PARI, 28 digits precision */
   
   test_ugaussian ();
+  test_tgaussian ();
   test_exponential ();
   test_exppow ();
   test_tdist (); 
@@ -407,6 +410,7 @@ main (void)
   test_beta (); 
 
   test_ugaussianinv ();
+  test_tgaussianinv ();
   test_exponentialinv ();
   test_gammainv (); 
   test_chisqinv (); 
@@ -535,6 +539,63 @@ void test_ugaussianinv (void) {
 }
 
 
+  /* Test values from John Burkardt, The Truncated Normal Distribution
+ Page 23. */
+
+void test_tgaussian 

[bug #59624] Buffer overflow in gsl_stats_quantile_from_sorted_data

[Patrick Alken]

Follow-up Comment #1, bug #59624 (project gsl):

from msunet =at= shellblade =dot= net

This one is a "fix" for bug 59624, simply checking that |f| is in the expected
range.
Though I don't think there is ultimately a way to check in C whether the
given
array has the right size, so a buggy application can still trigger the bug
described in the email thread.

---
 statistics/quantiles_source.c | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/statistics/quantiles_source.c b/statistics/quantiles_source.c
index e2956d9d..bf93a1a3 100644
--- a/statistics/quantiles_source.c
+++ b/statistics/quantiles_source.c
@@ -17,6 +17,7 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
USA.
  */
 
+#include 
 
 double
 FUNCTION(gsl_stats,quantile_from_sorted_data) (const BASE sorted_data[], 
@@ -24,6 +25,8 @@ FUNCTION(gsl_stats,quantile_from_sorted_data) (const BASE
sorted_data[],
const size_t n,
const double f)
 {
+  assert(0.0 <= f && f <= 1.0);
+
   const double index = f * (n - 1) ;
   const size_t lhs = (int)index ;
   const double delta = index - lhs ;
-- 
2.27.0



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59834] Misaligned memory access error in deque.c

[Patrick Alken]

URL:
  

 Summary: Misaligned memory access error in deque.c
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 06 Jan 2021 03:50:58 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

reported by zhoulai.fu =at= gmail =dot= com

Running doc/examples/filt_edge.c with the Undefined Sanitizer enabled
reveals many unaligned memory access errors in deque.c. (See
https://www.kernel.org/doc/Documentation/unaligned-memory-access.txt
for unaligned memory access.)

deque.c:58:11: runtime error: member access within misaligned address
0x024010f4 for type 'struct deque', which requires 8 byte alignment
0x024010f4: note: pointer points here
  00 00 00 00 00 00 00 00  00 00 00 00 00 00 00 00  00 00 00 00 00 00 00 00
 00 00 00 00 00 00 00 00
  ^
deque.c:59:11: runtime error: member access within misaligned address
0x024010f4 for type 'struct deque', which requires 8 byte alignment
0x024010f4: note: pointer points here
  00 00 00 00 ff ff ff ff  00 00 00 00 00 00 00 00  00 00 00 00 00 00 00 00
 00 00 00 00 00 00 00 00
  ^
deque.c:60:11: runtime error: member access within misaligned address
0x024010f4 for type 'struct deque', which requires 8 byte alignment
0x024010f4: note: pointer points here
  00 00 00 00 ff ff ff ff  00 00 00 00 00 00 00 00  00 00 00 00 00 00 00 00
 00 00 00 00 00 00 00 00
  ^
deque.c:61:12: runtime error: member access within misaligned address
0x024010f4 for type 'struct deque', which requires 8 byte alignment
0x024010f4: note: pointer points here
  00 00 00 00 ff ff ff ff  00 00 00 00 05 00 00 00  00 00 00 00 00 00 00 00
 00 00 00 00 00 00 00 00
  ^
...

To reproduce the error, compile GSL and filt_edge.c with
"-fsanitize=undefined" flag and run the binary.

See the attached error.txt for the full error messages.

*Preliminary diagnosis:*

Gdb shows the following backtrace at the place where the error occurs:

(gdb) bt
#0  deque_init (n=5, d=0x6d20f4) at deque.c:59
#1  0x0040cd71 in mmacc_init (n=4, vstate=0x6d2060) at mmacc.c:88
#2  0x00401d59 in rmedian_init (n=4, vstate=0x6d2050) at
rmedian.c:175
#3  0x00404aeb in gsl_movstat_apply_accum
(endtype=GSL_MOVSTAT_END_PADVALUE, x=0x7fffd960, accum=0x66d9c0
,
accum_params=0x7fffd950, y=0x7fffd990, z=0x0, w=0x6d2010) at
apply.c:74
#4  0x00401bba in gsl_filter_rmedian
(endtype=GSL_FILTER_END_PADVALUE, x=0x6d4100, y=0x6d8040, w=0x6d1eb0) at
rmedian.c:147
#5  0x0040108d in main () at filt_edge.c:38


That is, the error occurs at in the function

 static int
 deque_init(const size_t n, deque * d)

of deque.c. In the function, the compiler expects the  pointer of type
deque, d, to be 8-aligned. Yet d got an address 0x6d20f4 that is not
8-aligned.

The cause of the error should closely relate to mmann.c:84

  state->maxque = (deque *) ((unsigned char *) state->minque + deque_size(n
+ 1));

In this operation of pointer casting, if we assume that the first term
above, state->minque, is a multiple of 8, then, the second term above,
deque_size(n + 1), should be a multiple of 8 as well, in order to meet the
8-alignment requirement. Now, take a look at the function deque_size in
deque.c:

 static size_t
 deque_size(const size_t n)
 {
   size_t size = 0;

   size += sizeof(deque);
   size += n * sizeof(deque_type_t); /* array */

   return size;
 }

The function returns 24 + n * 4, since sizeof(deque) is determined to be
24 at compile time, and sizeof(deque_type_t) is 4 as deque_type_t is a
typedef of int. This number, 24 + n * 4, is a multiple of 8 if and only if
n is even. In other words, the misalignment should occur whenever n
is odd.

Note that the kind of
misaligned memory access can cause undefined behavior, triggering
performance
penalty, crash or incorrect results. See
https://developer.ibm.com/technologies/systems/articles/pa-dalign/.



___

File Attachments:


---
Date: Wed 06 Jan 2021 03:50:58 PM UTC  Name: error.txt  Size: 13KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59781] Eigenvalue failure

[Patrick Alken]

URL:
  

 Summary: Eigenvalue failure
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Tue 29 Dec 2020 11:38:25 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

The attached code demonstrates a failure of the non-symmetric eigenvalue
solver. Reported by Dmitry Cheshkov (dcheshkov =at= gmail =dot= com)



___

File Attachments:


---
Date: Tue 29 Dec 2020 11:38:25 PM UTC  Name: eig.c  Size: 2KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: gsl: francis.c:209: ERROR: maximum iterations reached without finding all eigenvalues

[Patrick Alken]

Dmitry,

  I confirm I get the same error using the matrix transpose. I will log
this in the bug tracker, since it might take me some time to diagnose
and fix the problem. In the meantime, perhaps you can use the
non-transposed matrix, since A and A^T have the same eigenvalues.

Patrick

On 12/29/20 3:15 PM, Dmitry Cheshkov wrote:
> Dear Patrick,
>
> I'm also using gsl 2.6. And the difference between your code and my is
> that I'm filled the matrix in transposed form. Then I tried to
> calculate eigenvalues for transposed matrix by your code, the error
> come back.
>
> Try to calculate eigenvalues of the folowing matrices by your code:
>
>   const double data[] = {
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 3, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 4,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 4, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 2, 0, 2,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 3, 0, 0, 0,
>    0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 2, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 2, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    1, 0, 1, 0, 0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
>
>
> another similar matrix:
>
>   const double data[] = {
>    0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0,
>    0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0,
>    0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0,
>    0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0,
>    0, 9, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 9, 0, 0, 0, 9,
>    9, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 9, 0, 0, 0, 9, 0,
>    0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0,
>    0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0,
>    0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>    0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
>
> Wolfram Mathematica:
> In[4]:= Eigenvalues[data]
> Out[4]= {-9, -9, -9, -9, 9, 9, 9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
>
>
> With best regards,
> Dmitry Cheshkov
>
> 29.12.2020 20:11, Patrick Alken пишет:
>> Dear Dmitry,
>>
>>    I successfully found the eigenvalues of your matrix using GSL 2.6
>> (see
>> attached code). The output of the program is:
>>
>> ---
>> $ ./eig
>> eval[1] = 0.e+00 + i 0.e+00
>> eval[2] = 0.e+00 + i 0.e+00
>> eval[3] = 0.e+00 + i 0.e+00
>> eval[4] = 4.300443739734e-16 + i 0.e+00
>> eval[5] = -6.e+00 + i 0.e+00
>> eval[6] = 6.e+00 + i 0.e+00
>> eval[7] = 6.e+00 + i 0.e+00
>> eval[8] = -6.e+00 + i 0.e+00
>> eval[9] = 5.377642775528e-17 + i 0.e+00
>> eval[10] = 3.e+00 + i 0.e+00
>> eval[11] = -3.e+00 + i 0.e+00
>> eval[12] = 6.e+00 + i 0.e+00
>> eval[13] = -6.e+00 + i 0.e+00
>> eval[14] = 3.e+00 + i 0.e+00
>> eval[15] = -2.912706738243e-16 + i 0.e+00
>> eval[16] = 6.e+00 + i 0.e+00
>> eval[17] = -6.e+00 + i 0.e+00
>> eval[18] = -3.e+00 + i 0.e+00
>> ---
>>
>> Can you clarify which version of GSL you are using? Also you might want
>> to print out your matrix m to make sure it is initialized correctly. I
>> did not try to debug your code.
>>
>> Best,
>> Patrick
>>
>> On 12/29/20 7:44 AM, Dmitry Cheshkov wrote:
>>> Dear Sirs!
>>>
>>> I have to calculate the eigenvalue

Re: gsl: francis.c:209: ERROR: maximum iterations reached without finding all eigenvalues

[Patrick Alken]

Dear Dmitry,

  I successfully found the eigenvalues of your matrix using GSL 2.6 (see
attached code). The output of the program is:

---
$ ./eig
eval[1] = 0.e+00 + i 0.e+00
eval[2] = 0.e+00 + i 0.e+00
eval[3] = 0.e+00 + i 0.e+00
eval[4] = 4.300443739734e-16 + i 0.e+00
eval[5] = -6.e+00 + i 0.e+00
eval[6] = 6.e+00 + i 0.e+00
eval[7] = 6.e+00 + i 0.e+00
eval[8] = -6.e+00 + i 0.e+00
eval[9] = 5.377642775528e-17 + i 0.e+00
eval[10] = 3.e+00 + i 0.e+00
eval[11] = -3.e+00 + i 0.e+00
eval[12] = 6.e+00 + i 0.e+00
eval[13] = -6.e+00 + i 0.e+00
eval[14] = 3.e+00 + i 0.e+00
eval[15] = -2.912706738243e-16 + i 0.e+00
eval[16] = 6.e+00 + i 0.e+00
eval[17] = -6.e+00 + i 0.e+00
eval[18] = -3.e+00 + i 0.e+00
---

Can you clarify which version of GSL you are using? Also you might want
to print out your matrix m to make sure it is initialized correctly. I
did not try to debug your code.

Best,
Patrick

On 12/29/20 7:44 AM, Dmitry Cheshkov wrote:
> Dear Sirs!
>
> I have to calculate the eigenvalues of the following 18x18 real
> non-symmetric square matrix:
>
> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
> 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0
> 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 0 0 0
> 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 3
> 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 6 0
> 0 0 0 0 0 0 0 0 0 6 0 3 0 0 0 6 0 3
> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
> 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0
> 0 0 0 0 0 0 0 0 0 0 0 0 4 0 2 0 0 0
> 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0 0 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0 3 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0
> 0 6 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0
> 6 0 3 0 0 0 6 0 3 0 0 0 0 0 0 0 0 0
> 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0 0 4 0 2 0 0 0 0 0 0 0 0 0 0 0 0
>
> I wrote the folowing code:
>
> #include 
> #include 
> #include 
> #include 
>
> using namespace std;
>
> int main(void)
> {
>
>   ifstream istr("matrix");
>   int N = 18;
>   gsl_matrix* m = gsl_matrix_alloc(N, N);
>   gsl_vector_complex* eval = gsl_vector_complex_alloc(N);
>   gsl_eigen_nonsymm_workspace* w = gsl_eigen_nonsymm_alloc(N);
>   for (int i = 0; i < N; i++)
>  for (int j = 0; j < N; j++)
>     istr >> m->data[N * j + i];
>   gsl_eigen_nonsymm(m, eval, w);
>   for (int i = 0; i < N; i++)
>    cout << eval->data[2 * i] << '\t' << eval->data[2 * i + 1] << endl;
>   cout << endl;
>   return 0;
>
> }
>
> If the file "matrix" contains the above mentioned matrix, I have the
> following error:
>
> gsl: francis.c:209: ERROR: maximum iterations reached without finding
> all eigenvalues
> Default GSL error handler invoked.
>
> While Wolfram Mathematica works well:
> In[3]:= Eigenvalues[matrix]
> Out[3]= {-6, -6, -6, -6, 6, 6, 6, 6, -3, -3, 3, 3, 0, 0, 0, 0, 0, 0}
>
> Is it possible to calculate this eigenvalues using GSL library and how?
>
>
> With best regards,
> Dmitry Cheshkov
>

#include 
#include 
#include 

#include 
#include 
#include 
#include 

int
main()
{
  const double data[] = {
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 1, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 3,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 6, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0, 0, 0, 6, 0, 3,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 3, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 6, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 6, 0, 3, 0, 0, 0, 6, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  const size_t N = 18;
  gsl_matrix_const_view m = gsl_matrix_const_view_array(data, N, N);
  gsl_matrix * A = gsl_matrix_alloc(N, N);
  gsl_vector_complex * eval = gsl_vector_complex_alloc(N);
  gsl_eigen_nonsymm_workspace * w = gsl_eigen_nonsymm_alloc(N);
  size_t i;

  gsl_matrix_memcpy(A, );
  gsl_eigen_nonsymm(A, eval, w);

  for (i = 0; i < N; ++i)
{
  gsl_complex ei = gsl_vector_complex_get(eval, i);

  printf("eval[%zu] = %.12e + i %.12e\n",
 i + 1, GSL_REAL(ei), GSL_IMAG(ei));
}

  gsl_matrix_free(A);
  

[bug #59759] spmatrix test fails on x86_64

[Patrick Alken]

Follow-up Comment #1, bug #59759 (project gsl):

from noloader =at= gmail =dot= com

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59759] spmatrix test fails on x86_64

[Patrick Alken]

URL:
  <https://savannah.gnu.org/bugs/?59759>

 Summary: spmatrix test fails on x86_64
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 23 Dec 2020 10:08:37 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

On Fri, Aug 28, 2020 at 9:47 PM Patrick Alken
 wrote:
>
> Hello, I just committed a fix to the git repository. Would you be able
> to test the git version to see if its resolved on your platform?

Hi Patrick.

Here's the result for the tip of master:

make  check-TESTS
make[2]: Entering directory '/home/jwalton/Build-Scripts/gsl-2.6/specfunc'
make[3]: Entering directory '/home/jwalton/Build-Scripts/gsl-2.6/specfunc'
FAIL: test

Testsuite summary for gsl 2.6+

# TOTAL: 1
# PASS:  0
# SKIP:  0
# XFAIL: 0
# FAIL:  1
# XPASS: 0
# ERROR: 0

See specfunc/test-suite.log


$ cat gsl-2.6/specfunc/test-suite.log
===
   gsl 2.6+: specfunc/test-suite.log
===

# TOTAL: 1
# PASS:  0
# SKIP:  0
# XFAIL: 0
# FAIL:  1
# XPASS: 0
# ERROR: 0

.. contents:: :depth: 2

FAIL: test
==

FAIL: res[60]+0.0 [2351]
  expected: 6.7753780053831863e+71
  obtained: 6.7753780053832010e+71
  fracdiff: 1.0856927230202910e-15
  value not within tolerance of expected value
  6.775378005383201023e+71
FAIL: res[100]+0.0 [2353]
  expected: 2.9579660004912027e+118
  obtained: 2.9579660004911946e+118
  fracdiff: 1.3629435007875587e-15
  value not within tolerance of expected value
  2.957966000491194605e+118
FAIL: gsl_sf_hermite_array_deriv(23, 100, 0.75) [2354]
FAIL: Hermite Functions [2498]
FAIL: gsl_sf_laguerre_n_e(1e6+1, 2.5, 2.5, ) [2991]
  expected: -4.8017961545391276e+05
  obtained: -4.8017961546181823e+05 ± 3.6891004271829821e-02
(rel=7.68275e-08)
  fracdiff: 8.2317878690294596e-12
 tolerance: 3.6379788070917130e-12
  value not within tolerance of expected value
  -4.801796154618182336e+05  3.689100427182982062e-02
FAIL: Laguerre Polynomials [3030]
FAIL test (exit status: 1)

==

On Fri, Aug 28, 2020 at 9:47 PM Patrick Alken
 wrote:
>
> Hello, I just committed a fix to the git repository. Would you be able
> to test the git version to see if it's resolved on your platform?

My bad... I should have built with 'make -k' to get back to the spmatrix
test...

Here is the result of master for spmatrix. Things have improved:

$ cat gsl-2.6/spmatrix/test-suite.log
===
   gsl 2.6+: spmatrix/test-suite.log
===

# TOTAL: 1
# PASS:  0
# SKIP:  0
# XFAIL: 0
# FAIL:  1
# XPASS: 0
# ERROR: 0

.. contents:: :depth: 2

FAIL: test
==

FAIL: gsl_spmatrix_complex_float_scale[53,107](COO) [21671]
FAIL: gsl_spmatrix_complex_float_scale[53,107](CSC) [28476]
FAIL: gsl_spmatrix_complex_float_scale[53,107](CSR) [35281]
FAIL: gsl_spmatrix_complex_float_scale[40,20](COO) [65775]
FAIL: gsl_spmatrix_complex_float_scale[40,20](CSC) [66416]
FAIL: gsl_spmatrix_complex_float_scale[40,20](CSR) [67057]
FAIL: gsl_spmatrix_complex_float_scale[30,30](COO) [76375]
FAIL: gsl_spmatrix_complex_float_scale[30,30](CSC) [78176]
FAIL: gsl_spmatrix_complex_float_scale[30,30](CSR) [79977]
FAIL test (exit status: 1)

I see the same results when I backport 76f5f97e3389 to the GSL 2.6 tarball.

Attached are the test-suite logs from a 'make -k' from master. All of
the test results should be present.

Jeff


==

Hi Everyone/Patrick,

A quick note since I'm making a testing pass on GSL...

Ubuntu 18.04 LTS x86_64 (fully patched) fails to link out of the box
due to a missing -lm. The symptom is undefined references for symbols
like sincos.

It may be a good idea to add an Autoconf test for:

* -lm, to avoid the link error on Linux
* -Wl,--as-needed to discard unneeded libraries

-lm goes in LIBS, and -Wl,--as-needed goes in AM_CFLAGS when available.

Autoconf will probably be the best in this case since -lm may be
missing on the target platform. I'm not sure how AIX, BSDs, OS X and
Solaris are going to act. I have not gotten that far in testing.

Jeff


==

Hi Everyone/Patrick,

A quick testing pass with CFLAGS that include '-fsanitize=undefined
-fno-sanitize-recover=all' shows some undefined behavior in the
library or test suite. Several of them are listed below.

Some of them, like "runtime error: left shift of 1 by 31

[bug #59758] Doc bug

[Patrick Alken]

URL:
  

 Summary: Doc bug
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 23 Dec 2020 10:01:14 PM UTC
Category: Documentation
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from alzhou10 =at= hotmail =dot= com

I think in https://www.gnu.org/software/gsl/doc/html/multimin.html, there is a
typo in the Fletcher-Reeves formula. It should be
p'=-g'-βp. 




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59710] Undefined behavior error in Nonlinear Least-Squares Fitting

[Patrick Alken]

Update of bug #59710 (project gsl):

 Open/Closed:Open => Closed 

___

Follow-up Comment #1:


  Thanks for your report. The problem is that on the first iteration,
the lambda value jumps to -12.8, and the value of exp(12.8 * 60) is too
large to represent in double precision, so infs and nans appear in the
residual vector and Jacobian matrix.

For this particular problem, the value of lambda should be > 0, however
there is no constraint placed to assure that, so when the first
iteration jumps to -12.8, bad numerical things happen.

I think the best solution for now is to adjust the model parameters a
little to prevent this type of behavior, and bring TMAX closer to 1
(instead of 40). A better solution would be to implement a constraint so
that lambda > 0, but this would require a lot of development effort.

Fixed in commit 7caac3346

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59710] Undefined behavior error in Nonlinear Least-Squares Fitting

[Patrick Alken]

URL:
  

 Summary: Undefined behavior error in Nonlinear Least-Squares
Fitting
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Thu 17 Dec 2020 04:17:55 PM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from zhoulai =dot= fu =at= gmail =dot= com

Running the attached program nlfit_bug.c triggers undefined behavior,
calculating unreasonable
results. The program "nlfit_debug.c" is adapted from the example program
provided
by GSL's doc/examples/nlfit.c. Their difference lies in only one parameter
40.0
versus 60.0:

1d0
<
12c11,12
< #define TMAX   (60.0) /* time variable in [0,TMAX] */ //40.0 instead of
60.0 in nlfit.c
---
> //#define TMAX   (40.0) /* time variable in [0,TMAX] */
> #define TMAX   (60.0) /* time variable in [0,TMAX] */


The undefined behavior can be better visualized if GCC's Undefined
Behavior Sanitizer is enabled (-fsanitize=undefined). Below is the
execution results:

data: 59.3939 1.03153 0.101317
data: 60 1.14181 0.101239
iter  0: A = 1., lambda = 1., b = 0., cond(J) =  inf,
|f(x)| = 101.0200
iter  1: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  92.8216,
|f(x)| = -nan
iter  2: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter  3: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter  4: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan

*nielsen.c:98:7: runtime error: left shift of 4611686018427387904 by 1
places cannot be represented in type 'long int'*iter  5: A = 3.5110, lambda
= -12.8820, b = 1.2364, cond(J) =  nan, |f(x)| = -nan
iter  6: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter  7: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter  8: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter  9: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 10: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 11: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 12: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 13: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 14: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 15: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 16: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 17: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 18: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 19: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 20: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 21: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 22: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 23: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 24: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 25: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 26: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 27: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 28: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 29: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 30: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 31: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 32: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 33: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 34: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 35: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 36: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 37: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 38: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 39: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 40: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
|f(x)| = -nan
iter 41: A = 

Re: Undefined behavior error in Nonlinear Least-Squares Fitting

[Patrick Alken]

Hello,

  Thanks for your report. The problem is that on the first iteration,
the lambda value jumps to -12.8, and the value of exp(12.8 * 60) is too
large to represent in double precision, so infs and nans appear in the
residual vector and Jacobian matrix.

For this particular problem, the value of lambda should be > 0, however
there is no constraint placed to assure that, so when the first
iteration jumps to -12.8, bad numerical things happen.

I think the best solution for now is to adjust the model parameters a
little to prevent this type of behavior, and bring TMAX closer to 1
(instead of 40). A better solution would be to implement a constraint so
that lambda > 0, but this would require a lot of development effort.

I will modify the test program to make it less likely to run into this
problem in the future.

Thanks,
Patrick

On 12/7/20 3:18 PM, Zhoulai Fu@Gmail wrote:
> Running the attached program nlfit_bug.c triggers undefined behavior,
> calculating unreasonable
> results. The program "nlfit_debug.c" is adapted from the example program
> provided
> by GSL's doc/examples/nlfit.c. Their difference lies in only one parameter
> 40.0
> versus 60.0:
>
> 1d0
> <
> 12c11,12
> < #define TMAX   (60.0) /* time variable in [0,TMAX] */ //40.0 instead of
> 60.0 in nlfit.c
> ---
>> //#define TMAX   (40.0) /* time variable in [0,TMAX] */
>> #define TMAX   (60.0) /* time variable in [0,TMAX] */
>
> The undefined behavior can be better visualized if GCC's Undefined
> Behavior Sanitizer is enabled (-fsanitize=undefined). Below is the
> execution results:
>
> data: 59.3939 1.03153 0.101317
> data: 60 1.14181 0.101239
> iter  0: A = 1., lambda = 1., b = 0., cond(J) =  inf,
> |f(x)| = 101.0200
> iter  1: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  92.8216,
> |f(x)| = -nan
> iter  2: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter  3: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter  4: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
>
> *nielsen.c:98:7: runtime error: left shift of 4611686018427387904 by 1
> places cannot be represented in type 'long int'*iter  5: A = 3.5110, lambda
> = -12.8820, b = 1.2364, cond(J) =  nan, |f(x)| = -nan
> iter  6: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter  7: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter  8: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter  9: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 10: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 11: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 12: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 13: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 14: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 15: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 16: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 17: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 18: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 19: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 20: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 21: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 22: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 23: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 24: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 25: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 26: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 27: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 28: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 29: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 30: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 31: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 32: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 33: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 34: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = -nan
> iter 35: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) =  nan,
> |f(x)| = 

[bug #59624] Buffer overflow in gsl_stats_quantile_from_sorted_data

[Patrick Alken]

URL:
  

 Summary: Buffer overflow in
gsl_stats_quantile_from_sorted_data
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Fri 04 Dec 2020 10:06:37 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from zhoulai.fu =at= gmail =dot= com

Running the following code (also attached as a file) triggers a segmentation
error.

#include 
#include 
#include 

int main(void)
{
  double upperq;
  double data[5] = {17.2, 18.1, 16.5, 18.3, 12.6};
  gsl_sort (data, 1, 5);
  upperq = gsl_stats_quantile_from_sorted_data (data, 1, 5, 675);
  return 0;
}
// gcc statsort_bug.c -lgsl -lgslcblas; ./a.out

The error points to statistics/quantiles_source.c:41:

  result = (1 - delta) * sorted_data[lhs * stride] + delta *
sorted_data[(lhs + 1) * stride] ;

The segmentation error is due to a stack buffer overflow (where
lhs=2700, strid=1 as shown in GDB). The bug could be exploited for
security attack, knowing that it occurs when the quantile "f" is
beyond the expected [0,1] range (f=675 in this case).



___

File Attachments:


---
Date: Fri 04 Dec 2020 10:06:37 PM UTC  Name: statsort_bug.c  Size: 316B   By:
psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #59494] Bug in gsl_rstat

[Patrick Alken]

URL:
  

 Summary: Bug in gsl_rstat
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Fri 20 Nov 2020 04:21:29 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

Bug in gsl_rstat_median for small number of samples (code example attached)



___

File Attachments:


---
Date: Fri 20 Nov 2020 04:21:29 PM UTC  Name: bug.cpp  Size: 1KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: bug in gsl_rstat_workspace

[Patrick Alken]

Thanks, I've logged the issue into the bug tracker, and will investigate
when I find time

On 11/20/20 3:40 AM, Wilhelm Braunschober wrote:
> Yours Wilhelm Braunschober

Re: spmatrix tests fail on x86_64

[Patrick Alken]

Hello, I just committed a fix to the git repository. Would you be able
to test the git version to see if its resolved on your platform?

Patrick

On 8/28/20 6:26 PM, Jeffrey Walton wrote:
> Hi Everyone,
>
> I'm building the GSL 2.6 release tarball on Ubuntu 18.04 LTS x86_64
> fully patched. 'make check' is failing at the spmatrix test suite:
>
> $ cat gsl-2.6/spmatrix/test-suite.log
> ==
>gsl 2.6: spmatrix/test-suite.log
> ==
>
> # TOTAL: 1
> # PASS:  0
> # SKIP:  0
> # XFAIL: 0
> # FAIL:  1
> # XPASS: 0
> # ERROR: 0
>
> .. contents:: :depth: 2
>
> FAIL: test
> ==
>
> FAIL: gsl_spmatrix_complex_scale_columns[53,107](COO) [1132]
> FAIL: gsl_spmatrix_complex_scale_rows[53,107](COO) [1133]
> FAIL: gsl_spmatrix_complex_scale_columns[53,107](CSC) [1135]
> FAIL: gsl_spmatrix_complex_scale_rows[53,107](CSC) [1136]
> FAIL: gsl_spmatrix_complex_scale_columns[53,107](CSR) [1138]
> FAIL: gsl_spmatrix_complex_scale_rows[53,107](CSR) [1139]
> FAIL: gsl_spmatrix_complex_float_scale[53,107](COO) [1211]
> FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](COO) [1212]
> FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](COO) [1213]
> FAIL: gsl_spmatrix_complex_float_scale[53,107](CSC) [1214]
> FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](CSC) [1215]
> FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](CSC) [1216]
> FAIL: gsl_spmatrix_complex_float_scale[53,107](CSR) [1217]
> FAIL: gsl_spmatrix_complex_float_scale_columns[53,107](CSR) [1218]
> FAIL: gsl_spmatrix_complex_float_scale_rows[53,107](CSR) [1219]
> FAIL: gsl_spmatrix_complex_scale_columns[40,20](COO) [2450]
> FAIL: gsl_spmatrix_complex_scale_rows[40,20](COO) [2451]
> FAIL: gsl_spmatrix_complex_scale_columns[40,20](CSC) [2453]
> FAIL: gsl_spmatrix_complex_scale_rows[40,20](CSC) [2454]
> FAIL: gsl_spmatrix_complex_scale_columns[40,20](CSR) [2456]
> FAIL: gsl_spmatrix_complex_scale_rows[40,20](CSR) [2457]
> FAIL: gsl_spmatrix_complex_float_scale[40,20](COO) [2529]
> FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](COO) [2530]
> FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](COO) [2531]
> FAIL: gsl_spmatrix_complex_float_scale[40,20](CSC) [2532]
> FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](CSC) [2533]
> FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](CSC) [2534]
> FAIL: gsl_spmatrix_complex_float_scale[40,20](CSR) [2535]
> FAIL: gsl_spmatrix_complex_float_scale_columns[40,20](CSR) [2536]
> FAIL: gsl_spmatrix_complex_float_scale_rows[40,20](CSR) [2537]
> FAIL: gsl_spmatrix_complex_scale_columns[30,30](COO) [3768]
> FAIL: gsl_spmatrix_complex_scale_rows[30,30](COO) [3769]
> FAIL: gsl_spmatrix_complex_scale_columns[30,30](CSC) [3771]
> FAIL: gsl_spmatrix_complex_scale_rows[30,30](CSC) [3772]
> FAIL: gsl_spmatrix_complex_scale_columns[30,30](CSR) [3774]
> FAIL: gsl_spmatrix_complex_scale_rows[30,30](CSR) [3775]
> FAIL: gsl_spmatrix_complex_float_scale[30,30](COO) [3847]
> FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](COO) [3848]
> FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](COO) [3849]
> FAIL: gsl_spmatrix_complex_float_scale[30,30](CSC) [3850]
> FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](CSC) [3851]
> FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](CSC) [3852]
> FAIL: gsl_spmatrix_complex_float_scale[30,30](CSR) [3853]
> FAIL: gsl_spmatrix_complex_float_scale_columns[30,30](CSR) [3854]
> FAIL: gsl_spmatrix_complex_float_scale_rows[30,30](CSR) [3855]
>
> It looks like the issue was raised for Aarch64 and PowerPC at
> https://mail.gnu.org/archive/html/bug-gsl/2020-05/msg1.html. But I
> don't see a resolution.
>
> Attached are config.log and the test suite results.
>
> Jeff

[bug #58977] Any plans to make new release?

[Patrick Alken]

Update of bug #58977 (project gsl):

 Open/Closed:Open => Closed 


___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: Observed large error from spherical Bessel functions at high orders

[Patrick Alken]

Hello,

  Unfortunately some of the special functions do not produce the
expected accuracy for all inputs, and it seems you've found another
example. It would be great if you could look into algorithms for
computing these functions and possibly find a way to improve the
accuracy for these inputs. Unfortunately we don't have a lot of
developers working on GSL these days, and I don't have much expertise in
computing the bessel functions.

Please at least log this issue in the bug tracker:
http://savannah.gnu.org/bugs/?group=gsl

That way it won't be lost and someone may eventually look into it.

Thanks,
Patrick

On 7/14/20 11:05 AM, Ziqi Fan wrote:
> Dear developers and maintainers of GNU Scientific Library,
>
> My name is Ziqi Fan. I am a PhD candidate from University of Florida. I
> have been using GSL for at least 3 years for my numerical solver project.
> Recently, I met a bug in my solver and the bug should not exist from a
> theoretical perspective. So I checked very carefully my own implementation
> and found that the bug originated from using the routine for spherical
> bessel functions of the second kind: gsl_sf_bessel_yl.
>
> For instance, I tested the function using a large order n = 45, and got the
> following result: y_45(6.975948) = -726330209582507479571265224704.00,
> err = 19103761820604644.00, where the first term is the output and the
> second term is an estimated error provided by the routine
> "gsl_sf_bessel_yl_e". The error is relatively small compared to the output,
> but it can easily cause divergence of a numerical algorithm, considering
> the magnitude of its value.
>
> I once implemented spherical bessel functions myself using the numerical
> recipe, and I understand that large error occurs when x << n, where x is a
> positive input to the function, and n is the order of the function. I am
> wondering if it is possible to resolve the large absolute error at high
> order of the spherical bessel functions of the second kind. If I have an
> opportunity to communicate with an engineer or mathematician from you, I
> may be able to help contribute a routine with a higher precision. If so, a
> new routine with a higher precision can really contribute to the invention
> of many new algorithms for powerful numerical solvers.
>
> I really appreciate your effort in providing and maintaining this great
> software, and I look forward to hearing from you!
>
> Best,
> Ziqi

Re: Extrapolation returns NaN

[Patrick Alken]

On 6/19/20 6:13 AM, Yury Pakhomov wrote:
> Hello All!
>
> gsl_spline_eval returns NaN values when x is out of range.
>
> Since GSL-1.15 this situation remains unresolved
> http://git.savannah.gnu.org/cgit/gsl.git/tree/NEWS#n470
>
> BW, Yuri
>
This is the documented behavior ((see gsl_interp_eval description).

It would be nice to add extrapolation, but it has not been done yet.

[bug #52540] Complex QR decomposition

[Patrick Alken]

Update of bug #52540 (project gsl):

  Status:None => Fixed  


___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[bug #52540] Complex QR decomposition

[Patrick Alken]

Update of bug #52540 (project gsl):

 Open/Closed:Open => Closed 

___

Follow-up Comment #2:

Complex QR has been added to GSL

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: 'syrk' very slow compared to 'gemm'

[Patrick Alken]

Hello,

  The GSL CBLAS library is solely intended so that users can compile 
applications "out of the box" without needing to install an optimized 
BLAS library. The GSL CBLAS is not intended for serious work, and should 
never be used for performance critical applications. We always recommend 
installing an optimized BLAS library (such as ATLAS or MKL) and then 
never using gslcblas again.


Writing fast BLAS routines is non-trivial, and libraries such as ATLAS 
and MKL have spent many years (possibly decades) optimizing every single 
subroutine. Many of their routines are written in assembly language. The 
GSL project has no plans to try to reproduce this type of effort in the 
internal cblas library - as stated before our cblas library is only 
provided so that users can quickly install and start using GSL without 
installing additional dependencies. It should never be used for serious 
computations.


Hope this helps,
Patrick

On 2/19/20 12:47 AM, David Cortes wrote:

I'm comparing functions 'ssyrk' and 'dsyrk' in the CBLAS section
against 'sgemm' and 'dgemm' for the following operation:
C <- t(A)*A

And I'm finding that 'syrk' is seriously underperformant compared to
'gemm', to the point that it seems as if the program had crashed
instead. In particular, I've timed it with input sizes of A of
1,000,000x100 and 1,000x10,000, on an AMD Ryzen 7 2700 processor, with
these results:

1,000,000 x 100 sgemm: 1.63 s
1,000,000 x 100 ssyrk: 44.2 s  <- big problem here
1,000,000 x 100 dgemm: 3.26 s
1,000,000 x 100 dsyrk: 41.5 s  <- big problem here
1,000 x 10,000 sgemm: 48.3 s
1,000 x 10,000 ssyrk: 68   s   <- smaller problem
1,000,000 x 100 dgemm: 95  s
1,000,000 x 100 dsyrk: 88  s   <- smaller problem
(A naive 3-stage for loop over the *full* array calculating things
twice, with -O3 optimization, takes about the same as syrk)

Compare that against MKL running single threaded:
1,000,000 x 100 sgemm: 406 ms
1,000,000 x 100 ssyrk: 269 ms
1,000,000 x 100 dgemm: 864 ms
1,000,000 x 100 dsyrk: 532 ms
1,000 x 10,000 sgemm: 3.59 s
1,000 x 10,000 ssyrk: 1.84 s
1,000 x 10,000 dgemm: 7.61 s
1,000 x 10,000 dsyrk: 3.84 s

And against OpenBLAS:
1,000,000 x 100 sgemm: 449 ms
1,000,000 x 100 ssyrk: 326 ms
1,000,000 x 100 dgemm: 935 ms
1,000,000 x 100 dsyrk: 653 ms
1,000 x 10,000 sgemm: 4.32 s
1,000 x 10,000 ssyrk: 2.04 s
1,000 x 10,000 dgemm: 8.41 s
1,000 x 10,000 dsyrk: 3.96 s

Now, 'syrk' being slightly slower than 'gemm' is something that could
perhaps be expected, but being 25x slower (and 160x slower than MKL) is
too much I think.

Setup: GSL 2.5, installed from debian repositories.

Below are the function calls (A being an n-by-k matrix)

#include "cblas.h"
//#include "mkl.h"
#include 
void call_gemm(float *A, float *C, int n, int k, float alpha)
{
 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
 k, k, n,
 alpha, A, k, A, k,
 0., C, k);
}
void call_syrk(float *A, float *C, int n, int k, float alpha)
{
 cblas_ssyrk(CblasRowMajor, CblasUpper, CblasTrans,
 k, n, alpha,
 A, k, 0., C, k);
}
void for_loop(float *restrict A, float *restrict C, int n, int k, float
alpha)
{
 memset(A, 0, n*k*sizeof(float));
 for (int row = 0; row < k; row++)
 for (int col = 0; col < k; col++)
 for (int dim = 0; dim < n; dim++)
 C[col + row*k] += alpha*A[row + dim*n];
}

Re: Inverse Gamma Issue

[Patrick Alken]

Thanks Robert,

   There have been many bug reports sent in recently - I need to add 
them all the tracker. Hopefully I can find time to look at your work, 
though it may take a while.

Thanks,
Patrick

On 1/28/20 4:09 AM, Janes, Robert (Columbus) via Bug reports for the GNU 
Scientific Library wrote:
> I sent an e-mail a couple weeks ago about a fix I came up with for a problem 
> with the inverse Chi-squared function.  The problem was in the inverse gamma 
> function and I suggested increasing the iteration limit to 50.  That fixed 
> the immediate problem but then we found another problem case which this did 
> not fix.
>
> The IMSL documentation for this stated that for large degrees of freedom 
> after 100 iterations the code returns the present iteration value as the best 
> answer it can supply.  I implemented this in the gamminv.c code and found 
> excellent agreement with IMSL when running side-by-side comparisons for 
> degrees of freedom from 1,000 to 32,000 in steps of 1,000.  You might 
> consider this change over the one I suggested earlier.
>
> See the attached, updated source code.
>
>
> Rob Janes
>
> SGS Transportation
>
> Senior Software Engineer
>
> SGS – CMX
>
> 2860 N. National Road Suite A
>
> US – 47201 – Columbus, Indiana
>
> Phone: +1 - 812 - 378 - 7966
>
> Fax: +1 812 - 378 - 3393
>
> Email: robert.ja...@sgs.com
>
>
>
> Information in this email and any attachments is confidential and intended 
> solely for the use of the individual(s) to whom it is addressed or otherwise 
> directed. Please note that any views or opinions presented in this email are 
> solely those of the author and do not necessarily represent those of the 
> Company. Finally, the recipient should check this email and any attachments 
> for the presence of viruses. The Company accepts no liability for any damage 
> caused by any virus transmitted by this email. All SGS services are rendered 
> in accordance with the applicable SGS conditions of service available on 
> request and accessible at https://www.sgs.com/en/terms-and-conditions

Re: [bug #57173] Feature request: zeta function for complex arguments

[Patrick Alken]

This brings up an issue which I've wondered about for a while. GSL 
currently conforms to the c89 standard. I believe the previous 
maintainer (Brian) wanted this so that users on older/legacy systems 
could continue to use GSL. I have continued this tradition, and for 
every release, I make sure GSL will compile with -std=c89.

If we begin using features of later C standards (c99, c11, etc), then 
GSL will no longer compile with -std=c89.

I guess the question is whether we want to continue to support c89, or 
abandon it in favor of newer features of C.

Patrick

On 11/5/19 1:20 PM, Gerard Jungman wrote:
> On 11/5/19 11:26 AM, Mark Galassi wrote:
>>>>>>> "Patrick" == Patrick Alken  writes:
>>
>>
>>  Patrick> I don't think any progress was made on this
>>
>> One update to the discussion from back then is that C++11 has formalized
>> data compatibilty between the C11 _Complex type and the C++ complex
>> class.  This also raises the architectural question of whether GSL might
>> want to use C's complex type.
>
> Although, strangely enough, they made support for
> complex types optional. Does anybody know the
> rationale for this decision? This kind of
> back-peddling seems very strange.
>
> I have been using the new types myself in C code,
> for several years. Of course, gcc and clang are not
> going to drop support for the types. But I worry that
> this "optionality" means they may not get the love
> and attention they need. Or worse, that something
> weird might happen leading to a divergence of
> implementations.
>
>
> It would be good to move to the standard, if we know
> it is not going to deteriorate in some way, and that
> there are no hidden issues with the implementations
> as they exist now. There are some weird hidden
> issues, like the meaning and consequences of
> the CX_LIMITED_RANGE pragma.
>
> At the very least, it would be good to read the normative
> parts carefully. Annex G of either the c17 draft standard
> (n2176) or the c11 draft standard (n1570). The two seem
> to be identical, but I have not taken out a microscope.
>

Re: [bug #57173] Feature request: zeta function for complex arguments

[Patrick Alken]

I don't think any progress was made on this

On 11/5/19 9:49 AM, Mark Shoulson wrote:
> URL:
>
>
>   Summary: Feature request: zeta function for complex 
> arguments
>   Project: GNU Scientific Library
>  Submitted by: clsn
>  Submitted on: Tue 05 Nov 2019 04:49:32 PM UTC
>  Category: Accuracy problem
>  Severity: 3 - Normal
>  Operating System:
>Status: None
>   Assigned to: None
>   Open/Closed: Open
>   Release:
>   Discussion Lock: Any
>
>  ___
>
> Details:
>
> It's a little surprising that there is no support in GSL for the Riemann zeta
> functions (and related functions) for complex arguments.  This was apparently
> discussed back in 2008; web searching found me a thread at
> https://www.mail-archive.com/help-gsl@gnu.org/msg02130.html; was anything done
> with this?
>
> Thanks.
>
>
>
>
>  ___
>
> Reply to this item at:
>
>
>
> ___
>Message sent via Savannah
>https://savannah.gnu.org/
>
>

Re: [Bug-gsl] Possible error handling/return bug in vector/copy_source.c

[Patrick Alken]

Yes you're right, thanks for catching this. I put the vector length
checks first now on the git. This code was just recently changed for
v2.6 to add the BLAS calls, so I overlooked that.

Patrick

On 9/22/19 12:55 PM, Diego via Bug-gsl wrote:
> I've been playing with a wrapper for GSL for Chicken Scheme and I've found 
> what I believe might be a bug in the file vector/copy_source - do let me know 
> if I've missed something, because I'm extremely new to this code base. In 
> FUNCTION(gsl_vector, memcpy), the logic that checks that the length of the 
> vectors is equal falls under the #else clause.
>
> This means that (especially with no error handler set, as is recommended for 
> production code) if one of the other #if clauses is triggered during 
> compilation, memcpy will call a blas copy function instead, and can (and 
> will) fail silently - it calls the blas version /without/ returning, and then 
> goes on to unconditionally return GSL_SUCCESS. (In contrast, the memcpy in 
> matrix/copy_source.c checks the lengths before going into the #if/blas logic).
>
> For the proper errno codes to be returned, I would propose that either: a) 
> the vector length check and GSL_ERROR call should be before the first #if, as 
> in matrix/copy_source.c or b) each clause that ends up calling a blas 
> function should be using e.g. `return gsl_blas_dcopy(...)` so that the proper 
> error codes are returned. I think option a would be preferable for 
> consistency and so as to not have to rely on the internals of the gsl_blas 
> API.
>
> - Diego
>
>
>

Re: [Bug-gsl] About the web page of GSL - GNU Scientific Library

[Patrick Alken]

Thanks, I have updated the link

On 8/20/19 9:20 PM, mrkmh...@tmu.ac.jp wrote:
> In the web page of 
>GSL - GNU Scientific Library
>   URL 
>
> You write
>
>> A Japanese translation is also available online (may not be the most 
> recent version).
>>GSL Reference Manual - Japanese Translation (by Daisuke Tominaga, 
> AIST Computational Biology Research Center)
>
>
> However, the linked URL 
>
> http://www.cbrc.jp/~tominaga/translations/gsl/gsl-1.8/index.html
>
> on your Web page seems not working (You can try).
>
> But, I suppose it is currently changed to the following URL, I guess
>
>   http://cbrc3.cbrc.jp/~tominaga/translations/gsl/gsl-1.15/index.html
>
>
>
> ---
> Hiroshi Murakami
> Department of Mathematical Sciences
> Tokyo Metropolitan University
> Email: mrkmh...@tmu.ac.jp
> 
>
>
>

Re: [Bug-gsl] the multifit_nlinear subpackage is mucked up

[Patrick Alken]

On 8/8/19 3:30 AM, JohnT wrote:
> It fails to test ok.
>
>
>
>
Can you give details, or post the output from your test*.log files?

Re: [Bug-gsl] Hang caused by Inf inputs in hermite special function

[Patrick Alken]

Hello, some work was done on the hermite functions in the current branch
of GSL on the git repository. Could you test this issue by checking out
the latest git version?

Thanks,
Patrick

On 7/15/19 7:20 PM, Jackson Vanover wrote:
> Using version 2.5 of GSL from ftp://ftp.gnu.org/gnu/gsl/gsl-2.5.tar.gz
>
> OS is Ubuntu 18.04.2 LTS
>
> Hardware is a Dell XPS13 with Intel Corporation Xeon E3-1200 v6/7th Gen
> Core Processor Host Bridge/DRAM Registers [8086:5914] (rev 08)
>
> Compiler is gcc 7.4.0. Compulation options shown below
> gcc -O0 -I/path/to/include hangTest.c -o hangTest -L/path/to/lib -lgsl
> -lgslcblas -lm
>
> Description: When feeding gsl_sf_hermite_phys a positive or negative
> infinity, it hangs. The while condition at line 366 of hermite.c never
> evaluates to false.
>
> #include 
> #include 
>
> void main(int argc, char * argv){
> gsl_sf_hermite_phys(2,GSL_POSINF);
> }
>
> Best,
> Jackson Vanover

Re: [Bug-gsl] Building DLLs on Windows

[Patrick Alken]

Philip,

   I forwarded your email to Brian (cc'd) who has done a lot of work on windows 
builds of GSL. His response is pasted below.


To build GSL as it is now requires that WIN32 is defined somewhere in
the GSL code.

The _WIN32 define is always defined by the MSVC compiler so it this were
used it would mean that we did not have to define WIN32 oursleves.

But there are other (i.e. non Microsoft) compilers and it is not clear
to me that they all define _WIN32 so it seems better not to rely on this
and simply use and define WIN32 in the GSL source code.

I didn't put the WIN32 define in the code so I don't know if they had
the same rationale but I don't think it is a mistake since all it means
is that we define it within the GSL build if appropriate.

I am not subscribed to the GSL lists so I can't currently post there but
you could point him at my Githib repository at:

https://github.com/BrianGladman/gsl

where there is a fully working GSL DLL build (and where I am willing to
help out if necessary).

best regards,

  Brian



On 6/19/19 1:59 AM, P wrote:
> Hi there,
>
> I believe there is a slight typo in gsl_types.h
>
> Where it should be "#ifdef _WIN32"
> Source
> https://docs.microsoft.com/en-us/previous-versions/visualstudio/visual-studio-2008/b0084kay(v=vs.90)
>
> Without this change DLLs do not build.
>
> Secondly, it seems a large number of functions are missing the prefix
> "GSL_VARS" and cannot be used in DLLs built by MSVC.
>
> More information on this can be found here:
> https://github.com/ampl/gsl/issues/24#issuecomment-503330162
>
> Kind regards
> Philip Deegan
> https://github.com/Dekken

Re: [Bug-gsl] Minor inaccuracy in comment

[Patrick Alken]

Thank you, I have corrected this on the git

On 2/19/19 1:46 AM, Chris Hiszpanski wrote:
> Hello,
>
> Noticed a comment in rng/taus113.c is inaccurate -- z_1 is computed from z_0, 
> not the other way around. Patch is provided below.
>
> Chris
>
> diff --git a/rng/taus113.c b/rng/taus113.c
> index 04da6a89..a46248bc 100644
> --- a/rng/taus113.c
> +++ b/rng/taus113.c
> @@ -39,7 +39,7 @@
>      work on architectures where integers are 64-bit.
>
>      The generator is initialized with
> -   zi = (69069 * z{i+1}) MOD 2^32 where z0 is the seed provided
> +   z{i+1} = (69069 * z{i}) MOD 2^32 where z0 is the seed provided
>      During initialization a check is done to make sure that the initial seeds
>      have a required number of their most significant bits set.
>      After this, the state is passed through the RNG 10 times to ensure the

Re: [Bug-gsl] gsl 2.3.0 bug?

[Patrick Alken]

Hello,

   I converted your main.cpp to C (I couldn't make it compile with my 
g++). Also I linked it against GSL v2.5, but the CQUAD code hasn't 
changed since 2.3 as far as I know. The program works as expected (C 
code attached).

$ ./main
result = 4.390157623791e-164
abserr = 0.e+00
nevals = 33

Mathematica agrees with the result to about 11 significant digits.

I am unable to diagnose problems with fgsl since I have never used it.

Patrick

On 11/20/18 9:36 AM, Rico Buchholz wrote:
>
> Hello!
>
> I am working with a fortran project which uses fgsl (v1.2.0, together 
> with gsl v2.3.0).
>
> The code reports an arithmetic exception in an integration routine, 
> when trying to integrate exp(-x^2) from x=19.3 to 20.0.
> I would expect the code to work without problems (maybe underflow or 
> denormal).
> Trying to create a minimal working example i came up with the attached 
> file (the cmake input i used for compilation is also attached, sorry I 
> could not boil this further down).
> As a corresponding version for c++ works (also attached), this might 
> relate to a bug in fgsl, but for me (a person without knowledge of the 
> internals), that seems questionable according to what causes the 
> problem and what not.
>
>
> The output is as follows:
>
> Program received signal SIGFPE: Floating-point exception - erroneous 
> arithmetic operation.
>
> Backtrace for this error:
> #0  0x7fa39acc5d1d in ???
> #1  0x7fa39acc4f7d in ???
> #2  0x7fa39a1e405f in ???
> #3  0x560ac54c9eff in ???
> #4  0x560ac54860cf in neo2
>     at 
> /afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:31
> #5  0x560ac54861d8 in main
>     at 
> /afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:4
>
> It does not occur with the alternative integration routine commented 
> in the code.
> This does not occur, if one reduces the lower bound to 19.2 or smaller 
> (which is a reason why I suspect that this might be a gsl bug - a bug 
> in the interface should have no influence on this).
>
>
> With kind regards
> Rico Buchholz
>
>
> PS: I did also send the bug report to fgsl.
>

/**
 *
 * g++ main.cpp -o main.x -L`gsl-config --libs`
 *
 */

#include 
#include 
#include 

#include 
#include 

double f(double x, void * params) {
  return exp(-x*x);
}

int main (void)
{
  size_t const n = 100;
  double const a = 19.3;
  double const b = 20.0;
  double const epsabs = 1.0e-13;
  double const epsrel = 1.0e-13;
  double result = 0.0;
  double abserr = 0.0;
  size_t nevals = 0;

  gsl_integration_cquad_workspace * workspace;
  gsl_function F;

  F.function = 
  F.params = NULL;

  workspace = gsl_integration_cquad_workspace_alloc(n);

  gsl_integration_cquad(, a, b, epsabs, epsrel, workspace, , , );

  fprintf(stderr, "result = %.12e\n", result);
  fprintf(stderr, "abserr = %.12e\n", abserr);
  fprintf(stderr, "nevals = %zu\n", nevals);

  gsl_integration_cquad_workspace_free(workspace);

  return 0;
}

[Bug-gsl] [bug #54925] Reduce cache misses for source_gemm_r

[Patrick Alken]

URL:
  

 Summary: Reduce cache misses for source_gemm_r
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 31 Oct 2018 01:57:19 PM UTC
Category: Performance
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

>From max =dot= bruce12 =at= gmail =dot= com

Hey guys,

I noticed that your matrix multiplication code had bad cache
performance due to a misordering of a loop. In a replicated version of
my change, I saw about 20% performance gains on my AMD FX CPU.

-Max



___

File Attachments:


---
Date: Wed 31 Oct 2018 01:57:19 PM UTC  Name:
0001-Reduce-cache-misses-for-source_gemm_r.patch  Size: 2KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[Bug-gsl] [bug #54140] Updating email

[Patrick Alken]

Update of bug #54140 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

fixed in commit 358d77e71fcd93248d6fae7c6997d1bd0dcfb1f8

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[Bug-gsl] [bug #54921] delete keyword clash in gsl_movstat.h

[Patrick Alken]

Update of bug #54921 (project gsl):

  Status:None => Fixed  
 Open/Closed:Open => Closed 

___

Follow-up Comment #1:

Fixed in commit d9a747a9d2819899e436bb546edd61bfcab0558e

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: [Bug-gsl] bug with gsl_integration

[Patrick Alken]

Thank you for finding a solution Johannes. I've never worked with
mex/GSL myself.

If the default handler is used, then the program will abort/segfault
anyway, so I don't think it should matter if you are calling from mex or
from a command line.

On 10/18/2018 11:26 PM, Johannes Hendriks wrote:
> Dear sir/madam,
>
> I would like to close this thread as it has been resolved. The problem
> ended up being the combination of the default error handler along with the
> mex compiler. The default error handler would do use 'printf' to output
> error message to screen and abort the function. Unfortunately this would
> cause matlab to crash as output to display requires mexPrintf when using
> mex compiler.
>
> So I recognise this is not a bug in your functions. However, maybe it would
> be good to highlight somewhere that if mex is being used then the default
> error handler needs to be disabled / replaced. I am sure there are quite a
> few academic and research users who might come across this problem.
>
> Regards,
> Johannes
>
> On Wed, Oct 17, 2018 at 1:18 PM Johannes Hendriks <
> johannes.n.hendr...@gmail.com> wrote:
>
>> Dear sir/madam
>>
>> I have found a possible issue with the gsl_integration_qng
>> and gsl_integration_qag functions. Under particular conditions when these
>> functions fail to meet the specified relative or absolute tolerance,
>> instead of failing gracefully they are causing a seg fault.
>>
>> I am running on osx 10.12.6 and compiling using the matlab mex compiling.
>> I have matlab 2017b, and am compiling with the command
>>
>> mex intTwoK_mex.c -I/usr/local/include/ -lmwlapack -lmwblas -lgsl
>>
>> I have attached a zip file containing a matlab script that will cause the
>> integration routine to fault.
>>
>> Let me know if you require any further infomation.
>>
>> Regards,
>> Johannes Hendriks
>>
>>
>>

Re: [Bug-gsl] [bug #54140] Updating email

[Patrick Alken]

Hello, it is done

On 06/18/2018 03:22 AM, Nicolas Darnis wrote:
> URL:
>   
>
>  Summary: Updating email
>  Project: GNU Scientific Library
> Submitted by: ndarnis
> Submitted on: Mon 18 Jun 2018 11:22:21 AM CEST
> Category: None
> Severity: 3 - Normal
> Operating System: 
>   Status: None
>  Assigned to: None
>  Open/Closed: Open
>  Release: 
>  Discussion Lock: Any
>
> ___
>
> Details:
>
> Hello, could you please update my email address in AUTHORS file please (and
> everywhere its necessary)? You must write 'ndar...@free.fr' instead of
> 'ndar...@cvf.fr'
> Thank you very much.
> Best regards,
> N. Darnis.
>
>
>
>
> ___
>
> Reply to this item at:
>
>   
>
> ___
>   Message sent via Savannah
>   https://savannah.gnu.org/
>
>

Re: [Bug-gsl] Small documentation bug

[Patrick Alken]

Hello, thanks for spotting this. I've fixed it on the git

On 06/14/2018 08:40 AM, Damien Desfontaines wrote:

The diagram for gsl_ran_binomial_pdf looks wrong:
 
https://www.gnu.org/software/gsl/doc/html/randist.html#c.gsl_ran_binomial_pdf
It looks like this is a binomial distribution with n=10, not with n=9.

Cheers,

Damien

[Bug-gsl] [bug #53919] handle large values correctly in (log_)erf(c) functions

[Patrick Alken]

Follow-up Comment #1, bug #53919 (project gsl):

The patch uses HUGE_VAL (C99) and isnan() checks - there is probably a cleaner
way to handle the problem, but will require further investigation

___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[Bug-gsl] [bug #53919] handle large values correctly in (log_)erf(c) functions

[Patrick Alken]

URL:
  

 Summary: handle large values correctly in (log_)erf(c)
functions
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Wed 16 May 2018 02:55:06 AM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from alexhenrie24 =at= gmail =dot= com



___

File Attachments:


---
Date: Wed 16 May 2018 02:55:06 AM UTC  Name: erf.patch  Size: 4KiB   By: psa



___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: [Bug-gsl] Bug in Hypergeometric function

[Patrick Alken]

Hello,

  I have logged your bug report into our bug tracker. It may take some 
time to develop a fix for this.


Thanks,
Patrick

On 05/01/2018 08:51 AM, Salehi Vaziri Kamran wrote:

Dear Sir or Madam,


Let me first thank you for your amazing package that is very well-written with 
a nice documentation.


I am using your package "gsl" for c++/c to calculate Integrals of 
Hypergeometric functions. I have trouble getting output for Hypergeometric function 
themselves.


To be precise your code gives error for the value of "gsl_sf_hyperg_2F1(-10, 2 , 0.5 
, 0.5)".


I tried to change the components a little bit and I understood that your code 
works down to a=-10 but after that it gives an error. I also check from Python 
and Mathematica and there is a finite value for this function(0.0778324).


Would you please respond to my e-mail since I was spending two days to debug 
this and now I get to this point. ?


Best regards,

Kamran Salehi Vaziri

[Bug-gsl] [bug #53904] Bug gsl_matrix_complex_set

[Patrick Alken]

URL:
  

 Summary: Bug gsl_matrix_complex_set
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Mon 14 May 2018 09:39:47 PM UTC
Category: None
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from daniel =dot= comparat =at= u-psud =dot= fr

Dear sir or Madam

I found a bug when using *gsl_matrix_complex_set.*

You will find the code below, the output of the data array is wrong if I use

/gsl_matrix_complex_set(m, 1, 1, gsl_complex_rect(1., 0.))/

in the code. So there is a bug. However if using

/gsl_complex x = gsl_complex_rect(1., 0.);/

/  gsl_matrix_complex_set(m, 1, 1, x);/


it is O.K.


Here is a code (I was testing using a real matrix the complex eigenvalue)




#include 
#include 
#include 
#include 
#include 
#include   /* printf */
#include  /* system, NULL, EXIT_FAILURE */
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 


using namespace std;



int main()
{


double S_VN =0. ;
int size = 4;
cout << " size = " << size << endl;
double data[16] = { -1.0, 1.0, -1.0, 1.0,
-8.0, 4.0, -2.0, 1.0,
27.0, 9.0, 3.0, 1.0,
64.0, 16.0, 4.0, 1.0
  };

for (int i=0; i<16; i++)
cout << " i = " << i << " data[i] " << data[i] << endl;


gsl_eigen_herm_workspace *w = gsl_eigen_herm_alloc(size); // This function
allocates a workspace for computing eigenvalues
gsl_matrix_complex *m = gsl_matrix_complex_alloc(size, size);
gsl_vector *eval = gsl_vector_alloc (size); // Valeurs propres

gsl_complex x = gsl_complex_rect(1., 0.);

for (int i=0; i<4; i++)
for (int j=0; j<4; j++)
{
cout << " i = " << i << " j = " << j << " j+4*i = " << j+4*i << "
data[i][j] " << data[j+4*i] << " data [9] " << data[9] << endl;
gsl_matrix_complex_set(m, 1, 1, x); // matrice pour rho_gg =
rho_11
}
gsl_eigen_herm(m, eval, w);// This function computes the eigenvalues
of the complex hermitian matrix
gsl_eigen_herm_free(w);   // This function frees the memory associated
with the workspace w.

return 1;
} 




___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

[Bug-gsl] [bug #53903] Test failure with gsl_sf_synchrotron_1_e on x86

[Patrick Alken]

URL:
  

 Summary: Test failure with gsl_sf_synchrotron_1_e on x86
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Mon 14 May 2018 09:38:11 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from awilfox =at= adelielinux =dot= org

gsl-2.4 fails its test suite on x86, when configured to build for
Pentium MMX:

FAIL: test
==

FAIL: gsl_sf_synchrotron_1_e(0.01, ) [54784666]
  expected: 4.4497250411421063e-01
  obtained: -1.8137993642342178e-02 ± 1.2082330897339797e-17
(rel=6.66134e-16)
  fracdiff: 1.0849884121238955e+00
 tolerance: 4.4408920985006262e-16
  value/expected not consistent within reported error
  value not within tolerance of expected value
  -1.813799364234217754e-02  1.208233089733979683e-17
FAIL: Synchrotron Functions [54784674]



Environment:

libc: musl libc (i386) Version 1.1.19

gcc: gcc (Adelie 6.4.0) 6.4.0

ld: GNU ld (GNU Binutils) 2.30

CPU: Intel(R) Pentium(R) 4 CPU 2.40GHz

uname -a: Linux zach 4.7.0-zach #1 SMP Sat Jul 30 20:37:18 CDT 2016 i686
GNU/Linux

CFLAGS: -O2 -march=pentium-mmx -mtune=pentium-m -fno-omit-frame-pointer





___

Reply to this item at:

  

___
  Message sent via Savannah
  https://savannah.gnu.org/

Re: [Bug-gsl] git source: doc/examples wants gaussfilt.c, but it does not exist

[Patrick Alken]

Hello,

  It should work now

Patrick

On 05/10/2018 10:46 AM, Hans-Bernhard Bröker wrote:
> Hello all,
>
> as of yesterday's commit 4cbb37d11189bc5bab3c8243bcfa45ff2ef75d22
> the git version of doc/examples/Makefile.am requests a non-existing
> source-file which was apparently supposed to be added as part of commit:
>
> $ LANG=C make -C doc/examples
> make: Entering directory
> '/home/hbbro/src/gnu/GSL/bld_cygwin/doc/examples'
> make: *** No rule to make target 'gaussfilt.c', needed by 'all-am'. 
> Stop.
> make: Leaving directory '/home/hbbro/src/gnu/GSL/bld_cygwin/doc/examples'
>
>

Re: [Bug-gsl] typo error

[Patrick Alken]

Thanks, I have fixed it on the git

On 04/21/2018 05:54 PM, Alex Henrie wrote:
> On Sat, Apr 21, 2018 at 4:42 PM Patrick Alken <al...@colorado.edu> wrote:
>> Err, sorry I was looking at the wrong PDF. On page 534, I see a
>> discussion of Autoconf macros, still no equation for 'e'
> He was counting from the beginning of the PDF, he meant page 520.
>
> -Alex

Re: [Bug-gsl] typo error

[Patrick Alken]

Err, sorry I was looking at the wrong PDF. On page 534, I see a
discussion of Autoconf macros, still no equation for 'e'

On 04/21/2018 04:39 PM, Patrick Alken wrote:
> Hello, can you give the chapter/section number? On page 534, I see an
> example program for sparse linear algebra (section 43.3) and I don't see
> any equation for 'e' on this page.
>
> On 04/21/2018 02:05 PM, Србислав Нешић wrote:
>> In PDF document
>>
>> https://www.gnu.org/software/gsl/doc/latex/gsl-ref.pdf
>>
>> page 534
>>
>> e = *1* + 1/2! + 1/3! + 1/4! + ...
>>
>> should be
>>
>> e = *2* + 1/2! + 1/3! + 1/4! + ...
>>
>> or
>>
>> e = *1 + 1* + 1/2! + 1/3! + 1/4! + ...
>>
>> Regards, Srba
>
>

Re: [Bug-gsl] typo error

[Patrick Alken]

Hello, can you give the chapter/section number? On page 534, I see an
example program for sparse linear algebra (section 43.3) and I don't see
any equation for 'e' on this page.

On 04/21/2018 02:05 PM, Србислав Нешић wrote:
> In PDF document
>
> https://www.gnu.org/software/gsl/doc/latex/gsl-ref.pdf
>
> page 534
>
> e = *1* + 1/2! + 1/3! + 1/4! + ...
>
> should be
>
> e = *2* + 1/2! + 1/3! + 1/4! + ...
>
> or
>
> e = *1 + 1* + 1/2! + 1/3! + 1/4! + ...
>
> Regards, Srba

Re: [Bug-gsl] incorrect convergence check result

[Patrick Alken]

Hi Edwin,

  Interestingly, that Brown test problem is part of the test suite for 
nlinear (see multifit_nlinear/test_brown3.c). It looks like I set gtol = 
GSL_DBL_EPSILON^(0.9) for the test suite problems, so it is much 
stricter than the recommended value in the manual.


  As a general rule, determining convergence for optimization problems 
is more of an art than a science, which is why GSL offers 2 seperate 
rules (based on the step size and the gradient). I had tried to add a 
third rule based on the change in ||f|| from one iteration to the next, 
but I found this didn't work so well so I disabled it.


  If the gradient test isn't working well, you could disable it by 
setting gtol = 0, and then only the xtol test will be used.


  But in any case I will take a look at this Brown test problem to try 
to understand what is going on.


Thanks,
Patrick

On 03/22/2018 08:02 AM, Sarkissian, Edwin (398E) wrote:

Hello Patrick,

After I studied the gradient-test algorithm implemented in 
gsl_multifit_nlinear_test, I realized that the algorithm avoids some problems 
by scaling the gradient; however, the scaling method results in another 
problem.   Adding a constant to a cost function does not change the shape (ups 
and downs, minima and maxima, ...) of the cost function, and minima will stay 
where they are.   However, the gradient-base convergence test given at 
https://www.gnu.org/software/gsl/doc/html/nls.html#testing-for-convergence 
could produce false results if one simple adds a large constant to the cost 
function.  For example, for the same x and the same gtol, the test may indicate 
no-convergence for
( ||f(x)||^2 )/2
but indicate convergence for
( ||f(x)||^2 + ||C||^2 )/2
where C is constant.

I encountered such a problem when I was trying to find the minimum of a test 
nonlinear least square problem due to Brown.  It is the number 4 problem in 
Section 3 of the attached document.  For the gtol suggested at 
https://www.gnu.org/software/gsl/doc/html/nls.html#testing-for-convergence, the 
gradient-based convergence test gives a false positive (converged) result.

Regards,
Edwin



On 3/21/18, 12:26 PM, "Patrick Alken" <al...@colorado.edu> wrote:

 Hello,
 
Can you provide a minimal working example program which illustrates

 the issue?
 
 Thanks,

 Patrick
 
 On 03/21/2018 08:02 AM, Sarkissian, Edwin (398E) wrote:

 > Dear Madam/Sir,
 >
 > Sometimes, the function
 > gsl_multifit_nlinear_test
 > sets the argument info to 2 (converged based on gradient check) when the 
gradient is very large.  For example, in one case the gradient is [-99, 
0.179997].
 >
 > I have noticed the same issue with the function
 >  gsl_multifit_fdfsolver_test
 >
 > Thank you,
 > Edwin Sarkissian
 
 
 

Re: [Bug-gsl] incorrect convergence check result

[Patrick Alken]

Hello,

  Can you provide a minimal working example program which illustrates 
the issue?


Thanks,
Patrick

On 03/21/2018 08:02 AM, Sarkissian, Edwin (398E) wrote:

Dear Madam/Sir,

Sometimes, the function
gsl_multifit_nlinear_test
sets the argument info to 2 (converged based on gradient check) when the 
gradient is very large.  For example, in one case the gradient is [-99, 
0.179997].

I have noticed the same issue with the function
 gsl_multifit_fdfsolver_test

Thank you,
Edwin Sarkissian

Re: [Bug-gsl] please check lib gsl_multifit_nlinear.h for bugs

[Patrick Alken]

Hello,

  The website documentation and examples are for GSL v2.4, so they will 
likely not work with v2.1 since the nonlinear least squares code has 
changed significantly in recent versions. Can you please install v2.4 
and try the example program again? Also please send along any example 
code you cannot successfully compile.


Patrick

On 03/20/2018 05:30 AM, Dr. Jonny Birkhan wrote:

Dear colleagues,

after hours of trying to get your first multifit-nonlinear expamle 
running (https://www.gnu.org/software/gsl/doc/html/nls.html#examples), 
-including request for help along the help-gsl@...-,  I conclude that 
there is something wrong with the library.


I have also used a fresh ubuntu installation. Then I installed the gsl 
libs 2.1 from the synaptic package manager.


Copying the example code from the website and compiling it, always fails.

Compiling and running the linear fit example works well.

So, there seems to be something wrong.

Could you check the library in that naive way any non-expert user 
would do?


Currently I cannot recompile code written a year ago.

Thanks a lot in advance and best regards,

Jonny Birkhan

Re: [Bug-gsl] error in multifit/multiwlinear.c in gsl-2.3 and gsl-2.4

[Patrick Alken]

Hello,

  I responded to this yesterday. The issue was the
gsl_multifit_wlinear_svd function is deprecated (replaced by
gsl_multifit_wlinear_tsvd), and I hadn't tested the old function. The
issue should be fixed now in the git repository. Could you checkout the
latest git code and test your program?

Thanks,
Patrick

On 02/11/2018 12:44 PM, Vlad Koli wrote:
> Hi,
>
> There is an error in file multifit/multiwlinear.c. It is present at least  
> in sources for gsl-2.3, gsl-2.4. 
>
>   Line 69 should be replaced (I think) from:
>
>  else if (X->size1 != work->n || X->size2 != work->p)
>
>to
>
>  else if (X->size1 != work->nmax || X->size2 != work->pmax)  // vk:
>
> Values of "work->n", "work->p" for this file, for example, are nulls.
>
> After this replacement, re-compilation, and re-installation the tests 
> began to work correctly.
>
>
> Check, for example the test file "linfit-pseudoinv.c" quoted below.  
>
> ** compilation:
>
>  cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
>
> **  execution: 
>
>  ./a
>
> ** result:
>
>> gsl: multiwlinear.c:73: ERROR: size of workspace does not match size of 
>> observation matrix
>> Default GSL error handler invoked.
>> Aborted
>
> Regards,
> Vladimir Koliadin
>
>
> - begin file: linfit-pseudoinv.c:-
> // compile:
> // cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
> #include 
> #include 
>
>
> double tol = 1.e-08; // tolerance
>
> main()
> {
>  
> int n = 3; // observations
> int m = 2; // parameters 
>
> size_t rank;   // rank
> double chisq;  // sum of squares of residuals
>
> // allocate:
> gsl_vector *x = gsl_vector_alloc(m); // solution 
> gsl_vector *y = gsl_vector_alloc(n); // data 
> gsl_vector *w = gsl_vector_alloc(n); // weights
>
> gsl_matrix *A = gsl_matrix_alloc(n, m);
> gsl_matrix *cov = gsl_matrix_alloc(m, m); // covariances
>
> gsl_multifit_linear_workspace* wsp = gsl_multifit_linear_alloc (n, m);
>
>
> // set matrix
> gsl_matrix_set(A, 0,0, 1.);
> gsl_matrix_set(A, 0,1, 1.);
> gsl_matrix_set(A, 1,0, 2.);
> gsl_matrix_set(A, 1,1, 2.);
> gsl_matrix_set(A, 2,0, 3.);
> gsl_matrix_set(A, 2,1, 3.);
>
> // set data 
> gsl_vector_set(y, 0, 1.);
> gsl_vector_set(y, 1, 2.);
> gsl_vector_set(y, 2, 3.1);
>
> // set weights
> gsl_vector_set_all(w,1.);
> gsl_vector_set(w, 0, 0.);
> gsl_vector_set(w, 2, 100.);
>
>
> gsl_multifit_wlinear_svd (A, w, y,  tol, , x, cov, , wsp);
>  
>
> printf("solution: x=%.4f  y=%.4f   (rank=%d, chi2=%.3f)\n", 
> gsl_vector_get(x,0), gsl_vector_get(x,1), rank, chisq);
>
>
> }
> - end file: linfit-pseudoinv.c:-
>
>
>

Re: [Bug-gsl] multiwlinear.c:73: ERROR: size of workspace does not match size of observation matrix

[Patrick Alken]

Vlad,

  I have uploaded a fix for this issue to the git repository. Now when I
run your original program I get this output:

$ ./test
solution: x=0.5166  y=0.5166   (rank=1, chi2=0.004)

Thanks,
Patrick

On 02/10/2018 10:52 AM, Patrick Alken wrote:
> So as it turns out, the gsl_multifit_wlinear_svd function should be
> deprecated, and I removed it from the documentation, but kept it around
> for legacy reasons - without testing it properly!
>
> The new function is called gsl_multifit_wlinear_tsvd, and you can use
> this call in your code:
>
> gsl_multifit_wlinear_tsvd(A, w, y, tol, x, cov, , , wsp);
>
> This should make your program run correctly. In the meantime I will
> figure out what to do with gsl_multifit_wlinear_svd - maybe just make it
> a wrapper function for the new tsvd call.
>
> Sorry for the trouble,
> Patrick
>
> On 02/10/2018 10:28 AM, Patrick Alken wrote:
>> Thanks for reporting this, it is indeed a bug. My apologies, I had
>> rewritten the linear least squares routines recently and didn't make
>> proper tests for the wlinear_svd routine. I will work on a fix asap. In
>> the meantime, you can replace the gsl_multifit_wlinear_svd call with:
>>
>> gsl_multifit_wlinear(A, w, y, x, cov, , wsp);
>>
>> which should work as a temporary fix.
>>
>> On 02/10/2018 12:26 AM, Vlad Koli wrote:
>>> Hi,
>>>
>>> In Debian 9  (gsl 2.3) compilation/execution of the following test-file 
>>> causes 
>>> error messages. In Debian 7 (gsl ???)  it worked fine. The same error 
>>> message
>>> also appears in other my programs that worked well in Debian 7 and 
>>> were recompiled in Debian 9.
>>>
>>> ** compilation:
>>>
>>>  cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
>>>
>>> **  execution: 
>>>
>>>  ./a
>>>
>>> ** result:
>>>
>>>> gsl: multiwlinear.c:73: ERROR: size of workspace does not match size of 
>>>> observation matrix
>>>> Default GSL error handler invoked.
>>>> Aborted
>>> Thank you in advance.
>>>
>>> Vladimir
>>>
>>> - begin file: linfit-pseudoinv.c:-
>>> // compile:
>>> // cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
>>> #include 
>>> #include 
>>>
>>>
>>> double tol = 1.e-08; // tolerance
>>>
>>> main()
>>> {
>>>  
>>> int n = 3; // observations
>>> int m = 2; // parameters 
>>>
>>> size_t rank;   // rank
>>> double chisq;  // sum of squares of residuals
>>>
>>> // allocate:
>>> gsl_vector *x = gsl_vector_alloc(m); // solution 
>>> gsl_vector *y = gsl_vector_alloc(n); // data 
>>> gsl_vector *w = gsl_vector_alloc(n); // weights
>>>
>>> gsl_matrix *A = gsl_matrix_alloc(n, m);
>>> gsl_matrix *cov = gsl_matrix_alloc(m, m); // covariances
>>>
>>> gsl_multifit_linear_workspace* wsp = gsl_multifit_linear_alloc (n, m);
>>>
>>>
>>> // set matrix
>>> gsl_matrix_set(A, 0,0, 1.);
>>> gsl_matrix_set(A, 0,1, 1.);
>>> gsl_matrix_set(A, 1,0, 2.);
>>> gsl_matrix_set(A, 1,1, 2.);
>>> gsl_matrix_set(A, 2,0, 3.);
>>> gsl_matrix_set(A, 2,1, 3.);
>>>
>>> // set data 
>>> gsl_vector_set(y, 0, 1.);
>>> gsl_vector_set(y, 1, 2.);
>>> gsl_vector_set(y, 2, 3.1);
>>>
>>> // set weights
>>> gsl_vector_set_all(w,1.);
>>> gsl_vector_set(w, 0, 0.);
>>> gsl_vector_set(w, 2, 100.);
>>>
>>>
>>> gsl_multifit_wlinear_svd (A, w, y,  tol, , x, cov, , wsp);
>>>  
>>>
>>> printf("solution: x=%.4f  y=%.4f   (rank=%d, chi2=%.3f)\n", 
>>> gsl_vector_get(x,0), gsl_vector_get(x,1), rank, chisq);
>>>
>>>
>>> }
>>> - end file: linfit-pseudoinv.c:-
>>>
>>>
>>>
>

Re: [Bug-gsl] multiwlinear.c:73: ERROR: size of workspace does not match size of observation matrix

[Patrick Alken]

So as it turns out, the gsl_multifit_wlinear_svd function should be
deprecated, and I removed it from the documentation, but kept it around
for legacy reasons - without testing it properly!

The new function is called gsl_multifit_wlinear_tsvd, and you can use
this call in your code:

gsl_multifit_wlinear_tsvd(A, w, y, tol, x, cov, , , wsp);

This should make your program run correctly. In the meantime I will
figure out what to do with gsl_multifit_wlinear_svd - maybe just make it
a wrapper function for the new tsvd call.

Sorry for the trouble,
Patrick

On 02/10/2018 10:28 AM, Patrick Alken wrote:
> Thanks for reporting this, it is indeed a bug. My apologies, I had
> rewritten the linear least squares routines recently and didn't make
> proper tests for the wlinear_svd routine. I will work on a fix asap. In
> the meantime, you can replace the gsl_multifit_wlinear_svd call with:
>
> gsl_multifit_wlinear(A, w, y, x, cov, , wsp);
>
> which should work as a temporary fix.
>
> On 02/10/2018 12:26 AM, Vlad Koli wrote:
>> Hi,
>>
>> In Debian 9  (gsl 2.3) compilation/execution of the following test-file 
>> causes 
>> error messages. In Debian 7 (gsl ???)  it worked fine. The same error message
>> also appears in other my programs that worked well in Debian 7 and 
>> were recompiled in Debian 9.
>>
>> ** compilation:
>>
>>  cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
>>
>> **  execution: 
>>
>>  ./a
>>
>> ** result:
>>
>>> gsl: multiwlinear.c:73: ERROR: size of workspace does not match size of 
>>> observation matrix
>>> Default GSL error handler invoked.
>>> Aborted
>> Thank you in advance.
>>
>> Vladimir
>>
>> - begin file: linfit-pseudoinv.c:-
>> // compile:
>> // cc -o a -lm -lgsl -lgslcblas linfit-pseudoinv.c 
>> #include 
>> #include 
>>
>>
>> double tol = 1.e-08; // tolerance
>>
>> main()
>> {
>>  
>> int n = 3; // observations
>> int m = 2; // parameters 
>>
>> size_t rank;   // rank
>> double chisq;  // sum of squares of residuals
>>
>> // allocate:
>> gsl_vector *x = gsl_vector_alloc(m); // solution 
>> gsl_vector *y = gsl_vector_alloc(n); // data 
>> gsl_vector *w = gsl_vector_alloc(n); // weights
>>
>> gsl_matrix *A = gsl_matrix_alloc(n, m);
>> gsl_matrix *cov = gsl_matrix_alloc(m, m); // covariances
>>
>> gsl_multifit_linear_workspace* wsp = gsl_multifit_linear_alloc (n, m);
>>
>>
>> // set matrix
>> gsl_matrix_set(A, 0,0, 1.);
>> gsl_matrix_set(A, 0,1, 1.);
>> gsl_matrix_set(A, 1,0, 2.);
>> gsl_matrix_set(A, 1,1, 2.);
>> gsl_matrix_set(A, 2,0, 3.);
>> gsl_matrix_set(A, 2,1, 3.);
>>
>> // set data 
>> gsl_vector_set(y, 0, 1.);
>> gsl_vector_set(y, 1, 2.);
>> gsl_vector_set(y, 2, 3.1);
>>
>> // set weights
>> gsl_vector_set_all(w,1.);
>> gsl_vector_set(w, 0, 0.);
>> gsl_vector_set(w, 2, 100.);
>>
>>
>> gsl_multifit_wlinear_svd (A, w, y,  tol, , x, cov, , wsp);
>>  
>>
>> printf("solution: x=%.4f  y=%.4f   (rank=%d, chi2=%.3f)\n", 
>> gsl_vector_get(x,0), gsl_vector_get(x,1), rank, chisq);
>>
>>
>> }
>> - end file: linfit-pseudoinv.c:-
>>
>>
>>
>

[Bug-gsl] [bug #52927] make check fails on Bessel j2 test

[Patrick Alken]

URL:
  <http://savannah.gnu.org/bugs/?52927>

 Summary: make check fails on Bessel j2 test
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Thu 18 Jan 2018 11:13:07 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

from eshell =at= ucsc =dot= edu

subject: GSL 2.4 make check failing on Rocks cluster running CentOS 6.9

note: it looks like the error term (r) is causing the failure, not the
computed j2 value itself



Thanks much for the quick reply!  I incremented the tolerance value, but it
is still failing even at TEST_TOL6:

==
   gsl 2.4: specfunc/test-suite.log
==

# TOTAL: 1
# PASS:  0
# SKIP:  0
# XFAIL: 0
# FAIL:  1
# XPASS: 0
# ERROR: 0

.. contents:: :depth: 2

FAIL: test
==

FAIL: gsl_sf_bessel_j2_e(1048576.0, ) [168]
  expected: -3.1518539455252412e-07
  obtained: -3.1518539455252539e-07 ± 2.7994086564622246e-22
(rel=8.88178e-16)
  fracdiff: 2.0155588470164931e-15
 tolerance: 2.3283064365386963e-10
  value/expected not consistent within reported error
  -3.151853945525253879e-07  2.799408656462224591e-22
FAIL: Bessel Functions [407]

On Thu, Jan 18, 2018 at 3:01 PM, Patrick Alken <al...@colorado.edu> wrote:

> Hello,
>
>   It looks like its calculating the value correctly, but the test
> tolerance needs to be relaxed slightly. Can you locate this line in
> specfun/test_bessel.c:
>
>  186   TEST_SF(s,  gsl_sf_bessel_j2_e, (1048576.0, ),
> -3.1518539455252413111e-07, TEST_TOL3, GSL_SUCCESS);
>
> and change the TEST_TOL3 to TEST_TOL4, and let me know if the test passes?
> If not try TEST_TOL5 and then TEST_TOL6, and tell me which one allows the
> test to pass.
>
> Thanks,
> Patrick
>
>
> On 01/18/2018 03:56 PM, Eric Shell wrote:
>
>> Hello,
>>
>> I'm trying to install GSL 2.4 on a Rocks cluster running CentOS 6.9, gcc
>> version 4.4.7.  I am running configure with just a --prefix argument.
>> make
>> succeeds without errors, but make check is failing on the specfunc test.
>> Here are the contents of the log file:
>>
>> 
>> 
>>
>> ==
>> gsl 2.4: specfunc/test-suite.log
>> ==
>>
>> # TOTAL: 1
>> # PASS:  0
>> # SKIP:  0
>> # XFAIL: 0
>> # FAIL:  1
>> # XPASS: 0
>> # ERROR: 0
>>
>> .. contents:: :depth: 2
>>
>> FAIL: test
>> ==
>>
>> FAIL: gsl_sf_bessel_j2_e(1048576.0, ) [168]
>>expected: -3.1518539455252412e-07
>>obtained: -3.1518539455252539e-07 ± 2.7994086564622246e-22
>> (rel=8.88178e-16)
>>fracdiff: 2.0155588470164931e-15
>>   tolerance: 4.5474735088646412e-13
>>value/expected not consistent within reported error
>>-3.151853945525253879e-07  2.799408656462224591e-22
>> FAIL: Bessel Functions [407]
>>
>> 
>> 
>>
>> Is this a known issue?  How can I address the underlying issue?
>>
>> Thanks!
>>
>> - Eric
>>





___

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52927>

___
  Message sent via/by Savannah
  http://savannah.gnu.org/

Re: [Bug-gsl] FFT not generated

[Patrick Alken]

Hi,

  Did the library pass 'make check'. What happens if you do:

cd fft ; make test ; ./test

On 01/05/2018 05:25 PM, Peter Schacht wrote:
> Hello,
>
>  
>
> I installed the GSL in win10 64bit -system with minGw  and msys using the
> cofigure script. The installaion seemed to work fine and  many functions
> where
>
> generated and also the standard tests worked fine.  I wanted to test my old
> FFT programs and became the message "undefined reference to 
>
> `gsl_fft_real_radix2_transform'" checking both the libs libgsl.a and
> libgslcblas.a  showed that this function was not in there, only  fft.o and
> dft.o.
>
> A look in the build directory showed that only fft.o and dft.o have been
> produced in built dir fft. The header files where shifted to dir gsl in
> built but nothing built was built there.
>
>  
>
> Makefile.am  in dir  fft  looks not correct.
>
> attached you'lll find the config.log 
>
>  
>
> Best Regards
>
> Peter Schacht
>
>
>
> 258/127 Moo 3, Thepkrasatri Road
>
> Srisoonthorn, Thalang, Phuket 83110
>
>  
>
> Mobil +66(96)926 83 62
>
> Email    peter.scha...@schacht-co.de 
>
>  
>
>  
>
>  
>

Re: [Bug-gsl] Minor improvement to gsl_spmatrix.h suggested

[Patrick Alken]

Manuel,

  Sounds like a good suggestion to me - I've made the change on the git 
and hopefully didn't break anything. If you're currently using 
gsl_spmatrix in your work and can test the update I'd appreciate it.


Thanks,
Patrick

On 12/22/2017 09:44 AM, Schmitz Manuel (LBC) wrote:

Dear Sirs,

in gsl_spmatrix.h, the struct "gsl_spmatrix" gets defined:

typedef struct
{
 ...
 void* work;
 ...
} gsl_spmatrix;

I am suggesting to replace the member "work" by a union:

typedef struct
{
 ...
 union {
void * work;
size_t * work_sze;
double * work_dbl;
};
...
} gsl_spmatrix;

This has the following benefits:

-  This can avoid (nasty) casts when using the working memory later. 
Unions are considered to be safer than casts.

-  This can avoid compiler warnings about casts when using high warning 
levels.

-  It better expresses the intent: "work" shall be used as a buffer for "double" 
or for "size_t".

-  The change does not break the existing interface.

-  The runtime overhead is zero and the compile-time overhead is 
neglectable.

Best regards

i.A. Manuel Schmitz

TE-Krantechnik/Berechnung

Liebherr-Werk Biberach GmbH
PO Box 16 63
88396 Biberach an der Riss
Germany
Phone: +49 (7351) 414639
Fax: +49 (7351) 412879
manuel.schm...@liebherr.com
www.liebherr.com

Chairman of the supervisory board: Stefan Heissler
Managing directors: Dipl.-Ing. (FH) Marco Guariglia, Dipl.-Ing. Günther 
Hardock, Dr. Thomas Schwaninger, Dipl.-Ing. Dominique Tasch
Registered business address: 88400 Biberach an der Riss, Memminger Straße 120
Court of jurisdiction: Amtsgericht Ulm HRB 640075, USt-Id Nr. DE811120028

This email may contain confidential and/or privileged information. If you are 
not the intended recipient (or have received this email in error) please notify 
the sender immediately and destroy this email. Any unauthorized copying, 
disclosure or distribution of the material in this email is strictly forbidden.

[Bug-gsl] [bug #52570] Inaccuracy of the Airy function due to invocation of GSL's cosine function with large input parameters

[Patrick Alken]

URL:
  

 Summary: Inaccuracy of the Airy function due to invocation of
GSL's cosine function with large input parameters
 Project: GNU Scientific Library
Submitted by: psa
Submitted on: Fri 01 Dec 2017 05:32:49 PM UTC
Category: Runtime error
Severity: 3 - Normal
Operating System: 
  Status: None
 Assigned to: None
 Open/Closed: Open
 Release: 
 Discussion Lock: Any

___

Details:

zhoulai =dot= fu =at= gmail =dot= com

Dear GSL developers,



This is a follow-up email of my previous report of regarding GSL’s Airy and
cosine functions
(https://lists.gnu.org/archive/html/bug-gsl/2017-11/msg00011.html).


Here I would like to report another kind of unexpected results when using
the Airy function gsl_sf_Ai_e.  Code to reproduce the issue:


#include 

#include 

#include 

#include 



int main(){

  gsl_sf_result res;

  double   x = -1.14e34;

  gsl_mode_t mode = GSL_PREC_DOUBLE;

  int status = gsl_sf_airy_Ai_e(x, mode, );

  printf("status = %d, val = %g, err = %g\n", status,res.val,res.err);

}



Running the code gives:  status = 0, val = -inf, err = inf



The calculation result here, val = -inf, is incorrect: Airy functions are
damped oscillatory for negative x; using Mathematica you get AiryAi
[-1.1e34]=-1.36e-9.  Besides, the error estimate, inf, is obviously
imprecise, and the returned status, GSL_SUCCESS 0, may be somewhat
misleading.



Cause of the issue in my humble opinion: Airy_ai (-1.1e34) invokes
gsl_sf_cos_err_e (theta.val, theta.err, _result) where theta.val =
-8.1145794715437919e+50, theta.err = 7.4985953321244595e+35 (according to
gdb). After this invocation, cos_result.val and cos_result.err become -inf
and inf respectively, which is clearly wrong, which stays away from its
[-1,1] bound.  Since “it is known that the GSL gsl_sf_cos and gsl_sf_sin
functions fail for large inputs” (from https://lists.gnu.org/archive/
html/bug-gsl/2017-11/msg00012.html), maybe GSL's Airy functions, both the
first and the second kinds, should consider use libm’s cosine instead.


Attached is the same code as above.


Thanks.


Zhoulai Fu




___

Reply to this item at:

  

___
  Message sent via/by Savannah
  http://savannah.gnu.org/