Re: Deprecated Basis Splines functions

[Patrick Alken]

Hello, those functions should still be available in version 2.8, but as 
you note they are marked for deprecation in future releases.


Are you not able to use those functions in 2.8?

On 7/5/24 10:39, Andrew Forembski via Users list for GNU Scientific 
Library (GSL) help wrote:

[External email - use caution]


Hi,

Up until the release of gsl 2.8 I was using the functions:
gsl_bspline_eval_nonzero(...)
gsl_bspline_deriv_eval_nonzero(...)

to evaluate B-splines functions and their derivatives at specific points
along an axis in my program.
These functions seem to have been removed as of gsl 2.8, I have been able
to replace them with:
gsl_bspline_basis(...)
gsl_bspline_basis_deriv(...)

My issue is that there seems to be no note on deprecated functions in the
2.8 documentation and I don't seem to be able to find the documentation for
gsl < 2.8.
Are the functions gsl_bspline_basis and gsl_bspline_basis_deriv available
in earlier versions?
I'm asking this because I'm trying to keep my code portable and would
ideally prefer not requiring gsl to be updated to 2.8 as some users may be
on systems where this has to be done manually.

Would you have any suggestions regarding this problem?

Kind regards,
Andrew

*Andrew Forembski*

Government of Ireland Postgraduate Scholar GOIPG/2018/1070

Web: http://www.research.ie/

@IrishResearch

--
__

Séanadh Ríomhphoist/_

Email Disclaimer__
**

Tá an ríomhphost seo agus
aon chomhad a sheoltar leis faoi rún agus is lena úsáid ag an seolaí agus
sin amháin é. Is féidir tuilleadh a léamh anseo.

*
_

This e-mail and any
files transmitted with it are confidential and are intended solely for use
by the addressee. Read more here.
 _
*_

Re: how to debug multifit_nlinear

[Patrick Alken]

Hello, it looks like a problem with the finite-difference approximation 
to the Jacobian matrix. Can you provide a MWE we can test?


On 6/19/24 03:23, PICCA Frederic-Emmanuel wrote:

[External email - use caution]


Hello,

I am using

#include 

in order to do a least square fitting of this function

int model_f (const gsl_vector *x, void *data, gsl_vector *f)
{
 size_t n = ((struct data *)data)->n;
 float *t = ((struct data *)data)->t;
 float *y = ((struct data *)data)->y;
 int* lista = ((struct data *)data)->lista;
 size_t n_lista = ((struct data *)data)->n_lista;
 float* parameters = ((struct data *)data)->parameters;
 size_t n_parameters = ((struct data *)data)->n_parameters;
 FitFunction funcs = ((struct data *)data)->funcs;

 float dyda;

 size_t i;

 for (i=0; i

GSL use within NASA

[Patrick Alken]

Hello all,

  If anyone is a NASA employee using GSL, or if anyone knows of a 
current NASA mission using GSL, can you please contact me?


Thanks,

Patrick

GSL 2.8 Released

[Patrick Alken]

Dear GSL community,

  Version 2.8 of the GNU Scientific Library (GSL) has been released. 
Thank you to all who helped test the library prior to the release, and 
thank you to everyone for using the library and giving feedback and 
reports. The following changes have been added to the library:


* What is new in gsl-2.8:

** apply patch for bug #63679 (F. Weimer)

** updated multilarge TSQR method to store ||z_2|| and
   provide it to the user

** add routines for Hermite B-spline interpolation

** fix for bug #59624

** fix for bug #59781 (M. Dunlap)

** bug fix #61094 (reported by A. Cheylus)

** add functions:
   - gsl_matrix_complex_conjugate
   - gsl_vector_complex_conj_memcpy
   - gsl_vector_complex_div_real
   - gsl_linalg_QR_lssolvem_r
   - gsl_linalg_complex_QR_lssolvem_r
   - gsl_linalg_complex_QR_QHmat_r
   - gsl_linalg_QR_UR_lssolve
   - gsl_linalg_QR_UR_lssvx
   - gsl_linalg_QR_UR_QTvec
   - gsl_linalg_QR_UU_lssvx
   - gsl_linalg_QR_UD_lssvx
   - gsl_linalg_QR_UD_QTvec
   - gsl_linalg_complex_cholesky_{decomp2,svx2,solve2,scale,scale_apply}
   - gsl_linalg_SV_{solve2,lssolve}
   - gsl_rstat_norm

** add Lebedev quadrature (gsl_integration_lebedev)

** major overhaul to the B-spline module to add
   new functionality

Re: Test release for GSL 2.8

[Patrick Alken]

Thank you all for the testing reports you submitted. I think I have 
addressed the latest reports, though I may have missed some things. I 
uploaded a new test release here:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.102.tar.gz

I would appreciate if those who found warnings before could test again 
on this file.


Patrick

On 5/14/24 16:42, Patrick Alken wrote:
Thank you all again for your reports. I have (I think) fixed all the 
compiler warnings related to filter/movstat, as well as the fabs 
issues. I would greatly appreciate it if you could do another round of 
testing on this file:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.101.tar.gz

Thanks,

Patrick

On 5/12/24 13:34, Patrick Alken wrote:

[External email - use caution]


All, thank you for your testing reports. I believe I have addressed all
the compiler warnings and other issues raised. I have uploaded a new
test release here:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.100.tar.gz 



I would greatly appreciate if everyone can run the tests again:

./configure && make && make check

Thanks,

Patrick


On 5/10/24 16:04, Patrick Alken wrote:

[External email - use caution]


Dear all,

  It is time to make a new GSL release. I have uploaded a test release
to:

ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz 



ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz.sig 




All reports are welcome - anyone who can test on various platforms 
would

be appreciated (Linux, BSD, Mac OS, Windows).

Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need
to install the python-based sphinx software first:

1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than
system-wide)

3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: Test release for GSL 2.8

[Patrick Alken]

Thank you all again for your reports. I have (I think) fixed all the 
compiler warnings related to filter/movstat, as well as the fabs issues. 
I would greatly appreciate it if you could do another round of testing 
on this file:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.101.tar.gz

Thanks,

Patrick

On 5/12/24 13:34, Patrick Alken wrote:

[External email - use caution]


All, thank you for your testing reports. I believe I have addressed all
the compiler warnings and other issues raised. I have uploaded a new
test release here:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.100.tar.gz 



I would greatly appreciate if everyone can run the tests again:

./configure && make && make check

Thanks,

Patrick


On 5/10/24 16:04, Patrick Alken wrote:

[External email - use caution]


Dear all,

  It is time to make a new GSL release. I have uploaded a test release
to:

ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz 



ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz.sig 




All reports are welcome - anyone who can test on various platforms would
be appreciated (Linux, BSD, Mac OS, Windows).

Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need
to install the python-based sphinx software first:

1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than
system-wide)

3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: Test release for GSL 2.8

[Patrick Alken]

All, thank you for your testing reports. I believe I have addressed all 
the compiler warnings and other issues raised. I have uploaded a new 
test release here:



ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.100.tar.gz

I would greatly appreciate if everyone can run the tests again:

./configure && make && make check

Thanks,

Patrick


On 5/10/24 16:04, Patrick Alken wrote:

[External email - use caution]


Dear all,

  It is time to make a new GSL release. I have uploaded a test release 
to:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz 

ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz.sig 



All reports are welcome - anyone who can test on various platforms would
be appreciated (Linux, BSD, Mac OS, Windows).

Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need
to install the python-based sphinx software first:

1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than
system-wide)

3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: Test release for GSL 2.8

[Patrick Alken]

Mark,

  Thanks for looking at this. A while ago I began to rewrite the 
cspline interpolation module, to allow it to work for the N=2 case, and 
also to make it more cache-efficient when evaluating the spline. 
Unfortunately I never completely finished it and so 2 days ago I 
reverted everything back to the original v2.7 code (commit 
ac250bf26d80f723319acde73c4f42475d81fc98). So the test_cspline4() is 
gone now, and everything should be like it was at the v2.7 release.


Can you make sure you are working with the latest git code? After the 
2.8 release I will return to this and try to complete what I want to do!


Best,

Patrick

On 5/11/24 14:42, Mark Galassi wrote:

[External email - use caution]


Dear Patrick and help-gsl,

I had found, a few days ago, a test failure in interpolation/test.c.

Then just yesterday (May 10) there was a revert of interpolation/test.c to a 
previous version, which yanked out the test failure I had been investigating.

But the problem still exists, so let me outline here what I found and make some 
suggestions.

The test_cspline4() in git hash e647c6484702b638062d6e2154055f6535790cd7 has a 
failure.

The problem has two parts: what brings it out, and what

1. test_cspline4() calls make_xy_table() with an argument of 2, which is not OK 
- I guess you need at least 3 to do csplines, and things go wrong down the line 
when you start doing linear algebra.

The test should be updated to have an argument of 3, or the cspline subsystem 
should be updated to warn about that.

2. the real problem is that all the intervening function calls do not check for 
that size properly.  At one point two gets subtracted from it, and we end up 
with calls to malloc(0).  Note that if by chance we had called with 1 or 0, we 
might end up with malloc(-1) or malloc(-2).

The behavior of malloc(0) surprised me: I had thought it would just return 
NULL, but the man page states:

"""The memory is not initialized.  If size is 0, then malloc() returns either NULL, or a 
unique pointer value that can  later be successfully passed to free()."""

which means that this code in tridiag.c:

static
int
solve_tridiag(
   const double diag[], size_t d_stride,
   const double offdiag[], size_t o_stride,
   const double b[], size_t b_stride,
   double x[], size_t x_stride,
   size_t N)
{
   int status = GSL_SUCCESS;
   double *gamma = (double *) malloc (N * sizeof (double));
   double *alpha = (double *) malloc (N * sizeof (double));
   double *c = (double *) malloc (N * sizeof (double));
   double *z = (double *) malloc (N * sizeof (double));

   if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
 {
   GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
 }
   else

does not capture this kind of error, and we end up with bad access, since the 
interpolation index arithmetic is not up to handling a size of 2 (and hence N 
of 0).

I would suggest some things:

a. gsl vectors allow a size of zero, so we should not trap at the level of the 
gsl_vector subsystem

b. tridiag.c needs to require N > 0.  Right now it looks at the return values 
of malloc():

   double *gamma = (double *) malloc (N * sizeof (double));
   double *alpha = (double *) malloc (N * sizeof (double));
   double *c = (double *) malloc (N * sizeof (double));
   double *z = (double *) malloc (N * sizeof (double));

   if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
 {
   GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
 }
   else

and then it does unsigned arithmetic here:

   for (i = 1; i < N - 1; i++)
 {
   alpha[i] = diag[d_stride * i] - offdiag[o_stride*(i - 1)] * gamma[i 
- 1];

where N-1 is:

(gdb) print N-1
$50 = 18446744073709551615
(gdb)

so clearly a bug in allowing N to be zero.  Maybe something like:

   if (N <= 0)
 {
   GSL_ERROR("matrix size must be positive", GSL_EDOM);
 }
   else if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
 {
   GSL_ERROR("failed to allocate working space", GSL_ENOMEM);
 }
   else

and probably also do that kind of checking in the higher level function 
gsl_linalg_solve_symm_tridiag(), and probably many other parts of linear 
algebra.  This is a real project to examine the linear algebra package.

c. take cspline_init() and give it a check on size, insisting that it be >= 3.

Re: Test release for GSL 2.8

[Patrick Alken]

thank you for this, I upgraded autoconf to 2.72, which may fix the 
issue. I will prepare a new test release soon


On 5/11/24 13:05, Dave Allured - NOAA Affiliate via Users list for GNU 
Scientific Library (GSL) help wrote:

[External email - use caution]


MacPorts has a patch to 2.7.1 that seems related to the current dynamic
lookup warnings.  There is history and a trouble ticket behind the patch.
Perhaps this will help to fix this warning for 2.8.  See patch file in
today's MacPorts snapshot.
https://github.com/macports/macports-ports/tree/401b86c9f5c69025fb32cfbfd65413e235900ea7/math/gsl/files


On Sat, May 11, 2024 at 10:00 AM Dave Allured - NOAA Affiliate <
dave.allu...@noaa.gov> wrote:


Tested GSL 4.7.99 on Mac OS 12.7.4 Monterey, Intel mac.
Xcode 14.2 for compiler and developer tools.
C compiler is Apple clang 14.0.0.





The dynamic lookup and enum conversion warnings might need attention.





2 warning: -undefined dynamic_lookup may not work with chained fixups

Re: Test release for GSL 2.8

[Patrick Alken]

Hello, your suggestions for sparse matrices are good ones. I added a 
ticket to the bug tracker to track these - I will work on them when I 
have time


On 5/10/24 20:51, Brijesh Upadhaya wrote:

[External email - use caution]


Hi,

All tests passed for gsl-2.7.1+ in Debian 12 with GCC-12.2.0 and glibc-2.36.

Maybe you should update the author names in the source codes inside gsl/linalg. 
It seems that you have added/updated a few routines.

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.


Regards,
Brijesh


On 11. May 2024, at 1.18, Mohammad Akhlaghi  wrote:

Dear Patrick,

Thank you very much for sharing the good news about the new GSL release.

I just configured, built and checked the tarball on an Arch GNU/Linux with GCC 
14.1.1 and Glibc 2.39 (on 12 threads to speed it up):

./configure && make -j12 && make check -j12

There were no warnings or errors in the compilation or in the checks! 
Everything went smoothly all the way until the end :-).

Thanks for all the work on this wonderful package; GSL is a mandatory 
dependency of GNU Astronomy Utilities; and we rely heavily on it. I look 
forward to the release of GSL 2.8.

Cheers,
Mohammad

Test release for GSL 2.8

[Patrick Alken]

Dear all,

  It is time to make a new GSL release. I have uploaded a test release to:

ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz
ftp://alpha.gnu.org/gnu/gsl/gsl-2.7.99.tar.gz.sig

All reports are welcome - anyone who can test on various platforms would 
be appreciated (Linux, BSD, Mac OS, Windows).


Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need 
to install the python-based sphinx software first:


1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than 
system-wide)


3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: Nonlinear Least-Squares Fitting

[Patrick Alken]

Currently, only unconstrained nonlinear least squares is supported. You 
will need to look into other libraries for constrained methods


On 4/23/24 18:17, djvanpeur...@gmail.com wrote:

[External email - use caution]


Hello All,



I'm a new user and have a general question on the trust-region method.
There doesn't seem to be a direct way to restrict (bound) the solution space
to the non-negative domain.  Has anyone been able to do this by some
indirect mechanism?  I appreciate any suggestions.



Sincerely,

David

Re: basis splines for interpolation (not fitting)

[Patrick Alken]

code for this exists on the git repository, but not yet in an official 
release. If you can download from the git you can look at the examples


On 3/8/24 11:56, Simon Wiesheier wrote:

Dear all,

can I use basis splines for interpolation and not fitting?

Or, say I already computed the control points ai (f = sum_i B_i a_i) and
all what I want to do is
to evaluate the basis spline B_i, for which we need (i) a knot vector
(t_i,...,t_i+k) that describes the support of B_i and (ii) the order of the
spline k.

Are there classes / types that allow me to do that?

Best,
Simon

Re: Different eigenvectors obtained from gsl_eigen_hermv compared with scipy.linalg.eigh

[Patrick Alken]

Note that eigenvectors are unique only up to a constant factor. Since 
the output of hermv scales the eigenvectors to unit magnitude, they can 
be scaled by any complex number of unit magnitude and still be a valid 
eigenvector. It is likely that the GSL and scipy eigenvectors differ by 
a constant of unit magnitude.


On 1/3/24 23:39, Jiasen Guo wrote:

Dear all,
I am having an issue using gsl_eigen_hermv to reproduce a complex hermitian
eigenvalue problem previously solved with scipy.linalg.eigh. I have
included the cpp and python code in the attached zip file. Basically, the
matrix is created in cpp and solved using gsl. the same matrix (*real.csv*
for the real part and *imag.csv* for the imaginary part) is dumped out and
solved in Python again. While I get the same eigenvalues from the two
approaches, the eigenvalues are, however, different. And if I calculated
the module of the complex numbers in the eigenvector, then they are
identical from both methods.

I use g++ --std=c++11 test.cpp -o test.o -lgsl -lgslcblas for compilation.
Please help me understand what is the problem here. Thanks for your kind
help.

Best,
Jiasen

Re: question about the trust region in the gsl library.

[Patrick Alken]

I'm unclear what you are trying to do. Do you simply want to keep track 
of the trust radius at each step so you can plot it later?


On 5/2/23 07:29, Aleja Carmento wrote:

  Dear Sir/Madam,



I am currently working on a gsl least squares method. In this specific case
, i have it set to the levenberg-marquardt with geodesic acceleration. i
used the same example from the gsl website and added a few changes to it
like so :


#include #include #include
#include #include
#include #include
#include #include 


#define max_iterations 100


struct data
{
 double* t;
 double* y;
 size_t n;
};

