[Help-gsl] Numerical Precision

[Thilo Schneider]

Dear all,

right now I'm struggling with my first use of GSL. Probably it is just  
a stupid mistake of mine, but nonetheless I am not able to resolve it.
I implemented an algorithm consisting of about 80 operations, mostly  
using linear algebra routines. As the results are not what I desired,  
I tried to track down the error by stepwise comparison between MATLAB  
and the results produced by GSL. Surprisingly, it is already the first  
matrix multiplication that leads to inaccurate results. What I do:

gsl_matrix *A = gsl_matrix_alloc(600, 250);
// initalize A. At this point A == matlabA up to working precision
// MATLAB: matlabB = matlabA' * matlabA
gsl_matrix *B = gsl_matrix_alloc(A->size2, A->size2);
gsl_matrix_set_zero(B);
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, A, A, 0.0, B);
// Now: max(abs(B - matlabB))  > 3e-10

Now these differences get bigger after every operation. I fear it is  
the datatype double which offers not enough precision for my  
computations.

Is there any way to increase GSLs precision, eg a custom datatype  
offering more accurate results?

Thank you for your help!

Kind regards,
Thilo Schneider