[Bug-gsl] Possible bug in linalg/hh.c Householder solution to Ax = b (patch)

[Dave Wiltshire]

I think there is a bug in the Linear Algebra Householder solver for the
problem Ax = b.  The function gsl_linalg_HH_svx incorrectly calls an over
determined system under determined (and vice versa).  Looking in to the
code I see that the dimensions are mixed up such that the code will only
work correctly for a square matrix and not for the over determined case.

I've written a test case based on the expected output from case 4 from
here:
http://people.sc.fsu.edu/~jburkardt/f_src/qr_solve/qr_solve_prb_output.txt
that is attached.

I've also attached a patch file that will fix the problem for
gsl_linalg_HH_svx (apply with patch linalg/hh.c < gsl_linalg_hh.c.patch).

Looking at the code there is a bigger problem with the function
gsl_linalg_HH_solve which assumes that the length of the vector b and the
length of the vector x (in Ax = b) are equal, copying the data from b in to
x before passing this to the function gsl_linalg_HH_svx.  This is only true
in a square case and in the over determined case x will have less elements
than b, such that this copy can't be done.  The solution to this problem is
harder (probably allocate a working array of size M  the length of b, and
copy vector b there and then copy the first N elements out in to x ...
using standard notation of A being an M x N matrix).  Anyway, I've left
that alone.

Dave

gsl_hh_test.c
Description: Text Data

gsl_linalg_hh.c.patch
Description: Text Data