static double scaled_norm(const gsl_vector* D, const gsl_vector* a)
{
 const size_t n = a->size;
 double e2 = 0.0;
 size_t i;

 for (i = 0; i < n; ++i) {
 double Di = gsl_vector_get(D, i);
 double ai = gsl_vector_get(a, i);
 double u = Di * ai;

 e2 += u * u;
 }

 return sqrt(e2);
}
/* model function: a * exp( -1/2 * [ (t - b) / c ]^2 )
*/doublegaussian(const double a, const double b, const double c, const
double t)
{
 const double z = (t - b) / c;
 return (a * exp(-0.5 * z * z));
}
intfunc_f(const gsl_vector* x, void* params, gsl_vector* f)
{
 struct data* d = (struct data*)params;
 double a = gsl_vector_get(x, 0);
 double b = gsl_vector_get(x, 1);
 double c = gsl_vector_get(x, 2);
 size_t i;

 for (i = 0; i < d->n; ++i)
 {
 double ti = d->t[i];
 double yi = d->y[i];
 double y = gaussian(a, b, c, ti);

 gsl_vector_set(f, i, yi - y);
 }

 return GSL_SUCCESS;
}
intfunc_df(const gsl_vector* x, void* params, gsl_matrix* J)
{
 struct data* d = (struct data*)params;
 double a = gsl_vector_get(x, 0);
 double b = gsl_vector_get(x, 1);
 double c = gsl_vector_get(x, 2);
 size_t i;

 for (i = 0; i < d->n; ++i)
 {
 double ti = d->t[i];
 double zi = (ti - b) / c;
 double ei = exp(-0.5 * zi * zi);

 gsl_matrix_set(J, i, 0, -ei);
 gsl_matrix_set(J, i, 1, -(a / c) * ei * zi);
 gsl_matrix_set(J, i, 2, -(a / c) * ei * zi * zi);
 }

 return GSL_SUCCESS;
}
intfunc_fvv(const gsl_vector* x, const gsl_vector* v,
 void* params, gsl_vector* fvv)
{
 struct data* d = (struct data*)params;
 double a = gsl_vector_get(x, 0);
 double b = gsl_vector_get(x, 1);
 double c = gsl_vector_get(x, 2);
 double va = gsl_vector_get(v, 0);
 double vb = gsl_vector_get(v, 1);
 double vc = gsl_vector_get(v, 2);
 size_t i;

 for (i = 0; i < d->n; ++i)
 {
 double ti = d->t[i];
 double zi = (ti - b) / c;
 double ei = exp(-0.5 * zi * zi);
 double Dab = -zi * ei / c;
 double Dac = -zi * zi * ei / c;
 double Dbb = a * ei / (c * c) * (1.0 - zi * zi);
 double Dbc = a * zi * ei / (c * c) * (2.0 - zi * zi);
 double Dcc = a * zi * zi * ei / (c * c) * (3.0 - zi * zi);
 double sum;

 sum = 2.0 * va * vb * Dab +
 2.0 * va * vc * Dac +
 vb * vb * Dbb +
 2.0 * vb * vc * Dbc +
 vc * vc * Dcc;

 gsl_vector_set(fvv, i, sum);
 }

 return GSL_SUCCESS;
}
void callback(const size_t iter, void* params,
 const gsl_multifit_nlinear_workspace* w)
{
 gsl_vector* f = gsl_multifit_nlinear_residual(w);
 gsl_vector* x = gsl_multifit_nlinear_position(w);
 double avratio = gsl_multifit_nlinear_avratio(w);
 double rcond;


 (void)params; /* not used */


 gsl_multifit_nlinear_trust_state* state =
(gsl_multifit_nlinear_trust_state*)w->state;
 /* compute reciprocal condition number of J(x) */
 gsl_multifit_nlinear_rcond(&rcond, w);

 fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, c = %.4f, |a|/|v|
= %.4f cond(J) = %8.4f, |f(x)| = %.4f\n",
 iter,
 gsl_vector_get(x, 0),
 gsl_vector_get(x, 1),
 gsl_vector_get(x, 2),
 avratio,
 1.0 / rcond,
 gsl_blas_dnrm2(f));
}
voidsolve_system(gsl_vector* x, gsl_multifit_nlinear_fdf* fdf,
 gsl_multifit_nlinear_parameters* params)
{
 const gsl_multifit_nlinear_type* T = gsl_multifit_nlinear_trust;
 const size_t max_iter = 200;
 const double xtol = 1.0e-8;
 const double gtol = 1.0e-8;
 const double ftol = 1.0e-8;

 const size_t n = fdf->n;
 const size_t p = fdf->p;

 size_t iter = 0;
 int status, info;

 gsl_multifit_nlinear_workspace* work =
gsl_multifit_nlinear_alloc(T, params, n, p);



 gsl_vector* f = gsl_multifit_nlinear_residual(work);
 gsl_vector* y = gsl_multifit_nlinear_position(work);


 double chisq0, chisq, rcond;


 double trust_radius_list[max_iterations];
 double step_size_list[max_iterations];

 /* initialize solver */
 gsl_multifit_nlinear_init(x, fdf, work);

 /* store initial cost */
 gsl_blas_ddot(f, f, &chisq0);

 /* iterate until convergence */
 if (callback) ca

Re: Reg :: Non Uniform FFT

[Patrick Alken]

Have you tried the NFFT library:

https://github.com/NFFT/nfft

That library is excellent. Also, I think this application is a bit too 
specialized (and complicated) to include in GSL, especially when there 
is a good GPL library available (NFFT)


Patrick

On 3/27/23 08:31, Sumit Adhikari wrote:

Hello Developers,

Recently, I am seeing a requirement for non uniform fft for the analysis of
our simulation results.

Is there any chance in future to add this function to gsl library. Or there
is some work around?

Regards,

Sumit

Re: fulfill the (cubic) basis splines the partition of unity at all points in [a,b]?

[Patrick Alken]

https://git.savannah.gnu.org/cgit/gsl.git

see also: https://savannah.gnu.org/git/?group=gsl

On 3/8/23 03:18, Simon Wiesheier wrote:
Is this <https://github.com/ampl/gsl/tree/master/bspline> the relevant 
section of the github repository you are referring to?


Am Mi., 8. März 2023 um 06:21 Uhr schrieb Patrick Alken 
:


Its also worth mentioning that there has been a substantial
overhaul to
the B-splines routines since v2.7. The new routines are on the git
repository, along with documentation. If you have the ability to
clone
the git and build the documentation from there I highly encourage it.
There are many example programs in the git documentation which are
not
in the 2.7 docs which may help you.

On 3/7/23 22:07, Rhys Ulerich wrote:
> This time remembering to CC the mailing list...
>
> On Tue, Mar 7, 2023, 9:57 PM Rhys Ulerich
 wrote:
>
>> On Tue, Mar 7, 2023, 5:11 PM Simon Wiesheier

>> wrote:
>>
>>> After reading the manual, it is not clear to me how GNU internally
>>> constructs the knot vector.
>>> There are the functions,
>>> gsl_bspline_knots
>>> gsl_bspline_knots_uniform,
>>> that create the knot vector based on given breakpoints.
>>>
>> I encourage you to initialize a cubic workspace (k=4, pick
nbreak) then to
>> use gsl_bspline_knots_uniform to have the GSL construct the
knot vector for
>> you given some [a, b]. You will be able to observe the
multiplicity of the
>> various knots in the resulting w->knots. The multiplicity is a
consequence
>> of the chosen k meaning that if you opted for quadratic or
quintic k you
>> will see a different knot multiplicity. Play around a bit.
>>
>> You may (or may not) find the routines at
>>
https://github.com/RhysU/suzerain/blob/master/suzerain/bspline.h to be
>> useful worked examples. Those include forming linear
combinations of the
>> basis.
>>
>> - Rhys
>>

Re: fulfill the (cubic) basis splines the partition of unity at all points in [a,b]?

[Patrick Alken]

Its also worth mentioning that there has been a substantial overhaul to 
the B-splines routines since v2.7. The new routines are on the git 
repository, along with documentation. If you have the ability to clone 
the git and build the documentation from there I highly encourage it. 
There are many example programs in the git documentation which are not 
in the 2.7 docs which may help you.


On 3/7/23 22:07, Rhys Ulerich wrote:

This time remembering to CC the mailing list...

On Tue, Mar 7, 2023, 9:57 PM Rhys Ulerich  wrote:


On Tue, Mar 7, 2023, 5:11 PM Simon Wiesheier 
wrote:


After reading the manual, it is not clear to me how GNU internally
constructs the knot vector.
There are the functions,
gsl_bspline_knots
gsl_bspline_knots_uniform,
that create the knot vector based on given breakpoints.


I encourage you to initialize a cubic workspace (k=4, pick nbreak) then to
use gsl_bspline_knots_uniform to have the GSL construct the knot vector for
you given some [a, b]. You will be able to observe the multiplicity of the
various knots in the resulting w->knots. The multiplicity is a consequence
of the chosen k meaning that if you opted for quadratic or quintic k you
will see a different knot multiplicity. Play around a bit.

You may (or may not) find the routines at
https://github.com/RhysU/suzerain/blob/master/suzerain/bspline.h to be
useful worked examples. Those include forming linear combinations of the
basis.

- Rhys

Re: Sparse submatrix view

[Patrick Alken]

Hello,

  For dense matrices, submatrix views work because all the matrix 
elements are stored in a nice consecutive, linear manner. This is not 
the case for sparse matrices. GSL additionally supports three different 
storage formats for sparse matrices (COO, CSR, CSC). So it is not 
trivial to implement submatrix views for sparse matrices.


You could simply construct separate sparse matrices for each subsystem, 
and then add them all together in the final step.


Patrick

On 12/9/22 01:00, Brijesh Upadhaya wrote:

Hi,

I am trying to solve a coupled electromagnetic field problem. I need to
combine several subsystems into one large system. Is there any idea to get
the submatrix view of the gsl_spmatrix?

It seems to me that the submatrix view is not yet implemented in the
development version.

BR,
Brijesh

Re: First ever GSL Technical Report (ALFs)

[Patrick Alken]

I like the idea of a DOI. Also the gnu.org website is pretty stable and 
for citing purposes, I think using the technical report aspect of bibtex 
with a link to the gnu.org site should be fine.


On 2/17/22 11:33, Mark Galassi wrote:

I suggest archiving on Zenodo https://zenodo.org/

Zenodo works well for this and gives you a citable location and a DOI; I put 
the gsl design document on there, so if you create a gsl community then we can 
have 2 things :-)

Re: First ever GSL Technical Report (ALFs)

[Patrick Alken]

Thanks Mark,

  I don't think it is suitable for publication in a journal. There is 
nothing novel here, its just technical details of how to calculate ALFs 
efficiently.


I'll see if I can fix the formatting issue on eq. 36. Its a good idea to 
cite GSL, I will add that :)


Patrick

On 2/16/22 22:51, Mark Galassi wrote:

Patrick, this is a great paper.  It shows the same care you apply to 
maintaining gsl.  The writing is also very clear, and I love the table of 
acronyms :-).  Do you plan to submit for publication in a numerical analysis 
journal, or submit to the arXiv?

The only previous quasi-report had been the ongoing design document that I kept 
going with James and Brian in the early years, but it is a working design doc, 
not anything of the scope you have shown here.

Tiny suggestions: alignment of equation 36 on page 5: "l >= 1" could be moved quite a bit 
to the left.  Maybe an extra & to make it match the start of the l >= 1, m > 0.

You might also want to also cite the reference manual (as we ask people to do 
when they use gsl :-) ).  A recent bibtex skeleton on that is the 2019 Network 
Theory edition.

@book{gslteam2019gnuscientificlibrary,
   title={GNU Scientific Library Reference Manual},
   author={Galassi, Mark and Davies, Jim and Theiler, James and Gough, Brian 
and Jungman, Gerard and Alken, Patrick and Booth, Michael and Rossi, Fabrice 
and Ulerich, Rhys},
   year={2019},
   isbn={9780954612078},
   publisher={Network Theory Limited}
}

First ever GSL Technical Report (ALFs)

[Patrick Alken]

Hello all,

  I have recently re-written the array versions of the Associated 
Legendre functions (gsl_sf_legendre) to be much more efficient. The 
factors in the recurrence relations are now precomputed, and they are 
stored in a memory layout to maximize cache efficiency. The result is a 
significant improvement in speed over the previous versions. I decided 
to write up the implementation details in a technical report, which is 
now available at:


https://www.gnu.org/software/gsl/tr/tr001.pdf

This is effectively the very first GSL Technical Report. I envision 
these to be similar to the LAPACK working notes...useful documents to 
explain the details of complex algorithms in GSL for future developers 
and even users. Any comments on the report are welcome. I am labeling it 
as a draft for now, as I may update it in the future.


Best regards,

Patrick

Re: RANLUX++

[Patrick Alken]

Hello,

  Yes please do submit a patch. It would be best to submit a git diff 
against the latest master branch of GSL. I'll take a look when I can 
find the time. Thanks for your efforts!


Patrick

On 8/11/21 8:49 AM, Jonas Hahnfeld wrote:

Hello GSL developers,

over the past months, we (ROOT developers from CERN) created a portable 
implementation of RANLUX++, which was first proposed by Alexei Sibidanov in 
2017 and uses an LCG equivalent to achieve a much higher luxury level than the 
original Ranlux and provides some more advantages (better initialization and 
fast skipping). Results from this portable implementation are available in 
https://arxiv.org/abs/2106.02504 and already shipping with the latest release 
of ROOT.

For wider distribution and ease of use, we would like to also contribute this 
generator to GSL and have already written the necessary plumbing code to 
integrate with the GSL interfaces. It is currently available at my GitHub: 
https://github.com/hahnjo/ranluxpp

Could you let us know what the next steps would be? Shall we submit a patch for 
integration into GSL or do we have to be added to the list of 
"Extensions/Applications" first?

Kind regards
Jonas

Re: Problem making anonymous contribution

[Patrick Alken]

Hello,

  I don't know why a git diff would try to push changes upstream. Did 
you clone an anonymous copy of the gsl git repository? Perhaps you will 
need to google that error to try to understand what is going wrong.


Patrick

On 6/12/21 1:32 PM, Carlos Davalillo wrote:

Hello

I recently got in contact with Mr. Patrick from the GSL group to contribute
an alternative quadrature algorithm based on the tanh-sinh method. When I
try to send git diff I get the following error message

  Unable to push to 'origin' - unexpected http status code: 403

I tried to do it by the ssh way but got the same result.

So I ask you for some help to overcome this problem. I'm located in
Cabimas, Venezuela by the way.

Thank you in advance for your help!

Re: Copying C array to a GSL vector

[Patrick Alken]

First you need to make a "vector view" of your array, and then you can 
treat it as a vector. Example:


gsl_vector_view v = gsl_vector_view_array(my_array, n);

gsl_vector_memcpy(dest, &v.vector);

On 6/10/21 11:48 PM, Fritz Sonnichsen wrote:

I am about to use some filtering routines and I am new to GSL. My C routine
creates a vector (1x1024) and I want to convert this to a gsl_vector type.
(In particular I want to use the /int gsl_vector_memcpy(gsl_vector *dest,
constgsl_vector *src) function.

I see copying routines in your document for copying gsl_vector types
 (for example int gsl_vector_memcpy(gsl_vector *dest, constgsl_vector
*src) )
  but I do not see one to copy my C vector to the gsl type for further
processing. Is there such a routine?

Thanks
Fritz

GNU Scientific Library 2.7 released

[Patrick Alken]

Version 2.7 of the GNU Scientific Library (GSL) is now available. GSL
provides a large collection of routines for numerical computing in C.

This release introduces some new features and fixes several bugs. The
full NEWS file entry is appended below.

The file details for this release are:

ftp://ftp.gnu.org/gnu/gsl/gsl-2.7.tar.gz
ftp://ftp.gnu.org/gnu/gsl/gsl-2.7.tar.gz.sig

The GSL project homepage is http://www.gnu.org/software/gsl/

GSL is free software distributed under the GNU General Public License.

Thanks to everyone who reported bugs and contributed improvements.

Patrick Alken

---

* What is new in gsl-2.7:

|** fixed doc bug for gsl_histogram_min_bin (lhcsky at 163.com) ** fixed
bug #60335 (spmatrix test failure, J. Lamb) ** fixed bug #36577 **
clarified documentation on interpolation accelerators (V. Krishnan) **
fixed bug #45521 (erroneous GSL_ERROR_NULL in ode-initval2, thanks to M.
Sitte) ** fixed doc bug #59758 ** fixed bug #58202 (rstat median for
n=5) ** added support for native C complex number types in gsl_complex
when using a C11 compiler ** upgraded to autoconf 2.71, automake 1.16.3,
libtool 2.4.6 ** updated exponential fitting example for nonlinear least
squares ** added banded LU decomposition and solver (gsl_linalg_LU_band)
** New functions added to the library: - gsl_matrix_norm1 -
gsl_spmatrix_norm1 - gsl_matrix_complex_conjtrans_memcpy -
gsl_linalg_QL: decomp, unpack - gsl_linalg_complex_QR_* (thanks to
Christian Krueger) - gsl_vector_sum - gsl_matrix_scale_rows -
gsl_matrix_scale_columns - gsl_multilarge_linear_matrix_ptr -
gsl_multilarge_linear_rhs_ptr - gsl_spmatrix_dense_add (renamed from
gsl_spmatrix_add_to_dense) - gsl_spmatrix_dense_sub -
gsl_linalg_cholesky_band: solvem, svxm, scale, scale_apply -
gsl_linalg_QR_UD: decomp, lssolve - gsl_linalg_QR_UU: decomp, lssolve,
QTvec - gsl_linalg_QR_UZ: decomp - gsl_multifit_linear_lcurvature -
gsl_spline2d_eval_extrap ** bug fix in checking vector lengths in
gsl_vector_memcpy (dieg...@pm.me) ** made gsl_sf_legendre_array_index()
inline and documented gsl_sf_legendre_nlm()|

Re: Test release for GSL 2.7

[Patrick Alken]

Thanks Patrick, some responses below

On 5/29/21 4:50 AM, Patrick Dupre wrote:
> I run the test and generate the doc on
> Linux Teucidide 5.11.22-100.fc32.x86_64 #1 SMP Wed May 19 18:58:25 UTC 2021 
> x86_64 x86_64 x86_64 GNU/Linux
>
> I did not any issue.
>
> Questions
>
> gsl is based at least partially on blas (and lapack?)
GSL calls an external CBLAS library. It ships with a default (gslcblas)
library, but it is far better to use an optimized library such as ATLAS
which is parallelized
> Can we benefit of the parallelism (openmp)?
> If yes how?
> Did you test?
>
> The default c compiler is gcc.
> Did you test gsl with clang?
I have not used clang
> It something offers very interesting performances compared with gcc.
>
> Best Regards.
>
> ===
>  Patrick DUPRÉ | | email: pdu...@gmx.com
>  Laboratoire interdisciplinaire Carnot de Bourgogne
>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>  Tel: +33 (0)380395988| | Room# D114A
> =======
>
>
>> Sent: Friday, May 28, 2021 at 9:26 PM
>> From: "Patrick Alken" 
>> To: "help-gsl@gnu.org" 
>> Subject: Test release for GSL 2.7
>>
>> Dear all,
>>
>>   It is time to make a new GSL release. I have uploaded a test release to:
>>
>> ftp://alpha.gnu.org/gnu/gsl/gsl-2.6.90.tar.gz
>> ftp://alpha.gnu.org/gnu/gsl/gsl-2.6.90.tar.gz.sig
>>
>> All reports are welcome - anyone who can test on various platforms would be 
>> appreciated (Linux, BSD, Mac OS, Windows).
>>
>> Please try testing the build:
>>
>> ./configure && make && make check
>>
>> If you wish you can test building the documentation, but you will need to 
>> install the python-based sphinx software first:
>>
>> 1. pip install -U --user Sphinx
>>
>> 2. pip install -U --user sphinx_rtd_theme
>>
>> (the --user flag will install in your home directory rather than system-wide)
>>
>> 3. cd doc ; make html
>>
>> Please report any successes/failures.
>>
>> Thanks,
>> Patrick
>>
>>
>>

Test release for GSL 2.7

[Patrick Alken]

Dear all,

  It is time to make a new GSL release. I have uploaded a test release to:

ftp://alpha.gnu.org/gnu/gsl/gsl-2.6.90.tar.gz
ftp://alpha.gnu.org/gnu/gsl/gsl-2.6.90.tar.gz.sig

All reports are welcome - anyone who can test on various platforms would be 
appreciated (Linux, BSD, Mac OS, Windows).

Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need to 
install the python-based sphinx software first:

1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than system-wide)

3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: savitzky golay routines?

[Patrick Alken]

Currently this is not implemented in GSL

On 5/27/21 7:36 PM, Fritz Sonnichsen wrote:
> Looking over the routines available from GSL I don't see a lot about
> filters. There are some areas of curve fits covered.
>  I am working with Savitzky-Golay filters at this time--is there no
> routines for this type of filtering/fitting?
>
> Thanks
> Fritz

Re: physical constants that might need updating

[Patrick Alken]

I logged this email into the bug tracker so it won't be lost. I agree it
would be best to stay current with the latest values of the constants. I
am concerned that such a change could potentially break existing
software which might depend on certain values.

Patrick

On 5/5/21 11:17 PM, Mark Galassi wrote:
> 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 */
>

Re: QR and Pivoted QR for Complex Matrix

[Patrick Alken]

Hello,

  The complex QR decomposition is available in the git repository - we
haven't made a release since it was added. The algorithm is the Level-3
BLAS method of Elmroth and Gustavson, so you need to supply the T matrix
(instead of the tau vector) to the routine.

If you are willing/able to clone the git and install it, then the
function you will need is:

gsl_linalg_complex_QR_decomp_r(gsl_matrix_complex * A,
gsl_matrix_complex * T)

The other functions (solve, lssolve) follow the same calling format as
the real versions (with the _r suffix). Pivoted QR is only available for
real matrices at the moment.

On 4/13/21 3:16 PM, Songting Luo wrote:
> Hello GSL users and developers,
>
> Could you please let me know whether GSL has QR and pivoted QR for complex
> matrices? I couldn't find them in the online manual.  Thanks.
>
> -SL

Re: [PATCH] clarify accelerator usage for multiple interpolations

[Patrick Alken]

Thank you, I have added this text to the documentation

On 1/29/21 11:26 PM, vishn...@gmail.com wrote:
> From: "Vishnu V. Krishnan" 
>
> ---
>  doc/interp.rst | 10 ++
>  1 file changed, 6 insertions(+), 4 deletions(-)
>
> diff --git a/doc/interp.rst b/doc/interp.rst
> index 5a062d50..9ee31e09 100644
> --- a/doc/interp.rst
> +++ b/doc/interp.rst
> @@ -168,10 +168,12 @@ which is a kind of iterator for interpolation lookups.
>  
>  .. function:: gsl_interp_accel * gsl_interp_accel_alloc (void)
>  
> -   This function returns a pointer to an accelerator object, which is a
> -   kind of iterator for interpolation lookups.  It tracks the state of
> -   lookups, thus allowing for application of various acceleration
> -   strategies.
> +   This function returns a pointer to an accelerator object, which is a kind 
> of
> +   iterator for interpolation lookups.  It tracks the state of lookups, thus
> +   allowing for application of various acceleration strategies. When 
> multiple,
> +   concurrent interpolants are in use, the same accelerator object may be 
> used
> +   for all functions defined on the same domain (x_array), but different
> +   accelerators must be defined for data defined on different domains.
>  
>  .. function:: size_t gsl_interp_accel_find (gsl_interp_accel * a, const 
> double x_array[], size_t size, double x)
>  

Re: Using matrix decompositions

[Patrick Alken]

Hello,

  There are 2 interfaces for the QR decomposition, the original
interface which uses Level-2 BLAS, and a recently added interface using
Level-3 BLAS (_r suffix to function names). In both cases, there is an
"unpack" function to extract the Q matrix:

gsl_linalg_QR_unpack
gsl_linalg_QR_unpack_r

The R matrix is always stored in the upper triangle of the input matrix
(A) by the QR_decomp routines.

In most applications, you never need to actually construct the full Q
matrix itself. If you describe more about what you want to do, or post
your source code, we can probably help more.

Thanks,
Patrick

On 1/30/21 1:49 AM, Iaroslav Postovalov wrote:
> Hello,
>
> I'm quite confused with using QR decomposition with GSL. The library
> returns data in very unclear data format, and I don't understand how to get
> Q and R matrices from it.
>
> Certain users may have the same problem for other decompositions (since
> literally all of them require some effort to extract the resulting
> matrices), and I think that the GSL documentation has to provide some
> instructions for them.
>
> All the bests,
> Iaroslav

Re: Interpolation accelerator usage

[Patrick Alken]

The accelerator works by doing a binary search on the x array, and then
caches the index for future lookups. So you can indeed use it for
multiple data sets if they share a common x axis.

Patrick

On 1/28/21 9:54 AM, Vishnu V. Krishnan wrote:
> Hi!
>
> I have a question about allocating accelerators for interpolation look-ups.
>
> If I have multiple data-sets for which I need interpolating functions, do I 
> allocate separate accelerators for each of them? Or can all the data sets 
> that 
> share the same x-axes, share accelerators?
>
> If it is the latter case, it would be nice if this were made clearer in the 
> documentation: https://www.gnu.org/software/gsl/doc/html/interp.html
>
> Thanks a lot!
> Vishnu
>

Re: POTENTIAL SPAM: Speed of complex vs gsl_complex computations

[Patrick Alken]

Pavel,

  GSL was originally written to conform to the C89 standard, which did
not support native complex numbers, so GSL wrote its own implementation
of complex numbers and arithmetic. Native support for complex numbers
was introduced in C99 and later extended in C11. In the current master
branch of the git repository, GSL now supports native C complex numbers.

If you #include  prior to gsl_complex.h in your C program,
and you are using a C11 compliant compiler, GSL will automatically
define gsl_complex to be the same as complex double. This functionality
is backward compatible, since the original definition is binary
compatible with the native C type.

This means you can do any operation with gsl_complex as you can with
complex double. You are welcome to try it out and see if it does what
you need.

Patrick

