[Help-gsl] AX=B octave vs gsl

[Awhan Patnaik]

Hello all

I am trying to solve for X in AX=B where A, B, X are all matrices.
A is a symmetric 15x15 matrix
X is the unknown 15x30 matrix
B is of course 15x30 matrix.

condition number c of A calculated in Octave is 3.2532e+06.

A matrix data can be obtained here
http://dl.dropbox.com/u/842054/gsl%20vs%20octave%20LU%20decomposition/a.in
Can be loaded in octave using
load a.in

B matrix data can be obtained here
http://dl.dropbox.com/u/842054/gsl%20vs%20octave%20LU%20decomposition/b.in
Can be loaded in octave using
load b.in

I solved for X using LU method with refinement using this code:
http://dl.dropbox.com/u/842054/gsl%20vs%20octave%20LU%20decomposition/main.cpp

Please ignore the read() and write() function in my C++ file as well
as utils.h. I show the C++ file just to show the method used and get
verification if I am doing it right.

I also solve this problem in Octave and am surprised at the great
accuracy obtained in Octave. Let X_gsl be the solution obtained using
LU in GSL and let X_oct be the octave solution obtained using A\B

norm(A*X_oct - B) = 6.6081e-11

norm(A*X_gsl - B) = 0.13395

Aslo I find the residuals calculated in the LU refinement step are of
the order of 1e-7 so not much improvement actually occurs.

Octave solution can be obtained here
http://dl.dropbox.com/u/842054/gsl%20vs%20octave%20LU%20decomposition/x.oct


GSL+LU solution can be obtained here
http://dl.dropbox.com/u/842054/gsl%20vs%20octave%20LU%20decomposition/x.gsl


What can I do to get as much accuracy as the Octave method?