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
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_linalg.h"
#include <stdio.h>
#define M 5
#define N 3
int main(void)
{
double A[M][N] = {{1.0, 1.0, 1.0},
{1.0, 2.0, 4.0},
{1.0, 3.0, 9.0},
{1.0, 4.0, 16.0},
{1.0, 5.0, 25.0}};
double b[M] = {1.0, 2.3, 4.6, 3.1, 1.2};
gsl_vector *vb;
gsl_matrix *mA;
int rc;
vb = gsl_vector_alloc(M);
mA = gsl_matrix_alloc(M, N);
for (int i = 0; i < M; ++i) {
gsl_vector_set(vb, i, b[i]);
for (int j = 0; j < N; ++j)
gsl_matrix_set(mA, i, j, A[i][j]);
}
rc = gsl_linalg_HH_svx(mA, vb);
printf("Solution Vector x:");
for (int i = 0; i < N; ++i)
printf("\t%0.6f", gsl_vector_get(vb, i));
printf("\n");
return rc;
}
--- linalg/hh_old.c 2014-05-30 12:57:40.232099826 +1000
+++ linalg/hh.c 2014-05-30 12:58:53.603105257 +1000
@@ -59,20 +59,20 @@
int
gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x)
{
- if (A->size1 > A->size2)
+ if (A->size2 > A->size1)
{
/* System is underdetermined. */
GSL_ERROR ("System is underdetermined", GSL_EINVAL);
}
- else if (A->size2 != x->size)
+ else if (A->size1 != x->size)
{
GSL_ERROR ("matrix and vector sizes must be equal", GSL_EBADLEN);
}
else
{
- const size_t N = A->size1;
- const size_t M = A->size2;
+ const size_t M = A->size1;
+ const size_t N = A->size2;
size_t i, j, k;
REAL *d = (REAL *) malloc (N * sizeof (REAL));