On 1/11/21 2:16 PM, Pavel via Users list for GNU Scientific Library
(GSL) help wrote:
> Hello!
>
> I've installed GSL for use in Visual Studio Community 2017 following
> instructions of the blog
> https://secure-web.cisco.com/1RG6v1TMIMAP9F4BNHb-l2ekF6CbMsLcORRPqxI-nfy4wYhY-jPm9hb53F50xGHIwAi7jy-FQz1OHyLzoGzzlJoblX_BHz5y8IALUDZW1tIW6jd7jSFMcKZ0RGjursfcehDgambN9crgIlKwo2lv0HC9X_o5WZY9JAnw_DhuuZegvZgv5-tk3R2igAqUyB2ED-rAfiITcxYuszOSRfeoQ0gNli9NjdeCLR83rhEfaEYlCLGphD4WSDwaIMD7pxA1S3z_ax50iPAuFRx7ebxKedimatr2j6FNiLhsJ-AVzAdRRvoqZx3HePGeIYIVhN5XFzLM0rv3XSO2ebKCIO6fBUmLqAfZnZNaQ4_yF_Q31zMnQhP8vQxFmBpP0gG1qvrMzMPodWW5gF1AGwSVTtTOcsO01r0ZxPZMRM1kfUISHCj59oJh06sD1FUk6S1-eTv8JTDyvPRV_HlmgOTI44z4Hgw/https%3A%2F%2Fsolarianprogrammer.com%2F2020%2F01%2F26%2Fgetting-started-gsl-gnu-scientific
> -library-windows-macos-linux/.
>
> I've decided to measure speed of calculations of internal complex numbers
> realization using template complex (including header complex.h) and
> realization of complex numbers in GSL library. For this purpose I've
> calculated large vectors ( length 1000) of some complex function (that
> is actually frequency dependence of graphene conductivity). I've used speed
> optimization /O2 and /DHAVE_INLINE compiler options. The calculation time of
> the vector using gsl_complex numbers was 7 s, while using internal
> complex numbers took only 1s. Is it true that GSL complex numbers
> work slower than internal realization or maybe I've missed some optimization
> settings of GSL?
>
>

Re: eigensystem

[Patrick Alken]

It is roughly equivalent. GSL nonsymmv will return an error
(GSL_EMAXITER) if the algorithm fails to converge on all eigenvalues.
The number of converged eigenvalues is stored in
w->nonsymm_workspace_p->n_evals

Patrick

On 1/20/21 1:36 AM, Patrick Dupre wrote:
> Thank you Patrick,
>
> It seems that dgeev return INFO.
>  = 0: successful exit
> < 0:  if INFO = -i, the i-th argument had an illegal value.
> > 0:  if INFO = i, the QR algorithm failed to compute all the
> eigenvalues, and no eigenvectors have been computed; elements i+1:N
> of WR and WI contain eigenvalues which have converged.
>
> It is the same value returned by gsl_eigen_nonsymmv?
>
>
>> The GSL algorithm does not check if the matrix is diagonalizable, so it
>> assumes that the input matrix is. This is the same behavior as LAPACK
>> DGEEV - in fact I ran your matrix through lapack and found the same answer.
>>
>> In order to check if your matrix is diagonalizable, one strategy is to
>> compute the condition number of the eigenvector matrix V (after the
>> nonsymmv computation). If the condition number is large (say > 1e7) then
>> that is a bad sign.
>>
>> I don't think there is an explicit routine in GSL atm to compute the
>> condition number of a complex matrix, but a quick way to do it would be
>> to compute the LU decomposition of V, and then look at the diagonal
>> elements of the U factor. If any of them are tiny, its badly conditioned.
>>
>> I can add a function to compute a more accurate estimate of the
>> condition number if needed.
>>
>> Patrick
>>
>> On 1/19/21 1:51 PM, Patrick Dupre wrote:
>>> gsl_eigen_nonsymmv_workspace
>>> has no member n_evals
>>>
>>> issue:
>>>
>>> Diagonalizing
>>> double data_3 [] = { 0.0, 0.0, 1.0,
>>> 0.0, 0.0, 0.0,
>>> 0.0, 0.0, 0.0 } ;
>>>
>>> I get
>>> eigenvalue = 0 +0i
>>> eigenvector = 
>>> 1 +0i
>>> 0 +0i
>>> 0 +0i
>>> eigenvalue = 0 +0i
>>> eigenvector = 
>>> 0 +0i
>>> 1 +0i
>>> 0 +0i
>>> eigenvalue = 0 +0i
>>> eigenvector = 
>>> -1 +0i
>>> 0 +0i
>>> 3.00625e-292 +0i
>>>
>>>
>>> which is wrong.
>>> The last eigenvector is not correct because this matrix is not 
>>> diagonalizable.
>>>
>>> I need to identify such matrices.
>>>
>>>
>>> ===
>>>  Patrick DUPRÉ | | email: pdu...@gmx.com
>>>  Laboratoire interdisciplinaire Carnot de Bourgogne
>>>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>>>  Tel: +33 (0)380395988
>>> ===
>>>
>>>
>>>> Sent: Tuesday, January 19, 2021 at 6:56 PM
>>>> From: "Patrick Alken" 
>>>> To: help-gsl@gnu.org
>>>> Subject: Re: eigensystem
>>>>
>>>> What do you mean by handle it? According to the documentation, if the
>>>> function cannot compute all eigenvalues, an error code is returned. In
>>>> the case of gsl_eigen_nonsymm, the number of converged eigenvalues is
>>>> stored in w->n_evals.
>>>>
>>>> Patrick
>>>>
>>>> On 1/19/21 10:33 AM, Patrick Dupre wrote:
>>>>> Hello,
>>>>>
>>>>> Is there a way to handle the possible error of gsl_eigen_nonsymmv ?
>>>>>
>>>>> For example, when the matrix is not diagonalizable.
>>>>>
>>>>> Thanks
>>>>>
>>>>> ===
>>>>>  Patrick DUPRÉ | | email: pdu...@gmx.com
>>>>>  Laboratoire interdisciplinaire Carnot de Bourgogne
>>>>>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>>>>>  Tel: +33 (0)380395988
>>>>> ===
>>>>>
>>>>>
>>>>
>>

Re: eigensystem

[Patrick Alken]

The GSL algorithm does not check if the matrix is diagonalizable, so it
assumes that the input matrix is. This is the same behavior as LAPACK
DGEEV - in fact I ran your matrix through lapack and found the same answer.

In order to check if your matrix is diagonalizable, one strategy is to
compute the condition number of the eigenvector matrix V (after the
nonsymmv computation). If the condition number is large (say > 1e7) then
that is a bad sign.

I don't think there is an explicit routine in GSL atm to compute the
condition number of a complex matrix, but a quick way to do it would be
to compute the LU decomposition of V, and then look at the diagonal
elements of the U factor. If any of them are tiny, its badly conditioned.

I can add a function to compute a more accurate estimate of the
condition number if needed.

Patrick

On 1/19/21 1:51 PM, Patrick Dupre wrote:
> gsl_eigen_nonsymmv_workspace
> has no member n_evals
>
> issue:
>
> Diagonalizing
> double data_3 [] = { 0.0, 0.0, 1.0,
>   0.0, 0.0, 0.0,
>   0.0, 0.0, 0.0 } ;
>
> I get
> eigenvalue = 0 +0i
> eigenvector = 
> 1 +0i
> 0 +0i
> 0 +0i
> eigenvalue = 0 +0i
> eigenvector = 
> 0 +0i
> 1 +0i
> 0 +0i
> eigenvalue = 0 +0i
> eigenvector = 
> -1 +0i
> 0 +0i
> 3.00625e-292 +0i
>
>
> which is wrong.
> The last eigenvector is not correct because this matrix is not diagonalizable.
>
> I need to identify such matrices.
>
>
> ===
>  Patrick DUPRÉ | | email: pdu...@gmx.com
>  Laboratoire interdisciplinaire Carnot de Bourgogne
>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>  Tel: +33 (0)380395988
> =======
>
>
>> Sent: Tuesday, January 19, 2021 at 6:56 PM
>> From: "Patrick Alken" 
>> To: help-gsl@gnu.org
>> Subject: Re: eigensystem
>>
>> What do you mean by handle it? According to the documentation, if the
>> function cannot compute all eigenvalues, an error code is returned. In
>> the case of gsl_eigen_nonsymm, the number of converged eigenvalues is
>> stored in w->n_evals.
>>
>> Patrick
>>
>> On 1/19/21 10:33 AM, Patrick Dupre wrote:
>>> Hello,
>>>
>>> Is there a way to handle the possible error of gsl_eigen_nonsymmv ?
>>>
>>> For example, when the matrix is not diagonalizable.
>>>
>>> Thanks
>>>
>>> ===
>>>  Patrick DUPRÉ | | email: pdu...@gmx.com
>>>  Laboratoire interdisciplinaire Carnot de Bourgogne
>>>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>>>  Tel: +33 (0)380395988
>>> ===
>>>
>>>
>>
>>

Re: eigensystem

[Patrick Alken]

A singular matrix can still be diagonalizable. Take

[ 1 0 ]
[ 0 0 ]

for example. So it wouldn't make sense to check if a matrix is singular.
The best test is to compute the condition number of V as per my previous
email.

Patrick

On 1/19/21 3:30 PM, Patrick Dupre wrote:
> Thanks, but I think that gsl_eigen_nonsymmv must return an error or a warming
> when dealing with singular matrices.
> Of course, I just took a 3 x 3 matrix as an example, but I need to deal with
> any size of matrices, and some maybe singular.
>
> Furthermore, maxima detects the singular matrices. Why not GSL?
> Currently, we can "diagonalize" singular matrices, and deal with eigenvectors
> as they were relevant eigenvectors.
>
> Indeed, gsl_eigen_nonsymmv  could return at least an error if the matrix 
> determinant is zero
> or if eigenvalues are degenerated.
> I guess that the first thing that gsl_eigen_nonsymmv does is to calculate the 
> determinant.
>
>> Does this help?
>>
>> https://lists.gnu.org/archive/html/help-gsl/2005-02/msg00012.html
>>
>> If you get an error or a zero, the matrix cannot be inverted.
>>
>> Maybe you are asking "Isn't there a simpler way?" and I think the answer
>> is "No." I can look at your matrix and see that it does not have full
>> rank because columns 0 and 1 are perfectly correlated. But the condition
>> that must apply to have full rank is more complex: No linear combination
>> of columns can be perfectly correlated with any other linear
>> combination. That's not hard to test with a 3x3 matrix, but it becomes
>> impractical with a larger matrix.
>>
>> You could compute correlations between all pairs of columns and if any
>> correlation is 1.0 then the matrix does not have full rank (and if it's
>> very close to 1.0, then you will have stability issues). But there are
>> matricies that will pass that test that do not have full rank. So that
>> test is probably somewhat cheaper computationally, but has a known false
>> negative rate and is not a good test.
>>
>> -Alan
>>
>>
>> On 1/19/2021 4:05 PM, Patrick Dupre wrote:
>>> This is perfectly correct.
>>>
>>> My question is how do I detect singular matrix since
>>> gsl_eigen_nonsymmv
>>> does not do it?
>>>
>>>
>>>> Your matrix is singular (the first two columns are perfectly correlated,
>>>> so your matrix does not have full rank).
>>>>
>>>> The standard test of singularity is that the determinant cannot be
>>>> computed (you try to calculate it, and you get zero or a numerical
>>>> error). This post describes the problem and a shortcut:
>>>>
>>>> https://stackoverflow.com/questions/13145948/how-to-find-out-if-a-matrix-is-singular
>>>>
>>>> -Alan
>>>>
>>>>
>>>> On 1/19/2021 2:51 PM, Patrick Dupre wrote:
>>>>> gsl_eigen_nonsymmv_workspace
>>>>> has no member n_evals
>>>>>
>>>>> issue:
>>>>>
>>>>> Diagonalizing
>>>>> double data_3 [] = { 0.0, 0.0, 1.0,
>>>>>   0.0, 0.0, 0.0,
>>>>>   0.0, 0.0, 0.0 } ;
>>>>>
>>>>> I get
>>>>> eigenvalue = 0 +0i
>>>>> eigenvector = 
>>>>> 1 +0i
>>>>> 0 +0i
>>>>> 0 +0i
>>>>> eigenvalue = 0 +0i
>>>>> eigenvector = 
>>>>> 0 +0i
>>>>> 1 +0i
>>>>> 0 +0i
>>>>> eigenvalue = 0 +0i
>>>>> eigenvector = 
>>>>> -1 +0i
>>>>> 0 +0i
>>>>> 3.00625e-292 +0i
>>>>>
>>>>>
>>>>> which is wrong.
>>>>> The last eigenvector is not correct because this matrix is not 
>>>>> diagonalizable.
>>>>>
>>>>> I need to identify such matrices.
>>>>>
>>>>>
>>>>> ===
>>>>>  Patrick DUPRÉ | | email: pdu...@gmx.com
>>>>>  Laboratoire interdisciplinaire Carnot de Bourgogne
>>>>>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>>>>>  Tel: +33 (0)380395988
>>>>> ===
>>>>>
>>>>>
>>>>>> Sent: Tuesday, January 19, 2021 at 6:56 PM
>>>>>> 

Re: eigensystem

[Patrick Alken]

What do you mean by handle it? According to the documentation, if the
function cannot compute all eigenvalues, an error code is returned. In
the case of gsl_eigen_nonsymm, the number of converged eigenvalues is
stored in w->n_evals.

Patrick

On 1/19/21 10:33 AM, Patrick Dupre wrote:
> Hello,
>
> Is there a way to handle the possible error of gsl_eigen_nonsymmv ?
>
> For example, when the matrix is not diagonalizable.
>
> Thanks
>
> ===
>  Patrick DUPRÉ | | email: pdu...@gmx.com
>  Laboratoire interdisciplinaire Carnot de Bourgogne
>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>  Tel: +33 (0)380395988
> ===
>
>

Re: GSL Commercial License

[Patrick Alken]

Hello,

  I don't fully understand whether you want to use GSL in a
closed-source app? The answer to that is no, as GSL is released under
the GPL.

If your app has a GPL-compatible license, then you can use GSL as you wish.

Best,
Patrick

On 1/13/21 1:26 PM, Haris Ali wrote:
> Hi,
>
> I was curious if it's possible to purchase a commercial license to use GSL
> in a public application in the app store. If so, how much would it cost?
>
> Thank you,

New branch for native complex support

[Patrick Alken]

Hi all,

  FYI there is a new branch on the GSL git repository called 'cplx'. In
this branch you will find that you can use gsl_complex as a native C99
complex double type, if you #include  prior to including
. This is possible because the gsl_complex type is
binary compatible with complex double. So basically you can
add/multiply/etc gsl_complex variables without needing to call the
gsl_complex_math routines now.

Many people have requested something like this over the years, so I am
finally getting around to doing it. Right now, this branch is
experimental, I am still working on writing some unit tests, and also
documentation. But I have been using it myself for a while and things
are working well.

The main issue right now, is you cannot do things like:

gsl_complex z;
GSL_REAL(z) = a;

which you can do with the original implementation. I don't know how to
define a macro which can do this with the C99 standard, although there
are non-portable GNU extensions which can do this. But otherwise,
everything should work as expected.

If interested, please try out the new branch and report back with
successes or failures.

Best,

Patrick

Re: How to contribute

[Patrick Alken]

Hello,

  The easiest way would be to clone the git repository (see GSL
webpage), make your additions, and then send me a diff (git diff) file.
It would be highly useful to also document your additions (see the doc
directory).

Thanks,
Patrick

On 12/14/20 10:34 AM, Tim wrote:
> Hello,
> I would like to make a small contribution to gsl (discrete hankel
> transform algorithm).
> How to file a pull request for a small contribution to gsl ?
>
> kind regards,
>
> Tim Berberich
>
>

Re: Sum algorithm

[Patrick Alken]

Hello,

  GSL uses the Welford algorithm, which is based on recurrence relations
for mean and variance. See

https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford%27s_online_algorithm

Source code here:

https://git.savannah.gnu.org/cgit/gsl.git/tree/statistics/mean_source.c
https://git.savannah.gnu.org/cgit/gsl.git/tree/statistics/variance_source.c

Reference:

Welford, B. P. (1962). "Note on a method for calculating corrected sums
of squares and products". /Technometrics
/. *4* (3): 419–420. doi
:10.2307/1266577
. JSTOR
 1266577
.

Patrick

On 12/9/20 5:23 PM, Dante Doménech wrote:
> Hello.
> It's a fast question. When the library calculates mean value, variance and
> others, in order to sum, does it use the kahan ( if yes, normal or second
> order ), or pairwise summation or something else?
> Thanks for your help.

Re: Memory limit allocation

[Patrick Alken]

Hi Pablo,

  The function gsl_spmatrix_alloc needs to allocate space for the
non-zero entries of the matrix. By default, it assumes a density of 10%.
So if N = 100,000, then it will try to allocate N*N*0.1 = 1e9 elements,
each of size sizeof(double), the the total allocation will be 8GB. But
then in the triplet representation, it needs not only the data values,
but also array tracking the row and column indices, so it will need
another 2*1e9*sizeof(int) = 8GB. So about 16GB total.

If you don't have 16GB of ram available, and if your matrix has much
less than 10% density, then you can instead use the function:

gsl_spmatrix_alloc_nzmax(n, n, nzmax, GSL_SPMATRIX_TRIPLET)

which allows you to precisely state how many non-zero entries will be in
the matrix. Then you can reduce the amount of memory allocated.

Patrick

On 12/9/20 6:22 AM, Pablo wrote:
> Hi,
>
> I've been searching in the web to solve this problem but I haven't
> found any solutions. My problem is related to the allocation limit of
> a program using sparse matrices. My project needs very large sparse
> matrices, with dimensions up to, i.e., 100,000x100,000, and the
> program returns;
>
> gsl: init_source.c:389: ERROR: failed to allocate space for memory block
> Default GSL error handler invoked.
>
> I've read something like memory leak in a loop, but my code can't be
> more simple;
>
> gsl_spmatrix* m = gsl_spmatrix_alloc(10, 10);
>
> How could I remove the limitation that prevents me for allocating
> space for a large matrix such like that? With dimensions 10,000x10,000
> it still works, and with other libraries such as Eigen3, I'm also able
> to build large matrices.
>
> Pablo
>
>

Re: Benchmarking GSL

[Patrick Alken]

Ludovic,

  I'm not aware of anyone doing formal benchmarking of GSL routines. I
myself have done some simple comparisons between GSL and LAPACK for some
of the matrix factorizations, but I don't have any formal results to
provide you. But I would be interested to see your results.

Patrick

On 11/30/20 9:18 AM, Ludovic Courtès wrote:
> Hi Patrick,
>
> Patrick Alken  skribis:
>
>>   Is there a particular part of the GSL library you want to benchmark?
>> Some routines have been optimized to run as fast as possible while
>> others have not. The overall focus of GSL algorithms is on robustness
>> and reliability, not necessarily speed. So you will probably find quite
>> different results depending which algorithms you are benchmarking.
> I’m interested in benchmarking routines that benefit from vectorization
> (AVX2, AVX512, etc.).
>
> Based on GCC’s ‘-fopt-info-vec’, this could be any of the following
> components:
>
> --8<---cut here---start->8---
> bspline
> bst
> cdf
> cheb
> combination
> complex
> dht
> diff
> eigen
> fft
> filter
> fit
> histogram
> integration
> interpolation
> linalg
> matrix
> min
> monte
> movstat
> multifit
> multifit_nlinear
> multifit
> multilarge_nlinear
> multilarge
> multimin
> multiroots
> multiset
> ode-initval2
> ode-initval
> permutation
> poly
> qrng
> randist
> rng
> roots
> rstat
> sort
> spblas
> specfunc
> splinalg
> spmatrix
> statistics
> vector
> wavelet
> --8<---cut here---end--->8---
>
> This is quite a long list, but really benchmarks for a couple of these
> would already be helpful to me.
>
> For background, I’m currently looking into automatic “function
> multi-versioning” where code is generated for several ISA extensions and
> the right one is selected at load time¹.
>
> Thanks in advance!
>
> Ludo’.
>
> PS: Please keep me Cc’d as I’m not subscribed.
>
> ¹ https://gitlab.inria.fr/guix-hpc/function-multi-versioning

Re: Benchmarking GSL

[Patrick Alken]

Hello,

  Is there a particular part of the GSL library you want to benchmark?
Some routines have been optimized to run as fast as possible while
others have not. The overall focus of GSL algorithms is on robustness
and reliability, not necessarily speed. So you will probably find quite
different results depending which algorithms you are benchmarking.

Patrick


On 11/24/20 3:24 AM, Ludovic Courtès wrote:
> Hello,
>
> I’m looking at ways to test the effect of using SIMD extensions like
> AVX2 on the performance of GSL, and I was wondering if there are
> established benchmarks for the various components that I could use.
>
> Any suggestions?
>
> Thanks in advance,
> Ludo’.
>

Re: Question on regularized regression.

[Patrick Alken]

Hello,

  Unfortunately, as you note, support for underdetermined least squares
problems is a bit lacking in GSL. GSL's SVD is currently implemented
only for M >= N. There may be other ways to solve your problem however.
If I understand, you are solving,

min_x || b - Ax ||^2 + \lambda^2 || L x ||^2

You mention ridge regression, which I think usually implies L = I, but
for now we can keep a general L matrix. The above LS problem is
equivalent to:

min_x || [ b ] - [    A     ] x ||^2
      || [ 0 ]   [ lambda*L ]   ||


For your problem, does the augmented matrix [ A ; lambda*L ] have M >=
N? If so, you can use the SVD or QR routines on this matrix to get the
solution. Though it will be a little tedius if you want to compute an
L-curve to find the optimal lambda parameter.

