LU solver results differing from numpy's

[MaximeLee]

I want to use gsl LU solver. But when comparing with Scipy's LU solver the results are a different. The zip archive contains the source main.c file that solves the linear sytem Ax = b with A the matrix in the A.dat file written using gsl_matrix_fwrite and b the vector whose only non-zero component is the last with value -1. Finally the test.py file print scipy's solution.

The solutions I get are with
GSL:

-0.944006998359841698054140124441 
2.689851515690580452400126887369 
-0.745844517330738865368289225444 
-0.000000000000231038487010262759

Scipy:

-5.40764393e-01  
2.99927045e+00 
-1.45850606e+00 
-2.47812410e-13

gsl_lu.zip

[rdguerrerom]

gsl_linalg_LU_decomp and gsl_linalg_LU_solve are fine. The issue is almost certainly that the matrix is being read incorrectly from disk.

gsl_matrix_fwrite() writes the matrix in binary format, not text. If your A.dat was generated with:

gsl_matrix_fwrite(f, A);

then it must be read back with:

gsl_matrix_fread(f, A);

and the file must be opened in binary mode:

FILE *f = fopen("A.dat", "rb");

If instead you are using something like fscanf() or opening with "r", the matrix contents become corrupted, which explains why LU gives a completely different solution.

The correct pattern is:

FILE *f = fopen("A.dat", "rb");

gsl_matrix *A = gsl_matrix_alloc(n, n);

gsl_matrix_fread(f, A);

fclose(f);

Then solve with:

gsl_permutation *p = gsl_permutation_alloc(n);
int signum;

gsl_linalg_LU_decomp(A, p, &signum);
gsl_linalg_LU_solve(A, p, b, x);

Your SciPy result is likely the correct one:

-5.40764393e-01
 2.99927045e+00
-1.45850606e+00
-2.47812410e-13

because the GSL output you showed is consistent with solving a different matrix.

Also verify that:

• the matrix dimensions match exactly,
• you are using double precision everywhere,
• the vector b is initialized identically in both codes,
• and you are not accidentally overwriting A before solving.

For reference, LU decomposition solves:

Ax=b

and the decomposition computed internally is:

PA=LU

where P is the permutation matrix from pivoting.