If your augmented matrix still has M < N, then you could use the normal
equations approach,

x = (A^T A + lambda^2 L^T L)^{-1} A^T b

This is not ideal, but works well in practice if the regularized leads
to a reasonably well-conditioned normal equations matrix.

If neither of these solutions work for you, then you might try using the
LAPACK library, which has an SVD routine which works for M < N.

Best,
Patrick

On 11/20/20 9:35 AM, Raymond Salvador wrote:
> Dear colleagues of gsl,
>
> I've been trying to use the regularized regression functions provided by
> the gsl library with an X matrix with nrows < ncols (i.e. with more
> variables than observations, a usual scenario in ridge regression) but
> prior to model fitting a svd is requested (gsl_multifit_linear_svd()) but
> this function does no allow matrices with nrows < ncols. Should I use the
> gsl_multifit_linear_tsvd() instead as a previous step to run
> gsl_multifit_linear_solve()? Thanks a lot for your attention,
>
> Raymond Salvador
>
>

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

[Patrick Alken]

Moving this to help-gsl since it is not a bug.

There are few patches people have submitted which I need to include. I
am also working on updates to the B-spline module which is not quite
ready. Its hard to find the time to get everything ready.

Is there a particular feature you need which was done after the last
release? You can always use the git version.

Patrick

On 8/18/20 11:37 AM, Tomasz Kłoczko wrote:
> URL:
>   
>
>  Summary: Any plans to make new release?
>  Project: GNU Scientific Library
> Submitted by: kloczek
> Submitted on: Tue 18 Aug 2020 05:37:03 PM UTC
> Category: None
> Severity: 3 - Normal
> Operating System: 
>   Status: None
>  Assigned to: None
>  Open/Closed: Open
>  Release: 
>  Discussion Lock: Any
>
> ___
>
> Details:
>
> There is no to many commits since last release so it should be easy to flush
> current list of commits tear ago :)
>
>
>
>
> ___
>
> Reply to this item at:
>
>   
>
> ___
>   Message sent via Savannah
>   https://savannah.gnu.org/
>
>

Re: The maximum number of subdivisions was exceeded

[Patrick Alken]

I cannot find that error anywhere in romberg.c. It looks more like a QAG
error.

What type of function are you integrating? Are there singularities? Have
you tried increasing n?

On 6/19/20 9:35 AM, Patrick Dupre wrote:
> Hello,
>
> Using Romberg, I get: The maximum number of subdivisions was exceeded
> with n= 25.
>
> How can I avoid it?
>
> Using cquad, I also got troubles!
>
> Thanks
>
> ===
>  Patrick DUPRÉ | | email: pdu...@gmx.com
>  Laboratoire interdisciplinaire Carnot de Bourgogne
>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>  Tel: +33 (0)380395988
> ===
>
>

Re: [Help-gsl] Complex QR decomposition

[Patrick Alken]

Christian and all,

  Again sorry for taking so long to working on your patch. Your Level-2
based complex QR has been added to GSL, and I have also just added the
recursive Level-3 complex QR algorithm today. I recommend using this
interface instead of the Level-2 version for speed.

So now GSL has fast (Level-3) versions of LU, Cholesky and QR, and each
has real and complex interfaces.

Enjoy,
Patrick

On 6/3/20 4:52 AM, Christian Krueger wrote:
> Hi Patrick,
>
> I don't want to insist on getting my patches introduced. But as I have,
> meanwhile, learned a bit more about git and such, I thought I would make
> my patch more easily accessible. I've merged my changes onto the current
> git commit (as of today -- bb7f532c) and attach two clean patches (one
> for the Householder fix and one for the complex QR) that could be applied.
>
> Cheers,
> Christian
>
>
> On 11/27/17 10:34 PM, Patrick Alken wrote:
>> Hi Christian,
>>
>>   This should be very useful. I will take a look at your code as soon
>> as I can, though it may be a while as I am very busy with other things
>> at the moment. I will add this to the bug tracker so it isn't lost.
>>
>> Thanks,
>> Patrick
>>
>> On 11/27/2017 08:34 AM, Christian Krueger wrote:
>>> Hi,
>>>
>>> I was looking for the QR decomposition for complex-valued matrices. As
>>> GSL does not provide this yet, I have copied the code for the
>>> real-valued QR and modified it to handle complex-valued matrices.
>>>
>>> As I've seen in the archive that this feature has been requested a few
>>> times already, I am happy to share it with everyone who wants to use it.
>>> My github repository can be found at:
>>> https://github.com/mangroveck/gsl-qr-complex-devel.git
>>>
>>> The file linalg/qrc.c is essentially a copy of the file linalg/qr.c and
>>> provides a bunch of gsl_linalg_complex_QR_* functions. I've had to add
>>> "_complex" to the function calls and sometimes take the conjugate of tau
>>> when calling the householder-functions but the whole process was fairly
>>> straight-forward and the decomposition works.
>>>
>>> * My repository also includes my patch for the
>>> gsl_linalg_complex_householder_hv() function to handle the N=1 case.
>>> (I've just reported the bug on the bug-list).
>>> * I've added various tests to the file linalg/test.c to test the new
>>> gsl_linalg_complex_QR_* functions (in the same style as the tests for
>>> the real-valued QR decomp).
>>> * I've also adjusted the header file linalg/gsl_linalg.h and the file
>>> Makefile.am
>>>
>>> The only exception is the function gsl_linalg_complex_QR_update() which
>>> I haven't adjusted yet as it requires Givens rotations and they don't
>>> take complex values yet (within GSL). All other functionality from qr.c
>>> is now available for complex values, too.
>>>
>>> If it helps, I'm happy for this code to be included in GSL.
>>>
>>> Christian
>>>
>>>
>>
>>
>>
>>

Re: [Help-gsl] Complex QR decomposition

[Patrick Alken]

Dear Christian,

  Thank you for following up, and I apologize for taking so long to get
back to you. I have looked at your patch and added the complex QR
routines, with some modifications. The main one is I have deprecated the
complex_householder_hm function because it is based on Level-1 BLAS
operations which are quite slow. I have made a new function
complex_householder_left() which uses Level-2 BLAS, at the cost of
requiring extra workspace. The end result is that the tau parameter to
QR_decomp must be length N, instead of MIN(M,N) now.

I have made a similar change to the real QR_decomp function, in a way
which is backward compatible with the current interface.

I updated the docs and put everything on the git now.

Ultimately, I would like to implement the Level-3 BLAS QR algorithm from
Elmroth and Gustavson (which is already implemented for the real case),
since the current Level-2 code is quite slow for larger matrices which
are routinely encountered in research. I hope to get around to this
soon, or if you are interested, please go ahead and try to implement it.

Thanks again for the delay,
Patrick

On 6/3/20 4:52 AM, Christian Krueger wrote:
> Hi Patrick,
>
> I don't want to insist on getting my patches introduced. But as I have,
> meanwhile, learned a bit more about git and such, I thought I would make
> my patch more easily accessible. I've merged my changes onto the current
> git commit (as of today -- bb7f532c) and attach two clean patches (one
> for the Householder fix and one for the complex QR) that could be applied.
>
> Cheers,
> Christian
>
>
> On 11/27/17 10:34 PM, Patrick Alken wrote:
>> Hi Christian,
>>
>>   This should be very useful. I will take a look at your code as soon
>> as I can, though it may be a while as I am very busy with other things
>> at the moment. I will add this to the bug tracker so it isn't lost.
>>
>> Thanks,
>> Patrick
>>
>> On 11/27/2017 08:34 AM, Christian Krueger wrote:
>>> Hi,
>>>
>>> I was looking for the QR decomposition for complex-valued matrices. As
>>> GSL does not provide this yet, I have copied the code for the
>>> real-valued QR and modified it to handle complex-valued matrices.
>>>
>>> As I've seen in the archive that this feature has been requested a few
>>> times already, I am happy to share it with everyone who wants to use it.
>>> My github repository can be found at:
>>> https://github.com/mangroveck/gsl-qr-complex-devel.git
>>>
>>> The file linalg/qrc.c is essentially a copy of the file linalg/qr.c and
>>> provides a bunch of gsl_linalg_complex_QR_* functions. I've had to add
>>> "_complex" to the function calls and sometimes take the conjugate of tau
>>> when calling the householder-functions but the whole process was fairly
>>> straight-forward and the decomposition works.
>>>
>>> * My repository also includes my patch for the
>>> gsl_linalg_complex_householder_hv() function to handle the N=1 case.
>>> (I've just reported the bug on the bug-list).
>>> * I've added various tests to the file linalg/test.c to test the new
>>> gsl_linalg_complex_QR_* functions (in the same style as the tests for
>>> the real-valued QR decomp).
>>> * I've also adjusted the header file linalg/gsl_linalg.h and the file
>>> Makefile.am
>>>
>>> The only exception is the function gsl_linalg_complex_QR_update() which
>>> I haven't adjusted yet as it requires Givens rotations and they don't
>>> take complex values yet (within GSL). All other functionality from qr.c
>>> is now available for complex values, too.
>>>
>>> If it helps, I'm happy for this code to be included in GSL.
>>>
>>> Christian
>>>
>>>
>>
>>
>>
>>

Re: Using gsl_matrix_float and gsl_vector_float in for linalg and multiroots functions

[Patrick Alken]

Hello,

  This is a great idea and should be done at some point, but
unfortunately GSL does not currently support this. The gsl_linalg
routines are mainly for the double types (gsl_vector and gsl_matrix) and
some routines are also supported for complex double (such as Cholesky
and LU).

If you need to do linear algebra on floats, you'll likely need to use
LAPACK. Or better yet, implement them in GSL and send me a patch :)

Patrick

On 3/25/20 9:16 PM, Gregory Zajac wrote:
> Is there a recommended way to get linalg functions
> like gsl_linalg_QR_decomp or multiroots functions
> like gsl_multiroot_function_fdf for gsl_matrix_float and gsl_vector_float
> types? I have not yet found equivalent functions defined in the manual nor
> in the source code. Being able to use the float instead of double types
> could help me save a lot of RAM.

Re: Integration

[Patrick Alken]

You could also directly try the QAGI algorithm, though I have never used it:

https://www.gnu.org/software/gsl/doc/html/integration.html#qagi-adaptive-integration-on-infinite-intervals

On 3/5/20 6:49 AM, Patrick Alken wrote:
> Hello, did you try transforming the integral to have finite limits (i.e.
> https://www.youtube.com/watch?v=fkxAlCfZ67E). Once you have it in this
> form, I would suggest trying the CQUAD algorithm:
>
> https://www.gnu.org/software/gsl/doc/html/integration.html#cquad-doubly-adaptive-integration
>
> Patrick
>
> On 3/5/20 2:02 AM, Patrick Dupre wrote:
>> Hello,
>>
>>
>> Can I collect your suggestions:
>>
>> I need to make the following integration:
>>
>> int_a^b g(x) f(x) dx
>>
>> where a can be 0 of -infinity, and b +infinity
>> g(x) is a Gaussian function
>> f(x) = sum (1/((x-x0)^2 + g)) / (1 + S* sum (1 / ((x-x0)^2 + g)))
>>
>> Typically, f(x) is a fraction whose numerator is a sum of Lorentzians
>> and the denominator is 1 + the same sum of Lorentzians weighted by a factor.
>>
>> Thank for your suggestions
>>
>> ===
>>  Patrick DUPRÉ | | email: pdu...@gmx.com
>>  Laboratoire interdisciplinaire Carnot de Bourgogne
>>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>>  Tel: +33 (0)380395988
>> ===
>>
>>
>

Re: Integration

[Patrick Alken]

Hello, did you try transforming the integral to have finite limits (i.e.
https://www.youtube.com/watch?v=fkxAlCfZ67E). Once you have it in this
form, I would suggest trying the CQUAD algorithm:

https://www.gnu.org/software/gsl/doc/html/integration.html#cquad-doubly-adaptive-integration

Patrick

On 3/5/20 2:02 AM, Patrick Dupre wrote:
> Hello,
>
>
> Can I collect your suggestions:
>
> I need to make the following integration:
>
> int_a^b g(x) f(x) dx
>
> where a can be 0 of -infinity, and b +infinity
> g(x) is a Gaussian function
> f(x) = sum (1/((x-x0)^2 + g)) / (1 + S* sum (1 / ((x-x0)^2 + g)))
>
> Typically, f(x) is a fraction whose numerator is a sum of Lorentzians
> and the denominator is 1 + the same sum of Lorentzians weighted by a factor.
>
> Thank for your suggestions
>
> ===
>  Patrick DUPRÉ | | email: pdu...@gmx.com
>  Laboratoire interdisciplinaire Carnot de Bourgogne
>  9 Avenue Alain Savary, BP 47870, 21078 DIJON Cedex FRANCE
>  Tel: +33 (0)380395988
> ===
>
>

Re: [Help-gsl] "gsl_spmatrix_float_alloc_nzmax" undefined reference

[Patrick Alken]

Your MWE program compiles and runs ok for me. Did you check if you might
be linking to an old version of GSL with -lgsl? You may have the new gsl
library installed somewhere different.

On 10/14/19 4:48 AM, 2592680259 wrote:
> I was using "g++ test_gsl.cpp -I/usr/local/include -lgsl -lgslcblas -o 
> test_gsl.o", but got error message:
>
>
> undefined reference to `gsl_spmatrix_float_alloc_nzmax'
> undefined reference to `gsl_spmatrix_float_d2sp'
>
>
> those two functions are declared in "gsl/gsl_spmatrix_float.h".
>
>
> The MWE for "test_gsl.cpp" is listed below
>
>
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
> #include 
>
>
> using namespace std;
>
>
> const int N = 32 * 40 * 2;
> const int N_nz = 25530;
>
>
> int main(int argc, char **argv)
> {
> gsl_matrix *test = gsl_matrix_alloc(N, N);
> gsl_matrix_float *test_float = gsl_matrix_float_alloc(N, N);
>
>
> gsl_spmatrix *test_sp = gsl_spmatrix_alloc_nzmax(N, N, N_nz, 
> GSL_SPMATRIX_COO);
> gsl_spmatrix_d2sp(test_sp, test);
>
>
> gsl_spmatrix_float *test_float_sp = gsl_spmatrix_float_alloc_nzmax(N, N, 
> N_nz, GSL_SPMATRIX_COO);
> gsl_spmatrix_float_d2sp(test_float_sp, test_float);
>
>
> gsl_matrix_free(test);
> gsl_matrix_float_free(test_float);
> gsl_spmatrix_free(test_sp);
> gsl_spmatrix_float_free(test_float_sp);
>
>
>
>
> return 0;
> }

Re: [Help-gsl] Understanding 'ftol' in nonlinear-least-squares fitter

[Patrick Alken]

Hello. Only xtol and gtol are currently used by
gsl_multifit_nlinear_test(). You can see the code here:

http://git.savannah.gnu.org/cgit/gsl.git/tree/multifit_nlinear/convergence.c

I made a stub for ftol as well, because the book by Dennis and Schnabel
talks about a test for testing the residual vector (see comment (3) in
the header of that function). However, in practice I found that test
wasn't as useful as the xtol and gtol tests, so it is currently disabled.

Also, xtol does not measure distance to the initial guess. It measures
the distance from the previous iteration - so if the step vector is
small relative to the current position, the iteration will stop.

Patrick

On 10/4/19 12:51 PM, Spencer Fleming wrote:
> Hello!
>
> I'm an undergraduate research student working at Boise State University.
> Keeping things brief, my job is to work on a project that aims to use the
> GNU Scientific Library to fit a large amount of data for scientific
> analysis.
>
> We are using the gsl_multifit_nlinear module to fit each individual unit of
> data. There are three coefficients it uses that are of interest to me
> currently: xtol, ftol and gtol.
>
> To my knowledge, each of these are used in a separate calculation for
> error, each of which answer a different question about how close our
> current guess is to the actual data.
>
> To make things a little simpler, let's say the theoretical, ideal function
> for the data we are fitting is named f(x), and the function we are
> generating to approximate it is named g(x).
>
> - xtol roughly answers "How close are we to the original guess? Did we
> wander too far?"
> - gtol roughly answers "For every data point x, how close is f(x) to g(x)?"
> - ftol I'm not as certain about. I was hoping I could be provided with some
> insight on this coefficient.
>
> Also, I apologize if my vocabulary is a bit confusing here; as a Computer
> Science major, a lot of this is new to me. However, don't shy away from
> giving me a complex answer if that's what it takes to understand ftol, or
> anything about the fitter. I can decipher it with a textbook and some grit!
>
> Thank you for your time!
>
> -Spencer Fleming

Re: [Help-gsl] gsl-config –version bring the old version

[Patrick Alken]

When you run 'make install' in gsl-2.6, it should install gsl-config in 
your prefix/bin directory, so if you ran:

./configure --prefix=~/usr

then gsl-config should be in ~/usr/bin/gsl-config

On 8/29/19 12:58 PM, mohamed jabir wrote:
> Nothing
>
>
> De : Alan Mead [mailto:am...@alanmead.org]
> Envoyé : 29 août 2019 14:57
> À : mohamed jabir ; help-gsl@gnu.org
> Objet : Re: [Help-gsl] gsl-config –version bring the old version
>
> What is the output of:
>
> $ find ~ -name 'gsl-config'
>
> ?
>
> On 8/29/2019 1:49 PM, mohamed jabir wrote:
>
> Hi,
>
>
>
> I have installed gsl (make install …) following the step described in this 
> link : https://coral.ise.lehigh.edu/jild13/2016/07/11/hello/
>
>
>
> I don’t have root privileges. When I run gsl-config –version I am getting the 
> old version (1.15). I need the latest version to used it insidethe R software.
>
>
>
> How can I make the gsl-config use the new installed version in my home 
> directory?
>
>
>
> Best,
>
> Mohamed
>
>
>
> --
>
>
>
> Alan D. Mead, Ph.D.
>
> President, Talent Algorithms Inc.
>
>
>
> science + technology = better workers
>
>
>
> http://secure-web.cisco.com/1JmLX4-otOmQUm4NDAwGK473vffIa92sqftQH6BNmBLwHrbzIvhNmhGCIU6PgfsuzW8dnHBriUSUlXkDN0bAQ8NOBG-lEQKxQMXgBprvdpAWEdxwGCBxuqp2hBoWPPdhg15HFqf96RM_QFqwgvLv0i_tJ-W-LMc5qEKhGgNtfqyMygv_7kB7dQbFBCn-qLZ22231gU_DNAHk6E8CIFk15hqqmBvZ_vm2ehSYe2ftO7wQ7TcJPNxXyHUnlTBZD_z9pGHDLPqFtFPXEJYaXUtwJ1dxqpYhwUAbXlO8hOEfs030-5nluZE_1XqKUB8sN8QfAAQ58FOREc8fjERDop1kGV6I_HDTFrAh8EkBBb435fdY02WQGWVFlL0uWYtUvbBEQYEQl6i5B1-s4INnl9hamjktKgmOflnLJmL1sYIT-1NSfZ72fKVFk3Xnl9_pZCkdFbsgZbPIg-MXh_49DihXWoQ/http%3A%2F%2Fwww.alanmead.org
>
>
>
> "You're an interesting species. An interesting mix.
>
> You're capable of such beautiful dreams, and such
>
> horrible nightmares. You feel so lost, so cut off,
>
> so alone, only you're not. See, in all our
>
> searching, the only thing we've found that makes
>
> the emptiness bearable, is each other."
>
>
>
> -- Carl Sagan, Contact
>

Re: [Help-gsl] latest version of GSL available for Red Hat

[Patrick Alken]

GSL can also be easily installed in your home directory,

./configure --prefix=~/usr
make && make install

Then just make sure you link with the correct library in ~/usr/lib

On 8/29/19 10:12 AM, Alan Mead wrote:
> You might have better luck asking on a RHEL or CentOS list. In the repos
> I have enabled, I only see 1.15 is available
>
> $ yum search gsl
> Loaded plugins: fastestmirror
> Loading mirror speeds from cached hostfile
>   * base: bay.uchicago.edu
>   * extras: bay.uchicago.edu
>   * updates: bay.uchicago.edu
> =
> N/S matched: gsl
> =
> gsl-devel.i686 : Libraries and the header files for GSL development
> gsl-devel.x86_64 : Libraries and the header files for GSL development
> gsl.i686 : The GNU Scientific Library for numerical analysis
> gsl.x86_64 : The GNU Scientific Library for numerical analysis
>
>    Name and summary matches only, use "search all" for everything.
> $ yum info gsl
> Loaded plugins: fastestmirror
> Loading mirror speeds from cached hostfile
>   * base: bay.uchicago.edu
>   * extras: bay.uchicago.edu
>   * updates: bay.uchicago.edu
> Available Packages
> Name    : gsl
> Arch    : i686
> Version : 1.15
> Release : 13.el7
> Size    : 903 k
> Repo    : base/7/x86_64
> Summary : The GNU Scientific Library for numerical analysis
> URL : http://www.gnu.org/software/gsl/
> License : GPLv3 and GFDL and BSD
> Description : The GNU Scientific Library (GSL) is a collection of
> routines for
>      : numerical analysis, written in C.
>
>
> But Fedora has 2.6:
>
> http://rpmfind.net/linux/rpm2html/search.php?query=gsl
>
> You could try that on RHEL, but I suspect it will tell you that there
> are unmet dependencies. Try using the earliest Fedora gsl rpm available.
>
> I love the stability of RHEL, but if you want the latest versions of any
> software, RHEL is a poor choice. Fedora would be a much better choice.
>
> -Alan
>
>
> On 8/29/2019 9:28 AM, mohamed jabir wrote:
>> Hi,
>>
>> What is the latest version of GSL available for "Red Hat Enterprise Linux 
>> Server release 7.7 (Maipo)"?
>>
>> I mean the version that can be installed using the Yum command.
>>
>> How can uninstall  the 1.15 version and install the 2.5 version of gsl?
>>
>> Best,
>> Mohamed

[Help-gsl] GNU Scientific Library 2.6 released

[Patrick Alken]

Version 2.6 of the GNU Scientific Library (GSL) is now available. GSL 
provides a large collection of routines for numerical computing in C.

This release introduces some new features and fixes several bugs. The 
full NEWS file entry is appended below.

The file details for this release are:

ftp://ftp.gnu.org/gnu/gsl/gsl-2.6.tar.gz
ftp://ftp.gnu.org/gnu/gsl/gsl-2.6.tar.gz.sig

The GSL project homepage is http://www.gnu.org/software/gsl/

GSL is free software distributed under the GNU General Public License.

Thanks to everyone who reported bugs and contributed improvements.

Patrick Alken

---

* What is new in gsl-2.6:

** add BLAS calls for the following functions:
  - gsl_vector_memcpy
  - gsl_vector_scale
  - gsl_matrix_memcpy
  - gsl_matrix_transpose_memcpy
  - gsl_matrix_tricpy
  - gsl_matrix_transpose_tricpy

** deprecated functions gsl_linalg_complex_householder_hm and
    gsl_linalg_complex_householder_mh

** add unit tests for gsl_linalg_symmtd and gsl_linalg_hermtd

** multilarge TSQR algorithm has been converted to use the new Level 3 
QR decomposition

** nonlinear least squares Cholesky solver now uses the new Level 3 BLAS
    method; the old modified Cholesky solver is still available under
    gsl_multifit_nlinear_solver_mcholesky and 
gsl_multilarge_nlinear_solver_mcholesky

** implemented Level 3 BLAS versions of several linear algebra routines:
  - Triangular matrix inversion
  - Cholesky decomposition and inversion (real and complex)
  - LU decomposition and inversion (real and complex)
  - QR decomposition (courtesy of Julien Langou)
  - Generalized symmetric/hermitian eigensystem reduction to 
standard form

** removed deprecated function gsl_linalg_hessenberg()

** renamed gsl_interp2d_eval_e_extrap() to gsl_interp2d_eval_extrap_e()
    to match documentation (reported by D. Lebrun-Grandie)

** renamed some of the gsl_sf_hermite functions to be more consistent
    with rest of the library, and deprecated old function names

** updated gsl_sf_hermite_func() to use a newer algorithm
    due to B. Bunck which is more stable for large x; also added
    gsl_sf_hermite_func_fast() which uses the faster Cauchy integral
    algorithm in the same paper by Bunck

** add gsl_vector_axpby()

** add un-pivoted LDLT decomposition and its banded
    variant (gsl_linalg_ldlt_* and gsl_linalg_ldlt_band_*)

** add binary search tree module (gsl_bst); based on GNU libavl

** remove -u flag to gsl-histogram

** updated spmatrix module
    - added routines and data structures for all types (float,uint,char,...)
    - added gsl_spmatrix_scale_columns() and gsl_spmatrix_scale_rows()
    - added gsl_spmatrix_add_to_dense()
    - more efficient reallocation of COO/triplet matrices (no longer 
rebuilds binary tree)
    - enhanced test suite
    - added gsl_spmatrix_min_index()

** add routines for banded Cholesky decomposition 
(gsl_linalg_cholesky_band_*)

** documented gsl_linalg_LQ routines and added gsl_linalg_LQ_lssolve()

Re: [Help-gsl] test release for GSL 2.6

[Patrick Alken]

Hi Sergey,

  I don't know what you mean by "polynomial entropy" - can you explain more?

Patrick

On 8/17/19 11:26 AM, Sergey Shcherbina wrote:
> Hi, Patrick,
>
> Congrats on  2.5.90!
>
> I have successfully used your GSL for coordinates calculation - the exactness 
> is equal about 
>
> 0.8-1.2  meters on distance about 2 kilometres. How about to add a polynomial 
> entropy in your 
>
> new GSL scientific system?  
>
> Regards,
> Sergey.

Re: [Help-gsl] test release for GSL 2.6

[Patrick Alken]

Thanks Rhys,

  I had not thought about other directory builds. I will add that to the todo 
list for next time :)

Patrick

On 8/17/19 10:51 AM, Rhys Ulerich wrote:
Hey Patrick,

Congrats on 2.5.90!

./configure && make && make check

Works for me on Ubuntu 18.04 x86_64, additionally trying a VPATH build, using 
-j for parallelism, and running 'make -j install' for a non-standard prefix.

cd doc ; make html

'make html' works in a regular build.  I've not looked at the Sphinx docs 
before.  They're pretty.

Does not work from a VPATH build (i.e. unpack GSL sources, run autogen.sh, 
create a wholly separate directory, then run GSL's configure from inside that 
separate directory).  In particular, with messages...

make[1]: Leaving directory '/home/rhys/Build/gsl-gcc/doc/examples'
Error: Config directory doesn't contain a conf.py file.

VPATH building of documentation is an edge case.  Likely not worth pursuing a 
fix at all and certainly not worth holding up 2.6.  I just tried it because I 
tend to use VPATH for autoconfiscated projects.

- Rhys

[Help-gsl] test release for GSL 2.6

[Patrick Alken]

Dear all,

  I would like to release the next version of GSL (2.6) since there have
been significant updates and the usual bug fixes. I have uploaded a test
release to:

ftp://alpha.gnu.org/gnu/gsl/gsl-2.5.90.tar.gz
ftp://alpha.gnu.org/gnu/gsl/gsl-2.5.90.tar.gz.sig

All reports are welcome - anyone who can test on various platforms would
be appreciated (Linux, BSD, Mac OS, Windows).

Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need
to install the python-based sphinx software first:

1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than
system-wide)

3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: [Help-gsl] Interrupting and then resuming minimization

[Patrick Alken]

Can you just save the current position vector (state->x) and then start
from there the next time?

On 8/6/19 5:06 PM, Pedro Gardete wrote:
> Hi,
>
> I'm running a minimization that may need to be interrupted some times, due
> to computational constraints.
>
> Is there a way I can save the "gsl_multimin_minimizer" variable, so I can
> resume optimization? The tricky part is that this struct has a void pointer
> with unknown size... I assume the remaining fields of the struct simply
> depend on the chosen minimization method.
>
> Thank you!
> Pedro

Re: [Help-gsl] finite element assembly and calling UMFPACK for linear system solution

[Patrick Alken]

For sorting the indices, take a look at csr_sort_indices from SciPy:

https://github.com/scipy/scipy/blob/v1.3.0/scipy/sparse/sparsetools/csr.h (line 
358)

I'll probably end up doing something similar in GSL when I find the time. 
You'll need to convert everything to CSC (sorting row indices instead of 
column).

For the parallel assembly, you'll need to clone the git and compile/install it.

Then, you'll do something like this:

gsl_spmatrix * A = gsl_spmatrix_alloc(M, N);

// non-parallel assembly for GSL to learn sparse pattern
loop i, j
  gsl_spmatrix_set(A, i, j, Aij);

// at this point, GSL has built a binary tree structure of the sparse pattern 
of A; now set a flag to tell gsl_spmatrix_set not to allow modifications of 
that tree
// structure
A->spflags |= GSL_SPMATRIX_FLG_FIXED;

// now you can safely call gsl_spmatrix_set from multiple threads (taking care 
not to set the same element from different threads of course
#pragma omp parallel for
for (i = 0; i < M; ++i) {
  for (j = 0; j < N; ++j) {
gsl_spmatrix_set(A, i, j, Aij); // different threads are changing different 
elements of A at the same time
  }
}

Note that this feature is still experimental and may change slightly before the 
next release. I would be quite interested to hear your experiences with it, and 
whether you have any suggestions for improvement. I have not looked in detail 
at how other sparse matrix libraries do parallel assembly, so there may be 
better ways of doing it, but so far I have used this in my work with good 
success.

Patrick

On 6/19/19 2:12 PM, Brijesh Upadhaya wrote:
This would be awesome.

I have a nonlinear magnetic problem, so for now I have reduced the number of 
unknowns (Coarse FE mesh). As far as I know the non-zero pattern of the 
nonlinear part of stiffness matrix should remain fixed. The matrix gets 
assembled in every Newton iteration. I think the new multi-threaded feature 
will come handy. I can then increase the number of unknowns as well (use dense 
FE mesh :) ) .

I'll ask you when I get hold of few things for now. I need to make this UMFPACK 
work so I am trying to sort column indices and corresponding values.

BR,
Brijesh

On Wed, Jun 19, 2019 at 10:45 PM Patrick Alken 
mailto:al...@colorado.edu>> wrote:
A new feature on the git is the ability to parallelize the assembly of the 
sparse matrix. This feature is not yet documented, but I can give you more 
details if interested.

Essentially you do a first assembly (non-parallized) so GSL can learn the 
non-zero pattern of the matrix. Then, if you want to change the matrix entries 
later (without changing the non-zero pattern), you can essentially "lock" the 
binary tree structure so that no new nodes will be inserted. Then you can 
change existing matrix entries from multiple threads without worrying.

This is useful if you need to build many sparse matrices quickly which all 
share the same non-zero pattern. Let me know if you want more details.

Patrick

On 6/19/19 1:40 PM, Brijesh Upadhaya wrote:
Thank Patrick for your prompt response.

I saw lots of new update in git.gsl. And yes now I assemble the entire matrix 
and add boundary conditions before compressing it. For now I'll just clone and 
compile the newer version.

BR,
Brijesh

On Wed, Jun 19, 2019 at 7:41 PM Patrick Alken 
mailto:al...@colorado.edu>> wrote:
Hello,

On 6/19/19 9:15 AM, Brijesh Upadhaya wrote:
> Hi everyone,
>
> I am working on a finite element problem and wanted to ask if anyone of you
> have tried using UMFPACK for the linear algebra. I have following findings
> when using sparse matrix module from gsl-2.5.
>
> 1. In CCS the row indices are not sorted so UMFPACK gives error (see
> attached example)
Its on my todo list to add a "sorted CSC" and "sorted CSR" spmatrix
type. Hopefully it will get done for the next release.
> 2. Adding boundary conditions were also difficult once compressed, so I
> couldn't utilize gsl_spmatrix_add() to add mass matrix with stiffness
> matrix. Instead I had to do it with triplet and calling the get/set
> routines. Could only add boundary condition at triplet level and compressed
> at the end and passed the pointers to UMFPACK routines.
Yes the triplet (COO) format is designed for easy set/get operations,
using a binary tree for fast lookup and insertion. Once you compress the
matrix, the binary tree disappears, and you can only modify existing
matrix entries - you cannot add new ones. Can't you assemble the entire
matrix, including boundary conditions, before compressing it?
> 3. Another thing (not related to gsl) was that UMFPACK don't use (size_t
> *), instead it has (int *) and (long *) for column pointers and row
> indices. I tried to cast (size_t *) to (long *) to utilize the pointers of
> gsl_spmatrix object.
I have decided to change the size_t pointers to int in gsl_spmatrix, fo

Re: [Help-gsl] finite element assembly and calling UMFPACK for linear system solution

[Patrick Alken]

A new feature on the git is the ability to parallelize the assembly of the 
sparse matrix. This feature is not yet documented, but I can give you more 
details if interested.

Essentially you do a first assembly (non-parallized) so GSL can learn the 
non-zero pattern of the matrix. Then, if you want to change the matrix entries 
later (without changing the non-zero pattern), you can essentially "lock" the 
binary tree structure so that no new nodes will be inserted. Then you can 
change existing matrix entries from multiple threads without worrying.

This is useful if you need to build many sparse matrices quickly which all 
share the same non-zero pattern. Let me know if you want more details.

Patrick

On 6/19/19 1:40 PM, Brijesh Upadhaya wrote:
Thank Patrick for your prompt response.

I saw lots of new update in git.gsl. And yes now I assemble the entire matrix 
and add boundary conditions before compressing it. For now I'll just clone and 
compile the newer version.

BR,
Brijesh

On Wed, Jun 19, 2019 at 7:41 PM Patrick Alken 
mailto:al...@colorado.edu>> wrote:
Hello,

On 6/19/19 9:15 AM, Brijesh Upadhaya wrote:
> Hi everyone,
>
> I am working on a finite element problem and wanted to ask if anyone of you
> have tried using UMFPACK for the linear algebra. I have following findings
> when using sparse matrix module from gsl-2.5.
>
> 1. In CCS the row indices are not sorted so UMFPACK gives error (see
> attached example)
Its on my todo list to add a "sorted CSC" and "sorted CSR" spmatrix
type. Hopefully it will get done for the next release.
> 2. Adding boundary conditions were also difficult once compressed, so I
> couldn't utilize gsl_spmatrix_add() to add mass matrix with stiffness
> matrix. Instead I had to do it with triplet and calling the get/set
> routines. Could only add boundary condition at triplet level and compressed
> at the end and passed the pointers to UMFPACK routines.
Yes the triplet (COO) format is designed for easy set/get operations,
using a binary tree for fast lookup and insertion. Once you compress the
matrix, the binary tree disappears, and you can only modify existing
matrix entries - you cannot add new ones. Can't you assemble the entire
matrix, including boundary conditions, before compressing it?
> 3. Another thing (not related to gsl) was that UMFPACK don't use (size_t
> *), instead it has (int *) and (long *) for column pointers and row
> indices. I tried to cast (size_t *) to (long *) to utilize the pointers of
> gsl_spmatrix object.
I have decided to change the size_t pointers to int in gsl_spmatrix, for
the reason you cite. Many external libraries use integer arrays instead
of size_t. This has already been done on the git repository (you can
just clone the gsl git).
>
> I am bit curious to know if anyone have a better experience. I wanted to
> use direct solver instead of iterative solver.
I would also like to include a direct solver in GSL itself, though I
don't know when I'll have time to implement it.
>
> King Regards,
> Brijesh

Re: [Help-gsl] finite element assembly and calling UMFPACK for linear system solution

[Patrick Alken]

Hello,

On 6/19/19 9:15 AM, Brijesh Upadhaya wrote:
> Hi everyone,
>
> I am working on a finite element problem and wanted to ask if anyone of you
> have tried using UMFPACK for the linear algebra. I have following findings
> when using sparse matrix module from gsl-2.5.
>
> 1. In CCS the row indices are not sorted so UMFPACK gives error (see
> attached example)
Its on my todo list to add a "sorted CSC" and "sorted CSR" spmatrix 
type. Hopefully it will get done for the next release.
> 2. Adding boundary conditions were also difficult once compressed, so I
> couldn't utilize gsl_spmatrix_add() to add mass matrix with stiffness
> matrix. Instead I had to do it with triplet and calling the get/set
> routines. Could only add boundary condition at triplet level and compressed
> at the end and passed the pointers to UMFPACK routines.
Yes the triplet (COO) format is designed for easy set/get operations, 
using a binary tree for fast lookup and insertion. Once you compress the 
matrix, the binary tree disappears, and you can only modify existing 
matrix entries - you cannot add new ones. Can't you assemble the entire 
matrix, including boundary conditions, before compressing it?
> 3. Another thing (not related to gsl) was that UMFPACK don't use (size_t
> *), instead it has (int *) and (long *) for column pointers and row
> indices. I tried to cast (size_t *) to (long *) to utilize the pointers of
> gsl_spmatrix object.
I have decided to change the size_t pointers to int in gsl_spmatrix, for 
the reason you cite. Many external libraries use integer arrays instead 
of size_t. This has already been done on the git repository (you can 
just clone the gsl git).
>
> I am bit curious to know if anyone have a better experience. I wanted to
> use direct solver instead of iterative solver.
I would also like to include a direct solver in GSL itself, though I 
don't know when I'll have time to implement it.
>
> King Regards,
> Brijesh

Re: [Help-gsl] permutation matrix with lu decomposition problem

[Patrick Alken]

GSL defines the permutation p acting on a vector as
(https://www.gnu.org/software/gsl/doc/html/permutation.html#permutations):

vi' = P v = v_{pi}

So for p = (2,0,1) acting on [ v1 v2 v3]^T, we would obtain [ v3 v1 v2 ]

The matrix representation of this is:

[ 0 0 1 ]
[ 1 0 0 ]
[ 0 1 0 ]

so it seems things are working correctly. Of course the two matrices you
listed are just transposes (inverses) of each other, so it really just
depends on the convention of the permutation definition. But I think in
this case things are correct and consistent.

Patrick

On 5/31/19 11:39 PM, Elias Posoukidis wrote:
> Hi everyone,
>
> i have the following problem with the permutation matrix after a lu
> -decomposition. The following code
>
> //--
>
>    double a_data[] = { 25, 5, 1, 64, 8, 1, 144, 12, 1 };
>
>     gsl_matrix_view m = gsl_matrix_view_array(a_data, 3, 3);
>
>     int s;
>
>     gsl_permutation *p = gsl_permutation_alloc(3);
>
>     gsl_linalg_LU_decomp(&m.matrix, p, &s);
>
>     for (size_t i = 0; i < 3; i++) {
>     for (size_t j = 0; j < 3; j++) {
>     printf("%3.2f,", gsl_matrix_get(&m.matrix, i, j));
>     }
>     printf("\n");
>     }
>
>     gsl_permutation_fprintf(stdout, p, " %u");
>     printf("\n");
>     gsl_permutation_free(p);
>
> //--
>
> gives for the lu-matrix
>
> 144.00,12.00,1.00,
> 0.17,2.92,0.83,
> 0.44,0.91,-0.20,
>
> which is the correct result, but for the permutation vector the
> following  {2 0 1}, which means that
>
> the permutation matrix is
>
> 0,1,0
>
> 0,0,1
>
> 1,0,0
>
> but the correct permutation matrix should be
>
> 0,0,1
>
> 1,0,0
>
> 0,1,0
>
> What am i doing wrong?
>
> kind regards,
>
> eliasp
>
>
>

Re: [Help-gsl] GNU Scientific Libary

[Patrick Alken]

Hello, you can contact me or post to the gsl-discuss mailing list.

Patrick

On 4/7/19 9:08 PM, Bernard Widynski wrote:
> I would like to suggest that a new RNG be included in the GNU Scientific
> Library.
>
> Whom should I contact regarding this?
>
> Sincererly,
>
> Bernard Widynski

Re: [Help-gsl] fread error help

[Patrick Alken]

I don't see how anyone can help you if you don't provide a minimal 
working example code

On 3/27/19 10:07 AM, Damian Kulec wrote:
> Hi,
>
> I've downloaded the GNU library because i need to do some matrix arithmetic
> for a program I am writing. First I need to read in a preexisting binary
> file which I then need to process.
>
> The issue occurs when I call gsl_matrix_fread. Debugging the issue simply
> points me to  ..\block\fwrite_source.c:64, where I get the warning: "ERROR:
> fread failed. Default GSL error handler invoked". I've gone through the GNU
> documentation and any forum I was able to find online referencing this
> specific error. I followed the documentation and the example provided, but
> to no avail.
>
> The answers I found on the forums were more useless. One answer said to
> "close the pointer to the file after opening it", but that just gave me
> another error: Expression: fh >= 0 && (unsigned)fh < (unsigned)_nhandle.
> Other solutions directed me to make sure I opened the file specifying that
> it's "read-only" and "binary", which I already have done. I can't figure
> out why this function is not working for me.
>
> Please let me know how to solve this.
>
> Thank you,
> Damian

Re: [Help-gsl] gsl_multifit_nlinear_fdf with combined function and Jacobian?

[Patrick Alken]

Olaf,

  Ironically, I removed the combined function to make things easier for 
users. I can confirm that the multifit_nlinear module will always call 
the Jacobian function after calling the function. So if you want you can 
store intermediate results in the function call, and then reuse them in 
the Jacobian call.


Note that sometimes the function is called without calling the Jacobian.

Patrick

On 08/15/2018 07:27 AM, Olaf Wucknitz wrote:

Hi,

I wonder if there is a reason that gsl_multifit_nlinear_fdf does not 
include a pointer to a function that computes both the function and 
its Jacobian at the same time, because this is often more efficient 
than doing it separately.


gsl_multiroot_function_fdf does have this option.

Or can we assume that, e.g., whenever the function for the Jacobian is 
called, the function itself is called with the same parameters before 
so that one can reuse intermediate steps internally?


Cheers,
Olaf

Re: [Help-gsl] non-linear least square fitting

[Patrick Alken]

GSL does not currently support constrained optimization directly. You
could, however, introduce a large penalty if your parameters are outside
the desired range. In your residual function, check if x is outside
[x_min,x_max] and set the output f_i to a large value. This will
discourage the solver from exploring that part of parameter space.

Patrick

On 08/18/2018 05:14 AM, Brijesh Upadhaya wrote:
> Hi,
>
> I am new to GSL and I want to use non-linear least-square fitting to fit a
> data set with 2000 data points. The function has 14 parameters in total. I
> tried using gsl_multifit_nlinear_driver () to solve the problem.
>
> It works nice but I want the parameters to be in certain range: x_min <=
> xpar <= x_max. How to fix the lower and upper bounds (Or I missed something
> in the gsl-2.5 documentation!). Some of the parameters in my model don't
> have any physical meaning if not bounded.
>
> I have another question as well!. I want to do a multi-objective
> optimization using GSL-2.5 but I couldn't locate any routine in the
> documentation. I have 40 data set having 2000 points each and I need to
> obtain 14 parameters (Global minimum).
>
> with kind regards,
> Brijesh

Re: [Help-gsl] delete function generates an error in movstat.h

[Patrick Alken]

Daniel,

  Sorry for the slow reponse. The problem is that "delete" is a C++ 
keyword which g++ doesn't like in that header file. I always use gcc to 
compile so I didn't see this problem. Please edit your gsl_movstat.h 
header file and change "delete" to "delete_sample". Then try compiling 
your code.


I will make this update in the git repository also.

Patrick

On 08/15/2018 02:25 PM, Montezano, Daniel wrote:

Hello there,

Just checking if anyone has come across a similar problem or if anyone has been 
able to reproduce the error with this code.

On my side, so far no luck. Have reinstalled everything and 
recompiled/reinstalled GSL and still the issue remains.

Thanks in advance for any help or suggestions.

Daniel

Sent: Thursday, July 26, 2018 10:25 AM
To: help-gsl@gnu.org
Subject: Re: [Help-gsl] delete function generates an error in movstat.h

Hi Patrick,

Thanks for the reply.
The source code I am trying to compile follows below.

The file vector.dat contains the following list of numbers:
12, -20, 12, 24, 23, 8, -19, 9, -4, 2, 3, 4, 1, 7, 6, 7, 0, -1, -1, 2, 44, 1, 
-9, 7, 5, 5, 6, 6, 1, 1

Thanks!

Daniel



--BEGIN SOURCE CODE--
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

using namespace std;
int main(int argc, char** argv) {
 // window length
 const size_t K = 13;
 // time series length
 const size_t N = 30;


 // input time series
 gsl_vector *x = gsl_vector_alloc (N);
 {
 FILE * f = fopen ("vector.dat","r");
 gsl_vector_fscanf (f,x);
 fclose (f);
 }


 // output mean time series
 gsl_vector *xmean = gsl_vector_alloc(N);
 // workspace for the rolling computation
 gsl_movstat_workspace * wkspc = gsl_movstat_alloc (K);


 // compute the rolling mean for a gsl_vector
 gsl_movstat_mean (GSL_MOVSTAT_END_PADZERO, x, xmean, wkspc);

 // important to free these things to prevent crashes
 gsl_vector_free(x);
 gsl_vector_free(xmean);
 gsl_movstat_free (wkspc);

 return 0;
}
-- END SOURCE CODE-



Can you provide your source code?

On 07/13/2018 07:31 PM, Daniel wrote:

Hi,

I am trying to compile an example from the docs. Compute a rolling mean using
the new package Moving Window Statistics.

My code is producing the following error (Windows 10, MinGW.org GCC-6.3.0-1,
GSL 2.5):



g++ gsl_rolling_stats.cpp -lgsl -o gsl_rolling_stats.exe
  In file included from gsl_rolling_stats.cpp:11:0:
  c:\documents\progrs\mingw32\include\gsl\gsl_movstat.h:58:9: error: expected
  unqualified-id before 'delete'
  int (*delete) (void * vstate);
   ^~
  c:\documents\progrs\mingw32\include\gsl\gsl_movstat.h:58:9: error: expected
  ')' before 'delete'

I compiled GSL compiled from source using MinGW and msys. I have tested my
GSL installation with other GSL code and everything is working fine.

I appreciate any help.

Thank you!

Daniel

Re: [Help-gsl] delete function generates an error in movstat.h

[Patrick Alken]

Can you provide your source code?

On 07/13/2018 07:31 PM, Montezano, Daniel wrote:
> Hi,
>
> I am trying to compile an example from the docs. Compute a rolling mean using 
> the new package Moving Window Statistics.
>
> My code is producing the following error (Windows 10, MinGW.org GCC-6.3.0-1, 
> GSL 2.5):
>
>
>> g++ gsl_rolling_stats.cpp -lgsl -o gsl_rolling_stats.exe
>>  In file included from gsl_rolling_stats.cpp:11:0:
>>  c:\documents\progrs\mingw32\include\gsl\gsl_movstat.h:58:9: error: expected
>>  unqualified-id before 'delete'
>>  int (*delete) (void * vstate);
>>   ^~
>>  c:\documents\progrs\mingw32\include\gsl\gsl_movstat.h:58:9: error: expected
>>  ')' before 'delete'
> I compiled GSL compiled from source using MinGW and msys. I have tested my 
> GSL installation with other GSL code and everything is working fine.
>
> I appreciate any help.
>
> Thank you!
>
> Daniel

[Help-gsl] GNU Scientific Library 2.5 released

[Patrick Alken]

Version 2.5 of the GNU Scientific Library (GSL) is now available. GSL provides 
a large collection of routines for numerical computing in C.

This release introduces some new features and fixes several bugs. The full NEWS 
file entry is appended below.

The file details for this release are:

ftp://ftp.gnu.org/gnu/gsl/gsl-2.5.tar.gz 
<ftp://ftp.gnu.org/gnu/gsl/gsl-2.4.tar.gz>


ftp://ftp.gnu.org/gnu/gsl/gsl-2.5.tar.gz.sig 
<ftp://ftp.gnu.org/gnu/gsl/gsl-2.4.tar.gz.sig>


The GSL project homepage ishttp://www.gnu.org/software/gsl/

GSL is free software distributed under the GNU General Public License.

Thanks to everyone who reported bugs and contributed improvements.

Patrick Alken



* What is new in gsl-2.5:

** doc bug fix in binomial distribution figure (Damien Desfontaines)

** added Wishart distribution (Timothée Flutre)

** added new module for digital filtering (gsl_filter); current filters include:
 Gaussian filter
 median filter
 recursive median filter
 impulse detection filter

** added new module for moving window statistics (gsl_movstat)

** added statistics functions:
 gsl_stats_median()
 gsl_stats_select()
 gsl_stats_mad()
 gsl_stats_mad0()
 gsl_stats_Sn_from_sorted_data()
 gsl_stats_Qn_from_sorted_data()
 gsl_stats_gastwirth_from_sorted_data()
 gsl_stats_trmean_from_sorted_data()

** added Romberg integration (gsl_integration_romberg)

** bug fix in deprecated functions gsl_multifit_wlinear_svd and
   gsl_multifit_wlinear_usvd (reported by Vlad Koli)

** documention corrected to state that gsl_sf_legendre functions do
   not include Condon-Shortley phase by default

** bug fix in exponential fitting example when using larger number
   of points (reported by Anna Russo)

** changed internal workspace inside gsl_spmatrix to a union to
   avoid casting (suggested by Manuel Schmitz)

** bug fixes in ode-initval2 for very rare solver crashing cases:
   #52230 in msadams (reported by Michael Kaufman) and
   #52336 in msbdf (reported by Andrew Benson). As a fix,
   the maximum scaling of controlled step length was decreased
   from 5.0 to 4.9.

** add histogram2d figure to manual (was missing in 2.4)

** bug fix in gsl_spmatrix_add for duplicate input arguments
   (reported by Alfredo Correa)

** add support for negative arguments nu in gsl_sf_bessel_Jnu and
   gsl_sf_bessel_Ynu (Konrad Griessinger)

** better texinfo documentation for gsl_sf_hyperg functions

** fix vector and matrix fread/fwrite testing on windows systems
   when tmpfile() fails

** fix for rstat/test.c on PPC64 (reported by Adam Majer)

Re: [Help-gsl] test release for GSL 2.5

[Patrick Alken]

Hi,

  Thanks for testing. I can't reproduce those warnings here, but in any
case they don't look serious.

Patrick

On 06/04/2018 11:56 PM, Kenneth Geisshirt wrote:
> Hey,
>
> I am on MacOS 10.13 with  clang 902.0.39.2. Everything builds and all tests
> pass. I see a few warnings like the following two:
>
> ./test_source.c:584:53: note: use function 'abs' instead
> if (fabs(r - y) > 2 * GSL_FLT_EPSILON * fabs(y))
>
> ./test_rmedian.c:39:39: warning: implicit conversion from enumeration type
> 'const gsl_filter_end_t' to different enumeration type 'gsl_movstat_end_t'
> [-Wenum-conversion]
>   size_t wsize = gsl_movstat_fill(endtype, x, i, H, H, window);
>
> /Kenneth
>
> On Mon, Jun 4, 2018 at 11:18 PM Jeremy Theler  wrote:
>
>> Every test passed in Debian Sid
>>
>> gtheler@tom:~/tmp/gsl-2.4.90$ uname -a
>> Linux tom 4.16.0-2-amd64 #1 SMP Debian 4.16.12-1 (2018-05-27) x86
>> <20%2018%2005%2027>_64
>> GNU/Linux
>> gtheler@tom:~/tmp/gsl-2.4.90$ gcc --version
>> gcc (Debian 7.3.0-21) 7.3.0
>> Copyright (C) 2017 Free
>> Software Foundation, Inc.
>> This is free software; see the source for copying conditions.  There is NO
>> warranty
>> ; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
>>
>>
>> docs also compile smoothly.
>>
>> --
>> jeremy theler
>> http://secure-web.cisco.com/1FHtwb1UrpvjxNzDJaWkJkK8v-1caW5VWcUQYWjHOIUVCDHnfWdXBFTXDOxzaBn9bSYWMq7WrGoVgUU45u9aJ15ve0zYdPtcurDFB89kM7XbW6vRkk1WehxEWM65f13kMRO48A-bJbqt8gmOYMw07wgccceTCPhW4ueA_CY3-0FGBApMPTrAD4Susd6dG9N845dS9umZ0qxPIymlFSjesbQH_w0faxViTILm347bXYJuBnuDxlfxh4VIvB4JhZL56lBcnnIvQ6AIoChY0tcUpRmlngh_fNCHzaq84hVy4vmseGIWpajhVhAocAVdz0Oc9sBX9Zd3OhH3-n6HZEO4IKQ18pTmhQhIzaPACVOKPKhb7DSe1BCoDFXSbMthY79hiPWvqg2KUFwJNrn1zOmw324MNxj5EbSmI1sbWmtKPabjt1tsEZnbMy4sB859yN7EYgQWJgg0Z0AHx0uUnvw/http%3A%2F%2Fwww.seamplex.com
>>
>>
>>
>> On Mon, 2018-06-04 <20%2018%2006%2004> at 09:09 -0600, Patrick Alken
>> wrote:
>>> Dear all,
>>>
>>>I would like to release the next version of GSL (2.5) since there
>>> have been some new features added and the usual bug fixes. I have
>>> uploaded a test release to:
>>>
>>> ftp://alpha.gnu.org/gnu/gsl/gsl-2.4.90.tar.gz
>>> ftp://alpha.gnu.org/gnu/gsl/gsl-2.4.90.tar.gz.sig
>>>
>>> All reports are welcome - anyone who can test on various platforms would
>>> be appreciated (Linux, BSD, Mac OS, Windows).
>>>
>>> Please try testing the build:
>>>
>>> ./configure && make && make check
>>>
>>> If you wish you can test building the documentation, but you will need
>>> to install the python-based sphinx software first:
>>>
>>> 1. pip install -U --user Sphinx
>>>
>>> 2. pip install -U --user sphinx_rtd_theme
>>>
>>> (the --user flag will install in your home directory rather than
>>> system-wide)
>>>
>>> 3. cd doc ; make html
>>>
>>> Please report any successes/failures.
>>>
>>> Thanks,
>>> Patrick
>>>
>>>
>>>
>> --
> Kenneth Geisshirt, M.Sc., Ph.D.
> Majbøl Allé 18, DK-2770 Kastrup, +45 60 62 71 82
>

[Help-gsl] test release for GSL 2.5

[Patrick Alken]

Dear all,

  I would like to release the next version of GSL (2.5) since there 
have been some new features added and the usual bug fixes. I have 
uploaded a test release to:


ftp://alpha.gnu.org/gnu/gsl/gsl-2.4.90.tar.gz
ftp://alpha.gnu.org/gnu/gsl/gsl-2.4.90.tar.gz.sig

All reports are welcome - anyone who can test on various platforms would 
be appreciated (Linux, BSD, Mac OS, Windows).


Please try testing the build:

./configure && make && make check

If you wish you can test building the documentation, but you will need 
to install the python-based sphinx software first:


1. pip install -U --user Sphinx

2. pip install -U --user sphinx_rtd_theme

(the --user flag will install in your home directory rather than 
system-wide)


3. cd doc ; make html

Please report any successes/failures.

Thanks,
Patrick

Re: [Help-gsl] Fwd: Modifying adaptive integration function

[Patrick Alken]

Just one last note, I assume you found this already, but the
gsl_coerce_double() function is located in sys/coerce.c. You could try
calling this function directly.

On 05/11/2018 01:35 AM, Francesco Florian wrote:
> On Thursday, May 10, 2018 8:17:47 PM CEST lostbits wrote:
>> On 5/10/2018 8:08 AM, Patrick Alken wrote:
>>> On 05/10/2018 04:18 AM, Francesco Florian wrote:
>>>> Hello!
>>>> Since I have received no answer, I wonder whether this is the right 
>>>> mailing list for this question. If it is not, can you please point me 
>>>> to the right one?
>>>> Thank you
>>> Hello,
>>>
>>>   What exactly are you trying to do that QAG cannot? Perhaps there is 
>>> an alternate way without rewriting the function. My guess is you could 
>>> get away with just removing the GSL_COERCE_DBL wrappers from the code, 
>>> but it may not match the original QUADPACK results, according to the 
>>> comments in qag.c.
>>>
>>> Patrick
>>>
>> My guess is that gcc -DHAVE_EXTENDED_PRECISION_REGISTERS ... would 
>> define the pre-process variable, if that's what you want. If you do 
>> nothing then the variable is not defined. More to the point:
>>
>> #if HAVE_EXTENDED_PRECISION_REGISTERS
>> #define GSL_COERCE_DBL(x) (gsl_coerce_double(x))  // used with gcc 
>> --DHAVE_EXTENDED_PRECISION_REGISTERS
>> #else
>> #define GSL_COERCE_DBL(x) (x)  // used without defining the variable
>> #endif
> I was looking for a way to check if the compiler actually supports it. The 
> method above seems to rely on me knowing that, and then passing 
> -DHAVE_EXTENDED_PRECISION_REGISTERS to gcc if appropriate.
> Patrick Alken's answer seems to suggest it is ok not to match QUADPACK 
> results; if it is correct it would work for me.
>

Re: [Help-gsl] Fwd: Modifying adaptive integration function

[Patrick Alken]

On 05/10/2018 04:18 AM, Francesco Florian wrote:

Hello!
Since I have received no answer, I wonder whether this is the right mailing 
list for this question. If it is not, can you please point me to the right one?
Thank you


Hello,

  What exactly are you trying to do that QAG cannot? Perhaps there is 
an alternate way without rewriting the function. My guess is you could 
get away with just removing the GSL_COERCE_DBL wrappers from the code, 
but it may not match the original QUADPACK results, according to the 
comments in qag.c.


Patrick

Re: [Help-gsl] Linear fitting

[Patrick Alken]

Francisco,

  You have some errors in the way you handle your 2D arrays. I have 
made a small modification to define new 1D arrays (x and y) and store 
your fit data into those arrays. The program output now says:


Performing fit with 20 points...done! (status = 0)
== BET results ==
y = 0.22 + 0.005321 x
[2.54531e-11, -1.26944e-10
 -1.26944e-10, 8.07463e-10]
Chi Square = 0.00
=

My changes are attached - you can find the lines I changed with the "PA" 
keyword


On 05/03/2018 04:13 AM, Francisco M Neto wrote:

Greetings!

I have been trying to use gsl for some linear least squares
fitting, but for some reason I'm getting bad results.

Instead of returning the results of the fit, gsl_fit_linear()
returns -nan for all parameters. I'm not sure if I'm passing something
wrong to it, but I do know that there is no problem with the data (I've
ran it through other LS fitting procedures and results are consistent).

Attached are the source and an example data file. I've compiled
and ran with:

$ gcc -I/usr/include/gsl/ -lm -lgsl -lgslcblas -o teste_bet teste_bet.c

$ ./teste_bet data.dat



#include 
#include 
#include 
#include 
#include 

/* Maximum name size + 1. */

#define MAX_NAME_SZ 256
#define MAX_LINE_SIZE 256

int file_lines (FILE *file);

int file_read (FILE *file, int n);

int main (int argc, char *argv[]) {

  int lines, i;
  FILE *input;
  char *name, line[MAX_LINE_SIZE];
  double **data, **bet;
  double *x, *y; // PA

  double a, b, cov00, cov01, cov11, chisq; // Fit variables
  int n, status;

  if (!(name = malloc(sizeof(char)*MAX_NAME_SZ))) {
printf("No memory!\n");
return 1;
  }

  if (argc == 1) {  // No extra arguments on the command line
printf("Name of the input file: ");
fgets(name, MAX_NAME_SZ, stdin);
  } else if (argc == 2) {
strcpy(name, argv[1]);
  }
  else {
printf ("Too many arguments! Exiting.\n");
exit(1);
  }
  if ((strlen(name) > 0) && (name[strlen (name) - 1] == '\n'))
name[strlen (name) - 1] = '\0';
  
  //printf ("Loading file: %s\n", name);

  if (!(input=fopen(name, "r")))
exit(1);

  lines = file_lines (input);
  //  printf ("Lines in the file: %d\n", lines);

  if(!(data = malloc(2 * sizeof(double *
exit(1);
  for (i=0; i<2; i++)
if(!(data[i] = malloc((lines) * sizeof(double
  exit(1);
  
  for (i=0; i= 0.05 && data[0][i] <= 0.3) {
  // PA: set x and y vectors
  x[n] = data[0][i];
  y[n] = 1/(data[1][i]*((1/data[0][i])-1));

  bet[0][i] = data[0][i];
  bet[1][i] = 1/(data[1][i]*((1/data[0][i])-1)); // y' = 1/(y*(x-1))
  bet[2][i] = 1/(bet[1][i]*0.05);  // weight for the fits, based on 5% estimated error
  printf ("%.8e\t%.8e\t%.8e\n", bet[0][i], bet[1][i], bet[2][i]);
  n++;
}
  }

  // Linear fit to BET-transformed data from 0.05 < x < 0.3 (y' = a + bx')
#if 1
  {
printf ("Performing fit with %d points...", n);
status = gsl_fit_linear (x, 1, y, 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
if (status) printf ("not ");
printf ("done! (status = %d)\n", status);
  }
#else
  printf ("Performing fit with %d points...", n);
  status = gsl_fit_linear (bet[0], 1, bet[1], 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
  if (status) printf ("not ");
  printf ("done! (status = %d)\n", status);
#endif

  printf ("== BET results ==\n");
  printf ("y = %lf + %lf x\n", a, b);
  printf ("[%g, %g\n %g, %g]\n", 
  cov00, cov01, cov01, cov11);
  printf ("Chi Square = %lf\n", chisq);
  printf ("=\n");
  
  fclose (input);
  free (name);
  free (data);
  free (bet);
  return 0;
  
}


int file_lines (FILE *file) {  // Counts the lines in a text file

  int i=0, c;

  while (!feof(file)) {
c = fgetc(file);
if (c == '\n')
  i++;
  }

  rewind(file);
  return i;
  
}

Re: [Help-gsl] Problems with compiling Non Linear Fit Example

[Patrick Alken]

Hello,

  I just installed a clean version of GSL v2.4 in my home directory to 
try your program. Here are the commands I used:


> cd gsl-2.4
> ./configure --prefix=/home/palken/usr
> make
> make install
> gcc -g -Wall -o test test.c -lm /home/palken/usr/lib/libgsl.a -lgslcblas
> ./test
...
summary from method 'trust-region/levenberg-marquardt'
number of iterations: 16
function evaluations: 23
Jacobian evaluations: 17
reason for stopping: small step size
initial |f(x)| = 62.202928
final   |f(x)| = 5.454180
chisq/dof = 0.804002
A  = 5.17379 +/- 0.27938
lambda = 0.11104 +/- 0.00817
b  = 1.05283 +/- 0.05365
status = success

Note that when compiling test.c, I gave the full path to the new GSL 
v2.4 library (your issue is probably that the linker is trying to link 
with an old version of GSL). Hope this helps.


Patrick

On 02/16/2018 01:45 AM, Dr. Jonny Birkhan wrote:

Dear all,

I have installed gsl-2.4 on a

Linux 4.4.0-103-generic #126-Ubuntu SMP Mon Dec 4 16:23:28 UTC 2017 
x86_64.


The libraries are installed in /usr/local/lib.

My gcc version is: gcc (Ubuntu 5.4.0-6ubuntu1~16.04.6) 5.4.0 20160609


When I try to compile the Exponential Fitting Example I get the 
following messages:


gcc -LLIBDIR=/usr/local/lib -o test test.c -lm -lgsl -lgslcblas

/tmp/ccuzgEIN.o:: In function `callback':
test.c:(.text+0x2b1): undefined reference to 
`gsl_multifit_nlinear_residual'
test.c:(.text+0x2c1): undefined reference to 
`gsl_multifit_nlinear_position'

test.c:(.text+0x2d8): undefined reference to `gsl_multifit_nlinear_rcond'
/tmp/cclIWgcr.o: In function `main':
test.c:(.text+0x386): undefined reference to `gsl_multifit_nlinear_trust'
test.c:(.text+0x399): undefined reference to 
`gsl_multifit_nlinear_default_parameters'

test.c:(.text+0x640): undefined reference to `gsl_multifit_nlinear_alloc'
test.c:(.text+0x665): undefined reference to `gsl_multifit_nlinear_winit'
test.c:(.text+0x671): undefined reference to 
`gsl_multifit_nlinear_residual'
test.c:(.text+0x6f6): undefined reference to 
`gsl_multifit_nlinear_driver'

test.c:(.text+0x708): undefined reference to `gsl_multifit_nlinear_jac'
test.c:(.text+0x729): undefined reference to `gsl_multifit_nlinear_covar'
test.c:(.text+0x755): undefined reference to 
`gsl_multifit_nlinear_trs_name'

test.c:(.text+0x764): undefined reference to `gsl_multifit_nlinear_name'
test.c:(.text+0x78f): undefined reference to `gsl_multifit_nlinear_niter'
test.c:(.text+0xa9e): undefined reference to `gsl_multifit_nlinear_free'
collect2: error: ld returned 1 exit status

I guess it's a problem with a missing library file. I couldn't find 
any bug report on that problem.


The source code is appended.

Can you help me?

Thanks in advance and best regards,

Jonny Birkhan

Re: [Help-gsl] GSL 2.4 make check failing on Rocks cluster running CentOS 6.9

[Patrick Alken]

Ah so it looks like the problem is with the internally estimated error, 
which is quite wrong for this input value. I will need to think about a 
fix for this, and I will open a bug report to track the issue.


In the meantime, don't worry about this error since it is calculating 
the correct value to within 14 decimal places.


You can run make -k to test the rest of the library (-k won't quit when 
it encounters an error). Can you let me know if you find any more failures?


Thanks,
Patrick

On 01/18/2018 04:05 PM, Eric Shell wrote:

Hi Patrick,

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, &r) [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  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, &r),
-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, &r) [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

Re: [Help-gsl] GSL 2.4 make check failing on Rocks cluster running CentOS 6.9

[Patrick Alken]

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, &r), 
-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, &r) [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

Re: [Help-gsl] gsl_integration_qawf

[Patrick Alken]

Hello,

  A return code of 0 is equivalent to GSL_SUCCESS, so yes that is good.
I'm not sure what you mean by inner/outer integration - can you provide
a minimal working source code to demonstrate the problem you're having?

  The adaptive integration routines in GSL are based on QUADPACK - you
can find further documentation on QAWF here:
http://www.netlib.org/quadpack/doc.
There is also an example calling program for QAWF on that page (in
fortran), perhaps it will help you.

  Thanks for the whitespace issue - I have fixed it on the git.

Patrick

On 01/15/2018 01:35 AM, Richard Hartmann wrote:
> Hi
>
> I'm struggling with the numeric Fourier integration, in particular
> with the error handling of the 'gsl_integration_qawf' routine.
>
> As of the nested integration and if the error handling is default, my
> test aborts since the inner integration reports a roundoff error.
> (current Git 3bbe65ad5a07362e82103b311feda55d4355a2fe, qawf.c l. 133)
>
> So far so good, however, if the error handling is off, the outer
> integration seem to run through. The overall return code is 0,
> indicating success. Is that intensional?
>
> I tried to answer the question my self, but I could not understand
> what the 'gsl_integration_qawf' routine was doing.
>
>
> Secondly and very minor, there is a 'space' missing in the error
> message at qawo.c l. 104
>
>   GSL_ERROR ("cannot reach tolerance because of roundoff error"
>  "on first attempt", GSL_EROUND);
>
> should be
>   GSL_ERROR ("cannot reach tolerance because of roundoff error "
>  "on first attempt", GSL_EROUND);
>
>
> Thanks for this great library and thanks for you help.
> Richard
>

Re: [Help-gsl] Contribution, PhD thesis topic, Computer Architecture and Scientific Computing

[Patrick Alken]

What are your interests? There are plenty of things to work on in GSL,
including linear algebra, numerical integration, optimization, RNGs,
ODEs, special functions, bug fixing and numerous other topics.

On 01/12/2018 09:35 AM, Ali kia wrote:
> Hi !
>
> I am a first year Computer Architecture PhD student and I am working to write 
> my proposal and am interested to do something that is also useful for other. 
> So I am sending this message to ask about your requirement or any idea for 
> working on it.
>
>
> I am really interested to hear from you.
>
>
> Best regards,
>
> A. K.
>

Re: [Help-gsl] Nonlinear Least-Squares Fitting

[Patrick Alken]

Hello,

  Yes this is a bug, thanks for reporting it. The problem was I had 
defined the time variable to be in [0, N], where N is the number of 
points to fit. Then, for some values of N, after the first iteration the 
solver finds a negative value for lambda, and so exp(-lambda * N) gets 
very large, resulting in nan values in the expb_f function.


  I have now replaced the time variable with:

t_i = i * TMAX / (N - 1)

so that t \in [0,TMAX] and you can increase N to anything you want. I 
left TMAX = 40. I will update the example program in the repository. In 
the meantime I have attached the corrected source code. I tested it up 
to N = 1000. Please let me know if you have any further problems.


Thanks,
Patrick

On 01/11/2018 04:50 AM, Anna Russo Russo wrote:

Hi!

I've run the "Exponential Fitting Example", from the GSL 2.4 documentation,
and I noticed that, if the number of data points is greater than 56, the
solver terminates with status "exceeded max number of iterations". It
happens even if I allow it to perform a considerably larger number of
iterations: after the first iteration the parameters are so different from
the right ones that the value of the function to be minimized is 'nan' and
they don't change in the next iterations.

What may be the cause for this behavior?

Thanks for reading,

Anna



#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#define N  100/* number of data points to fit */
#define TMAX   (40.0) /* time variable in [0,TMAX] */

struct data {
  size_t n;
  double * t;
  double * y;
};

int
expb_f (const gsl_vector * x, void *data, 
gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *t = ((struct data *)data)->t;
  double *y = ((struct data *)data)->y;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
{
  /* Model Yi = A * exp(-lambda * t_i) + b */
  double Yi = A * exp (-lambda * t[i]) + b;
  gsl_vector_set (f, i, Yi - y[i]);
}

  return GSL_SUCCESS;
}

int
expb_df (const gsl_vector * x, void *data, 
 gsl_matrix * J)
{
  size_t n = ((struct data *)data)->n;
  double *t = ((struct data *)data)->t;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < n; i++)
{
  /* Jacobian matrix J(i,j) = dfi / dxj, */
  /* where fi = (Yi - yi)/sigma[i],  */
  /*   Yi = A * exp(-lambda * t_i) + b  */
  /* and the xj are the parameters (A,lambda,b) */
  double e = exp(-lambda * t[i]);
  gsl_matrix_set (J, i, 0, e); 
  gsl_matrix_set (J, i, 1, -t[i] * A * e);
  gsl_matrix_set (J, i, 2, 1.0);
}

  return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
 const gsl_multifit_nlinear_workspace *w)
{
  gsl_vector *f = gsl_multifit_nlinear_residual(w);
  gsl_vector *x = gsl_multifit_nlinear_position(w);
  double rcond;

  /* compute reciprocal condition number of J(x) */
  gsl_multifit_nlinear_rcond(&rcond, w);

  fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
  iter,
  gsl_vector_get(x, 0),
  gsl_vector_get(x, 1),
  gsl_vector_get(x, 2),
  1.0 / rcond,
  gsl_blas_dnrm2(f));
}

int
main (void)
{
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  gsl_multifit_nlinear_workspace *w;
  gsl_multifit_nlinear_fdf fdf;
  gsl_multifit_nlinear_parameters fdf_params =
gsl_multifit_nlinear_default_parameters();
  const size_t n = N;
  const size_t p = 3;

  gsl_vector *f;
  gsl_matrix *J;
  gsl_matrix *covar = gsl_matrix_alloc (p, p);
  double t[N], y[N], weights[N];
  struct data d = { n, t, y };
  double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
  gsl_vector_view x = gsl_vector_view_array (x_init, p);
  gsl_vector_view wts = gsl_vector_view_array(weights, n);
  gsl_rng * r;
  double chisq, chisq0;
  int status, info;
  size_t i;

  const double xtol = 1e-8;
  const double gtol = 1e-8;
  const double ftol = 0.0;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* define the function to be minimized */
  fdf.f = expb_f;
  fdf.df = expb_df;   /* set to NULL for finite-difference Jacobian */
  fdf.fvv = NULL; /* not using geodesic acceleration */
  fdf.n = n;
  fdf.p = p;
  fdf.params = &d;

  /* this is the data to be fitted */
  for (i = 0; i < n; i++)
{
  double ti = i * TMAX / (n - 1.0);
  double yi = 1.0 + 5 * exp (-0.1 * ti);
  double si = 0.1 * yi;
  double dy = gsl_ran_gaussian(r, si);

  t[i] = ti;
  y[i] = yi + dy;
  weights[i] = 1.0 / (si * si);
  printf ("data: %g %g %g\n", ti, y[i], si);
};

  /* allocate workspace with default parameters */
  w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);

  /* initialize solver with starting point and weights */
  gsl_multifit_nlinear_winit (&x.vector

Re: [Help-gsl] Suggested RNG for inclusion in the GNU Scientific Library

[Patrick Alken]

Hello,

  Do you have an implementation which can interface easily with GSL? If 
so I can put it on the webpage for people to try out.


Patrick

On 01/06/2018 08:27 PM, Bernard Widynski wrote:

A new RNG based on the middle-square method has been placed in the arXiv.

arXiv:1704.00358 
This is one of the fastest RNGs and it passes all the statistical tests.

It has been published for over nine months and there have been no reported
errors.

Perhaps you might consider including it in the GSL.

Sincerely,

Bernard Widynski

Re: [Help-gsl] [Bug-gsl] FFT not generated

[Patrick Alken]

It looks like something went wrong with the makefile generation, since 
there should be a -I.. -I../gsl flag passed to the compiler to locate 
the *.h files.


Did you try the methods here:

https://www.gnu.org/software/gsl/extras/native_win_builds.html

for windows builds? There is also a cygwin package prebuilt for gsl.

Patrick

On 01/08/2018 01:53 AM, Peter Schacht wrote:

Hi,
'make check' was not done.

I've unpacked the gsl.tar.gz in dir gsl-2.4. In this dir I built dir
gslbuild.
In gslbuild I startet

.././configure  --prefix=C:\Qt\Tools\mingw492_32\i686-w64-mingw32
make
-(I didn't make check because I didn't know that at this time)
make install

Starting from die gslbuild:
$ cd fft
$ make test
gcc -g -O2../../fft/test.c   -o test
../../fft/test.c:20:20: fatal error: config.h: No such file or directory
  #include 
 ^
compilation terminated.
make: *** [test] Error 1

so ./test could not be executed.

config.h is in dir gslbuild  but not in the minGW gsl include Path

C:\Qt\Tools\mingw492_32\i686-w64-mingw32\include\gsl

nor in

C:\Qt\Tools\mingw492_32\i686-w64-mingw32\include



Grüße Peter

258/127 Moo 3, Thepkrasatri Road
Srisoonthorn, Thalang, Phuket 83110

Mobil +66(96)926 83 62
Email  peter.scha...@schacht-co.de







-Ursprüngliche Nachricht-
Von: Patrick Alken [mailto:al...@colorado.edu]
Gesendet: Sonntag, 7. Januar 2018 20:37
An: bug-...@gnu.org; peter.scha...@kabelmail.de
Betreff: Re: [Bug-gsl] FFT not generated

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   <mailto:peter.scha...@kabelmail.de> peter.scha...@schacht-co.de

  

  

  

Re: [Help-gsl] help in compiling and linking this file

[Patrick Alken]

Sorry, the 2.1 docs are texinfo based, so do:

cd doc ; make dvi

assuming you have texinfo installed

On 01/05/2018 02:38 PM, Patrick Alken wrote:
> I don't keep them stored online, you can compile them yourself in the
> doc/ directory if you have Sphinx installed (try cd doc; make html)
>
> Alternatively the old texinfo based docs are still up here from v2.3:
>
> https://www.gnu.org/software/gsl/manual/gsl-ref.html
>
> These are v2.3 but I don't think the integration routines changed
> between 2.1 and 2.3
>
> On 01/05/2018 02:33 PM, Vasu Jaganath wrote:
>> Thanks Patrick, where can i read docs for 2.1? Because I think it is
>> difficult for me to migrate to 2.4 as yet.
>>
>>
>>
>> On Fri, Jan 5, 2018 at 2:27 PM, Patrick Alken > <mailto:al...@colorado.edu>> wrote:
>>
>> That workspace was added in GSL 2.4, so make sure you are using the
>> latest version of the library
>>
>> On 01/05/2018 01:40 PM, Mohammad Akhlaghi wrote:
>> > Hi Vasu,
>> >
>> > The "unknown type name" is not a linker error, its a compiler
>> error (which comes prior to linking).
>> >
>> > The `gsl_integration_fixed_workspace' type is indeed defined in
>> `gsl_integration.h'. So it is very strange that your compiler
>> doesn't recognize it (after the inclusion of this header)!
>> >
>> > Try having a look in the used header file and see if this type
>> is indeed defined there.
>> >
>> > Cheers,
>> > Mohammad
>> >
>> > On January 5, 2018 8:04:17 PM GMT+01:00, Vasu Jaganath
>> mailto:vadie.develo...@gmail.com>> wrote:
>> >> Hi Forum,
>> >>
>> >> I don't know what libraries to link to compiler this file
>> >>
>> >> I am doing
>> >>
>> >> gcc -Wall quadFixed.c -lgsl -lgslcblas -lm
>> >>
>> >> however it complains
>> >>
>> >> error: unknown type name ‘gsl_integration_fixed_workspace’
>> >>   gsl_integration_fixed_workspace * w;
>> >>
>> >> took this example directly from
>> >> https://www.gnu.org/software/gsl/doc/html/integration.html
>> <https://www.gnu.org/software/gsl/doc/html/integration.html>
>> >>
>> >> I have attached the file, I am able to compile and run qags example
>> >> just
>> >> fine.
>> >>
>> >> Also, how to know what .so libraries to link for what header in
>> gsl?
>> >>
>> >> Thanks,
>> >> Vasu
>>
>>
>>
>>

Re: [Help-gsl] help in compiling and linking this file

[Patrick Alken]

I don't keep them stored online, you can compile them yourself in the
doc/ directory if you have Sphinx installed (try cd doc; make html)

Alternatively the old texinfo based docs are still up here from v2.3:

https://www.gnu.org/software/gsl/manual/gsl-ref.html

These are v2.3 but I don't think the integration routines changed
between 2.1 and 2.3

On 01/05/2018 02:33 PM, Vasu Jaganath wrote:
> Thanks Patrick, where can i read docs for 2.1? Because I think it is
> difficult for me to migrate to 2.4 as yet.
>
>
>
> On Fri, Jan 5, 2018 at 2:27 PM, Patrick Alken  <mailto:al...@colorado.edu>> wrote:
>
> That workspace was added in GSL 2.4, so make sure you are using the
> latest version of the library
>
> On 01/05/2018 01:40 PM, Mohammad Akhlaghi wrote:
> > Hi Vasu,
> >
> > The "unknown type name" is not a linker error, its a compiler
> error (which comes prior to linking).
> >
> > The `gsl_integration_fixed_workspace' type is indeed defined in
> `gsl_integration.h'. So it is very strange that your compiler
> doesn't recognize it (after the inclusion of this header)!
> >
> > Try having a look in the used header file and see if this type
> is indeed defined there.
> >
> > Cheers,
> > Mohammad
> >
> > On January 5, 2018 8:04:17 PM GMT+01:00, Vasu Jaganath
> mailto:vadie.develo...@gmail.com>> wrote:
> >> Hi Forum,
> >>
> >> I don't know what libraries to link to compiler this file
> >>
> >> I am doing
> >>
> >> gcc -Wall quadFixed.c -lgsl -lgslcblas -lm
> >>
> >> however it complains
> >>
> >> error: unknown type name ‘gsl_integration_fixed_workspace’
> >>   gsl_integration_fixed_workspace * w;
> >>
> >> took this example directly from
> >> https://www.gnu.org/software/gsl/doc/html/integration.html
> <https://www.gnu.org/software/gsl/doc/html/integration.html>
> >>
> >> I have attached the file, I am able to compile and run qags example
> >> just
> >> fine.
> >>
> >> Also, how to know what .so libraries to link for what header in
> gsl?
> >>
> >> Thanks,
> >> Vasu
>
>
>
>

Re: [Help-gsl] help in compiling and linking this file

[Patrick Alken]

That workspace was added in GSL 2.4, so make sure you are using the
latest version of the library

On 01/05/2018 01:40 PM, Mohammad Akhlaghi wrote:
> Hi Vasu,
>
> The "unknown type name" is not a linker error, its a compiler error (which 
> comes prior to linking).
>
> The `gsl_integration_fixed_workspace' type is indeed defined in 
> `gsl_integration.h'. So it is very strange that your compiler doesn't 
> recognize it (after the inclusion of this header)!
>
> Try having a look in the used header file and see if this type is indeed 
> defined there.
>
> Cheers,
> Mohammad
>
> On January 5, 2018 8:04:17 PM GMT+01:00, Vasu Jaganath 
>  wrote:
>> Hi Forum,
>>
>> I don't know what libraries to link to compiler this file
>>
>> I am doing
>>
>> gcc -Wall quadFixed.c -lgsl -lgslcblas -lm
>>
>> however it complains
>>
>> error: unknown type name ‘gsl_integration_fixed_workspace’
>>   gsl_integration_fixed_workspace * w;
>>
>> took this example directly from
>> https://www.gnu.org/software/gsl/doc/html/integration.html
>>
>> I have attached the file, I am able to compile and run qags example
>> just
>> fine.
>>
>> Also, how to know what .so libraries to link for what header in gsl?
>>
>> Thanks,
>> Vasu

Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?

[Patrick Alken]

You might be able to used fixed quadrature - you will probably need a
large number of nodes to capture that sharp feature. I would compare
both the CQUAD and fixed-point methods to see which works better for you.

Patrick

On 12/31/2017 11:10 PM, Vasu Jaganath wrote:
> HI Martin, Yes one of my Q is very discontinuous with respect to my
> integration variable Z.
>
> I have attached the plots for you to see
>
> On Sun, Dec 31, 2017 at 9:20 PM, Vasu Jaganath
> mailto:vadie.develo...@gmail.com>> wrote:
>
> Yes, I will show you the plots soon, Q is actually 2 variable
> function but for my calculations I am treating one of the
> variables as a parameter, which is a physically valid assumption.
> Yes I do encounter alpha and beta values less than 1.
>
>
> On Sun, Dec 31, 2017, 9:13 PM Martin Jansche  <mailto:mjans...@gmail.com>> wrote:
>
> So you want to find E[f] = \int_0^1 f(x) dbeta(x | a, b) dx.
> Can you plot
> your typical f? And what are typical values of the parameters
> a and b? Do
> you encounter a<=1 or b<=1? If so, how does f(x) behave as x
>     approaches 0
> or 1?
>
> On Mon, Jan 1, 2018 at 3:37 AM, Patrick Alken
> mailto:al...@colorado.edu>> wrote:
>
> > The question is whether your Q contains any singularities,
> or is highly
> > oscillatory? Is such cases fixed point quadrature generally
> doesn't do
> > well. If Q varies fairly smoothly over your interval, you
> should give
> > fixed point quadrature a try and report back if it works
> well enough for
> > your problem. The routines you want are documented here:
> >
> > http://www.gnu.org/software/gsl/doc/html/integration.html#
> <http://www.gnu.org/software/gsl/doc/html/integration.html#>
> > fixed-point-quadratures
> >
> > Also, if QAGS isn't working well for you, try also the CQUAD
> routines.
> > I've found CQUAD is more robust than QAGS in some cases
> >
> > On 12/31/2017 05:28 PM, Vasu Jaganath wrote:
> > > I have attached my entire betaIntegrand function. It is a
> bit complicated
> > > and very boiler-plate, It's OpenFOAM code (where scalar =
> double etc.) I
> > > hope you can get the jist from it.
> > > I can integrate the Q using monte-carlo sampling integration.
> > >
> > > Q is nothing but tabulated values of p,rho, mu etc. I
> lookup Q using the
> > > object "solver" in the snippet.
> > >
> > > I have verified evaluating  when I am not using those
>  values back
> > in
> > > the solution, It works OK.
> > >
> > > Please ask me anything if it seems unclear.
> > >
> > >
> > >
> > >
> > >
> > >
> > > On Sun, Dec 31, 2017 at 3:55 PM, Martin Jansche
> mailto:mjans...@gmail.com>>
> > wrote:
> > >
> > >> Can you give a concrete example of a typical function Q?
> > >>
> > >> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath <
> > vadie.develo...@gmail.com <mailto:vadie.develo...@gmail.com>>
> > >> wrote:
> > >>
> > >>> Hi forum,
> > >>>
> > >>> I am trying to integrate moments, basically first order
> moments ,
> > i.e.
> > >>> averages of some flow fields like temperature, density
> and mu. I am
> > >>> assuming they distributed according to beta-PDF which is
> parameterized
> > on
> > >>> variable Z, whose mean and variance i am calculating
> separately and
> > using
> > >>> it to define the shape of the beta-PDF, I have a cut off
> of not using
> > the
> > >>> beta-PDF when my mean Z value, i.e  is less than a
> threshold.
> > >>>
> > >>> I am using qags, the adaptive integration routine to
> calculate the
> > moment
> > >>> integral, however I am restricted to threshold of  

Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?

[Patrick Alken]

Why not just plot Q(x) over your integration interval?

On 12/31/2017 09:00 PM, Vasu Jaganath wrote:
> As far as I know, my Q doesn't contain any singularities but I will
> recheck my tables again just to confirm.
>
> Thanks,
> Vasu
>
> On Sun, Dec 31, 2017 at 8:37 PM, Patrick Alken  <mailto:al...@colorado.edu>> wrote:
>
> The question is whether your Q contains any singularities, or is
> highly
> oscillatory? Is such cases fixed point quadrature generally doesn't do
> well. If Q varies fairly smoothly over your interval, you should give
> fixed point quadrature a try and report back if it works well
> enough for
> your problem. The routines you want are documented here:
>
> 
> http://www.gnu.org/software/gsl/doc/html/integration.html#fixed-point-quadratures
> 
> <http://www.gnu.org/software/gsl/doc/html/integration.html#fixed-point-quadratures>
>
> Also, if QAGS isn't working well for you, try also the CQUAD routines.
> I've found CQUAD is more robust than QAGS in some cases
>
> On 12/31/2017 05:28 PM, Vasu Jaganath wrote:
> > I have attached my entire betaIntegrand function. It is a bit
> complicated
> > and very boiler-plate, It's OpenFOAM code (where scalar = double
> etc.) I
> > hope you can get the jist from it.
> > I can integrate the Q using monte-carlo sampling integration.
> >
> > Q is nothing but tabulated values of p,rho, mu etc. I lookup Q
> using the
> > object "solver" in the snippet.
> >
> > I have verified evaluating  when I am not using those 
> values back in
> > the solution, It works OK.
> >
> > Please ask me anything if it seems unclear.
> >
> >
> >
> >
> >
> >
> > On Sun, Dec 31, 2017 at 3:55 PM, Martin Jansche
> mailto:mjans...@gmail.com>> wrote:
> >
> >> Can you give a concrete example of a typical function Q?
> >>
> >> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath
> mailto:vadie.develo...@gmail.com>>
> >> wrote:
> >>
> >>> Hi forum,
> >>>
> >>> I am trying to integrate moments, basically first order
> moments , i.e.
> >>> averages of some flow fields like temperature, density and mu.
> I am
> >>> assuming they distributed according to beta-PDF which is
> parameterized on
> >>> variable Z, whose mean and variance i am calculating
> separately and using
> >>> it to define the shape of the beta-PDF, I have a cut off of
> not using the
> >>> beta-PDF when my mean Z value, i.e  is less than a threshold.
> >>>
> >>> I am using qags, the adaptive integration routine to calculate
> the moment
> >>> integral, however I am restricted to threshold of  = 1e-2.
> >>>
> >>> It complains that the integral is too slowly convergent. However
> >>> physically
> >>> my threshold should be around 5e-5 atleast.
> >>>
> >>> I can integrate these moments with threshold upto 5e-6, using
> Monte-Carlo
> >>> integration, by generating random numbers which are
> beta-distributed.
> >>>
> >>> Should I look into fixed point integration routines? What
> routines would
> >>> you suggest?
> >>>
> >>> Here is the (very simplified) code snippet where, I calculate
> alpha and
> >>> beta parameter of the PDF
> >>>
> >>>                     zvar   = min(zvar,0.*zvar_lim);
> >>>                     alpha = zmean*((zmean*(1 - zmean)/zvar) - 1);
> >>>                     beta = (1 - zmean)*alpha/zmean;
> >>>
> >>>                     // inside the fucntion to be integrated
> >>>                     // lots of boilerplate for Q(x)
> >>>                     f = Q(x) * gsl_ran_beta_pdf(x, alpha, beta);
> >>>
> >>>                    // my integration call
> >>>
> >>>                    helper::gsl_integration_qags (&F, 0, 1, 0,
> 1e-2, 1000,
> >>>                                                   w, &result,
> &error);
> >>>
> >>> And also, I had to give relative error pretty large, 1e-2.
> However some of
> >>> Qs are pretty big order of 1e6.
> >>>
> >>> Thanks,
> >>> Vasu
> >>>
> >>
>
>
>

Re: [Help-gsl] fixed point or adaptive integration for calculating moments using beta PDF?

[Patrick Alken]

The question is whether your Q contains any singularities, or is highly
oscillatory? Is such cases fixed point quadrature generally doesn't do
well. If Q varies fairly smoothly over your interval, you should give
fixed point quadrature a try and report back if it works well enough for
your problem. The routines you want are documented here:

http://www.gnu.org/software/gsl/doc/html/integration.html#fixed-point-quadratures

Also, if QAGS isn't working well for you, try also the CQUAD routines.
I've found CQUAD is more robust than QAGS in some cases

On 12/31/2017 05:28 PM, Vasu Jaganath wrote:
> I have attached my entire betaIntegrand function. It is a bit complicated
> and very boiler-plate, It's OpenFOAM code (where scalar = double etc.) I
> hope you can get the jist from it.
> I can integrate the Q using monte-carlo sampling integration.
>
> Q is nothing but tabulated values of p,rho, mu etc. I lookup Q using the
> object "solver" in the snippet.
>
> I have verified evaluating  when I am not using those  values back in
> the solution, It works OK.
>
> Please ask me anything if it seems unclear.
>
>
>
>
>
>
> On Sun, Dec 31, 2017 at 3:55 PM, Martin Jansche  wrote:
>
>> Can you give a concrete example of a typical function Q?
>>
>> On Sat, Dec 30, 2017 at 3:42 AM, Vasu Jaganath 
>> wrote:
>>
>>> Hi forum,
>>>
>>> I am trying to integrate moments, basically first order moments , i.e.
>>> averages of some flow fields like temperature, density and mu. I am
>>> assuming they distributed according to beta-PDF which is parameterized on
>>> variable Z, whose mean and variance i am calculating separately and using
>>> it to define the shape of the beta-PDF, I have a cut off of not using the
>>> beta-PDF when my mean Z value, i.e  is less than a threshold.
>>>
>>> I am using qags, the adaptive integration routine to calculate the moment
>>> integral, however I am restricted to threshold of  = 1e-2.
>>>
>>> It complains that the integral is too slowly convergent. However
>>> physically
>>> my threshold should be around 5e-5 atleast.
>>>
>>> I can integrate these moments with threshold upto 5e-6, using Monte-Carlo
>>> integration, by generating random numbers which are beta-distributed.
>>>
>>> Should I look into fixed point integration routines? What routines would
>>> you suggest?
>>>
>>> Here is the (very simplified) code snippet where, I calculate alpha and
>>> beta parameter of the PDF
>>>
>>> zvar   = min(zvar,0.*zvar_lim);
>>> alpha = zmean*((zmean*(1 - zmean)/zvar) - 1);
>>> beta = (1 - zmean)*alpha/zmean;
>>>
>>> // inside the fucntion to be integrated
>>> // lots of boilerplate for Q(x)
>>> f = Q(x) * gsl_ran_beta_pdf(x, alpha, beta);
>>>
>>>// my integration call
>>>
>>>helper::gsl_integration_qags (&F, 0, 1, 0, 1e-2, 1000,
>>>   w, &result, &error);
>>>
>>> And also, I had to give relative error pretty large, 1e-2. However some of
>>> Qs are pretty big order of 1e6.
>>>
>>> Thanks,
>>> Vasu
>>>
>>

Re: [Help-gsl] Complex QR decomposition

[Patrick Alken]

Hi Christian,

  This should be very useful. I will take a look at your code as soon 
as I can, though it may be a while as I am very busy with other things 
at the moment. I will add this to the bug tracker so it isn't lost.


Thanks,
Patrick

On 11/27/2017 08:34 AM, Christian Krueger wrote:

Hi,

I was looking for the QR decomposition for complex-valued matrices. As
GSL does not provide this yet, I have copied the code for the
real-valued QR and modified it to handle complex-valued matrices.

As I've seen in the archive that this feature has been requested a few
times already, I am happy to share it with everyone who wants to use it.
My github repository can be found at:
https://github.com/mangroveck/gsl-qr-complex-devel.git

The file linalg/qrc.c is essentially a copy of the file linalg/qr.c and
provides a bunch of gsl_linalg_complex_QR_* functions. I've had to add
"_complex" to the function calls and sometimes take the conjugate of tau
when calling the householder-functions but the whole process was fairly
straight-forward and the decomposition works.

* My repository also includes my patch for the
gsl_linalg_complex_householder_hv() function to handle the N=1 case.
(I've just reported the bug on the bug-list).
* I've added various tests to the file linalg/test.c to test the new
gsl_linalg_complex_QR_* functions (in the same style as the tests for
the real-valued QR decomp).
* I've also adjusted the header file linalg/gsl_linalg.h and the file
Makefile.am

The only exception is the function gsl_linalg_complex_QR_update() which
I haven't adjusted yet as it requires Givens rotations and they don't
take complex values yet (within GSL). All other functionality from qr.c
is now available for complex values, too.

If it helps, I'm happy for this code to be included in GSL.

Christian

Re: [Help-gsl] legendre poly deriv

[Patrick Alken]

I think the older routine from 1.16 has the same issue. For Legendre
derivatives, I believe the current algorithm divides by sin(theta), so
at the poles there is a singularity, therefore the routine fails for x=1
or x=-1 as you found.

The following paper discusses the problem and proposes an algorithm to
fix it:

http://www.sciencedirect.com/science/article/pii/S146418951010

It is not currently implemented in GSL but if you're interested in
working on it I would be glad to accept a patch for this

Patrick

On 11/10/2017 07:09 AM, Li Dong wrote:
> Hi,
>
> First of all, thanks for this great library! This is my first post in this
> group. Sorry if there is any rule I fail to follow.
>
> I was using gsl 1.16 and recently upgraded to 2.3. I found
> gsl_sf_legendre_Plm_deriv_array is replaced by gsl_sf_legendre_deriv_array.
> Maybe in most cases, this replacement is just fine. However, I find that in
> certain cases there is no way to get a result from current
> gsl_sf_legendre_deriv_array.
>
> In 1.16, the following code works.
> double L[5], DL[5]; // here 5 is just an arbitrary large number to
> initialize the array
> int lmax = 3, m = 2;
> double x=-1.0;
> gsl_sf_legendre_Plm_deriv_array (lmax, m, x, L, DL); // If m=1, this line
> doesn't work since x=-1.
>
> In 2.3, it needs to be written as:
> double L[5], DL[5]; \\ 5 is again arbitrary
> int lmax = 3;
> double x = -1.0;
> gsl_sf_legendre_deriv_array (GSL_SF_LEGENDRE_NONE, lmax, x, L, DL); // this
> line won't work since we calculate all 0<=m<=lmax, it breaks when
> calculates m=1 and x=-1
>
> I wonder if there is a way around to implement this?
>
> Thanks,
> Li

Re: [Help-gsl] Nonlinea fitting with exit info = 27

[Patrick Alken]

Hello, there are some tips in the Troubleshooting section of the
Nonlinear least squares manual:

http://www.gnu.org/software/gsl/doc/html/nls.html#troubleshooting

In particular, your problem appears to be an error in your Jacobian
matrix. I've attached a modified version of your program where I
switched to the finite difference Jacobian calculation. I also lowered
the solution tolerance to get a more accurate solution. You'll see that
the numerical solution now agrees with your model to about 1e-7,
although the computed parameters are not exactly what you started with.
This could likely be improved further when you find the error in your
Jacobian function (which I have not done).

Here are the differences with your original program:

> diff main.c main.c.orig
> 188,189c188,189
> <   const double xtol = 1.0e-6;
> <   const double gtol = 1.0e-6;
> ---
> >   const double xtol = 1.0e-2;
> >   const double gtol = 1.0e-2;
> 276d275
> < #if 0
> 279,282d277
> < #else
> <   fdf.df = NULL;
> <   fdf.fvv = NULL;
> < #endif
> 295c290
> <   fdf_params.trs = gsl_multifit_nlinear_trs_dogleg;
> ---
> >   fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;

Note I also switched to the dogleg algorithm - in order to use lmaccel
you need an accurate Jacobian routine and fvv routine.

On 10/16/2017 07:56 PM, Jiawei Zhao Zhao wrote:
> Dear GSL mailing list.
>
>
> I am new to GSL and trying to fit a Non linear function. the function is 
> given below.
>
>
> y = A*( ( b / (x-r0) ) ^ m- B* ( b / (x-r0) ) ^ n)
>
>
> I modified an example of non linear fitting (code attached). The function 
> successfully complied  and first deviation gives correct value.
>
>
> However, when I test with generating data from:
>
>
> A=1, b= 2.5, m=12, n=6, r0=0 and B=1
>
>
> and try to fit it start from
>
>
> A=1, b= 2.4, m=12, n=6, r0=0 and B=1
>
>
> The iteration seems stopped at the first steps with info = 27.
>
>
> Can anyone help me on this?
>
>
> Regards.
>
>
> Jiawei Zhao
>
>

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

struct data
{
  double *t;
  double *y;
  size_t n;
};


double
gaussian(const double a, const double b, const double c, const double t)
{
  const double z = (t - b) / c;
  return (a * exp(-0.5 * z * z));
}

  //const double a = 1; //energy
  //const double b = 2.5;  //sigema,distance
  //const double m = 12;   /* m */
  //const double n = 6;/* n */
  //const double r0 = 0;   /* cutoff */
  //const double B = 1;	// Extra attractive


int
func_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double m = gsl_vector_get(x, 2);
  double n = gsl_vector_get(x, 3);
  double r0 = gsl_vector_get(x, 4);
  double B = gsl_vector_get(x, 5);
  size_t i;

  for (i = 0; i < d->n; ++i)
{
  double ti = d->t[i];
  double yi = d->y[i];
  double y = a*(pow(b/(ti-r0),m)-B*pow(b/(ti-r0),n));//gaussian(a, b, c, ti);

  gsl_vector_set(f, i, yi - y);
}

  return GSL_SUCCESS;
}

int
func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double m = gsl_vector_get(x, 2);
  double n = gsl_vector_get(x, 3);
  double r0 = gsl_vector_get(x, 4);
  double B = gsl_vector_get(x, 5);
  size_t i;

  for (i = 0; i < d->n; ++i)
{
  double ti = d->t[i];
  gsl_matrix_set(J, i, 0, pow((-b/(r0 - ti)),m) - B*pow((-b/(r0 - ti)),n));
  gsl_matrix_set(J, i, 1, -a*((m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti))  );
  gsl_matrix_set(J, i, 2, a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m));
  gsl_matrix_set(J, i, 3, -B*a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n));
  gsl_matrix_set(J, i, 4, a*((b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti)));  
  gsl_matrix_set(J, i, 5, -a*pow((-b/(r0 - ti)),n));


}

  return GSL_SUCCESS;
}

int
func_fvv (const gsl_vector * x, const gsl_vector * v,
  void *params, gsl_vector * fvv)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double m = gsl_vector_get(x, 2);
  double n = gsl_vector_get(x, 3);
  double r0 = gsl_vector_get(x, 4);
  double B = gsl_vector_get(x, 5);
  double va = gsl_vector_get(v, 0);
  double vb = gsl_vector_get(v, 1);
  double vm = gsl_vector_get(v, 2);
  double vn = gsl_vector_get(v, 3);
  double vr0 = gsl_vector_get(v, 4);
  double vB = gsl_vector_get(v, 5);
  size_t i;

  for (i = 0; i < d->n; ++i)
{
  double ti = d->t[i];


  //double Daa = 0;
  double Dab = (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti) - (m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti);
  double Dam = log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m);
  double Dan = -B*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n);
  double Dar0 = 

Re: [Help-gsl] Thin QR decomposition

[Patrick Alken]

On 10/09/2017 06:03 AM, Ines Arous wrote:
> Hi GSL team,
>
> I am trying to use the QR function that is provided by GSL (which is very 
> effective compared to other QR implementation I found online). My problem is 
> that it computes the full decomposition while I need the thin QR 
> decomposition.
> Is there a way to compute the thin QR decomposition using GSL? If so can you 
> please point me on how I can do that?
>
> Best regards,
> Thanks!
> Inès Arous

What exactly do you need to do? Are you trying to form the explicit Q1
and Q2 matrices? You almost never need to do this.

If you just want the product Q^T b, you can use gsl_linalg_QR_QTvec().
If you want to compute the residual norm || y - A*x ||, this can be done
without explicitly forming Q1 and Q2. 

If you really do need Q1, Q2 and R(NxN), you can use
gsl_linalg_QR_unpack() and just take submatrix views to get each factor.

Patrick

Re: [Help-gsl] contribute code for Wishart distribution

[Patrick Alken]

Timothée,

  Yes can you send me a patch against the git master branch?

Thanks,
Patrick

On 07/17/2017 09:45 PM, Timothée Flutre wrote:

Hello,

I wrote some code to sample random matrices from the Wishart distribution,
as well as compute the (log) probability density function, along with tests.

For that I also added some code to cholesky.c to solve AX=B (where X is a
matrix) where A is provided via its Cholesky factor (but I didn't wrote any
test for that as I don't really understand the existing tests.)

Patrick, can I send you the different modified files directly?

Best,
Timothée Flutre

Re: [Help-gsl] test release for GSL 2.4

[Patrick Alken]

I tried Mohammad's suggestion, and the tests now write two files in the
current working directory: "test.dat" and "test_static.dat". This
appears to be working correctly on my system with make check -j8. All
function calls to tmpfile() have been removed.

It would be great if a windows person can test the latest git repository
to see if the issue is fixed

Thanks,
Patrick

On 06/22/2017 10:41 AM, Brian Gladman wrote:
> I assume that by 'race condition' you mean that any tests running in
> parallel would both be using the same file (IIRC "test.dat").
>
> If so, you could possibly revert your change and use the tmpnam()
> function to set the filename to be used.
>
> I don't know where the resulting file would be placed with GCC but on
> Windows this will return a filename that is unique in the current
> working directory, which would be likely to have the right access
> permission.  However, the MS documentation simply says that
>
> "tmpnam returns a name unique in the current working directory"
>
> which suggests that this won't solve the race condition.
>
> We could possibly maintain full c89 compatibility for GSL itself but
> relax this requirement for these two tests and use mkstemp and the MSVC
> equivalent (_mktemp or _mktemp_s) that seem to avoid the problem.
>
> I cant test for these failures however as they don't occur for me as I
> am running the build with admin privileges.
>
>best regards,
>
>  Brian
>
>

Re: [Help-gsl] test release for GSL 2.4

[Patrick Alken]

Its not a race condition between the files written by matrix/vector. The
vector module itself has two test programs, which are identical except
one tests the non-inline versions of the functions and the other tests
the inline versions. So if these two tests are run in parallel, via make
check -j8, they will write the same test file name at the same time if
we use a static filename.

Patrick

On 06/22/2017 09:56 AM, Mohammad Akhlaghi wrote:
> Hi Patrick,
>
> How about using two separate files (file names) for vector and matrix
> tests?
>
> For example "test-vector.dat" and "test-matrix.dat".
>
> Cheers,
> Mohammad
>
> On June 22, 2017 8:12:06 AM GMT+01:00, Patrick Alken
>  wrote:
>
> In GSL 2.3 and earlier, the vector and matrix modules tested the
> fwrite/fread routines by creating temporary files with mkstemp() and
> fdopen(). Unfortunately neither of these routines conform to the ANSI
> C89 standard and so I converted these tests to write to a fixed filename
> "test.dat". However this caused the file race condition which was
> reported in the recent test candidates (i.e. the "make check -j8"
> error). So finally I switched to use tmpfile() which is C89 and seems to
> be the only ANSI-compatible method of using temporary files in C. But
> now it seems that some windows systems fail with this method due to
> permission restrictions.
>
> I can certainly modify the tests to check for a NULL pointer coming out
> of tmpfile, but this would then disable these tests on your windows
> systems, which is not ideal either. It appears we must either accept
> that the matrix/vector fread/fwrite routines cannot be properly tested
> on windows systems, or go back to a non-ANSI method of writing the
> temporary files. Neither of these choices seem good to me, but my
> preference would be to follow the C89 standard if possible.
>
> I thought about writing a routine inside GSL to perform the same task as
> mkstemp() does, but it seems this is not possible in C89 since mkstemp
> relies on calling open() with O_EXCL to raise an error if the file
> already exists, but open() is also not C89.
>
> I am open to suggestions if anyone knows of a solution to this problem.
>
> Patrick
>
> On 06/21/2017 03:26 PM, Max Gaspa wrote:
>
> Dear all, Sisyphus wrote
>
> Much the same situation here - Windows 7 64 bit, mingw-w64
> port of 64-bit gcc-6.3.0 (x86_64-posix-sjlj-rev1). And,
> like you (apparently), I didn't test the release
> candidates :-( 
>
> Yes. And Yes. I didn't check the release candidate because the
> short time between the availability of the RC and the
> availability of the official release. Anyway I think I found
> the reason of the issue just following the program flow with
> gdb. The issue is NOT related to malloc(0) (I' still thinking
> that is not a good programming practice but I will accept it
> :-) ) I'm discussing the issue for vector just to describe it
> in a simpler way. The reason of the segmentation fault is the
> latest modification made by Patrick 7 days ago into
> test_source.c ( described by "fix file race condition for
> 'make check -j8'" ). The lines involved are - char filename[]
> = "test.dat"; + FILE *f = tmpfile(); in few words the change
> imply using tmpfile(). Unfortunately in Windows (XP is fine! 7
> fails) that call is creating a temporary file in C:\ that is
> usually not writable for security reason. The call fails and
> the pointer f is NULL. Because there are no checking for NULL
> after the call to tmpfile() (It seems GSL developer love not
> checking for null pointer :-) :-) :-) ) and the next call
> to fprintf will use NULL as its stream you get the
> segmentation fault. Now I'm trying to revert the change (just
> using the release candidate version should be fine) to see if
> the issue is fixed, but I'm quite sure because gdb told me
> that f was NULL but it was used as a valid pointer for the
> stream For reference: Microsoft documentation of the C runtime
> library (used by MinGW) for tmpfile() *** The tmpfile()
> function creates a temporary file and returns a pointer to
> that stream. The temporary file is created in the root
> directory. To create a temporary file in a directory other
> than the root, use tmpnam or tempnam in conjunction with
> fopen.