[Bug-gsl] bug in gsl_blas_dsyrk for versions 1.4 and 1.6

[Igor Nazarov]

Hi,
I have encountered a bug in gsl_blas_dsyrk() in GSL versions 1.4 and 
1.6 ( I tried it on different computers). Please find a test program 
below for calculating x*x' for a simple 3*2 matrix x.
--------------
In version 1.6 it produced an error (changing uplo to  CblasTrans fixed 
the error, but still produced incorrect result)
gsl: blas.c:1661: ERROR: invalid length
Default GSL error handler invoked.
Abort trap
--------------
In version 1.4 it produced the following output
1.000000 0.000000 0.000000
0.000000 5.000000 0.000000
1.000000 2.000000 2.000000
--------------
Actually, x*x'=[
1 0 1
0 4 2
1 2 2
Note that (2,2) entry should be 4.0 not 5.0!

That is just a simple test example showing that 2*2=5 for 
gsl_blas_dsyrk()
Trying different matrices x, one can see that only the last row(column) 
is correct in x*x' calculated using the dsyrk routine.
Also using gsl_blas_dgemm() to multiply x and x^T produces a correct 
result.

I would like to know whether this is fixed in newer releases (in 
version 1.9?).

Best Regards,
Igor.

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

int main()
{
double x[]={
1,0,
0,2,
1,1
};
int n=3;
double xtx[n*n];
int i,j;
gsl_matrix_view A = gsl_matrix_view_array(x, 3, 2);
gsl_matrix_view B = gsl_matrix_view_array(xtx, 3, 3);

gsl_blas_dsyrk (CblasUpper, CblasNoTrans,1.0, &A.matrix, 0.0, 
&B.matrix);
                for(i=0;i<n;i++)
                {
                        for(j=0;j<n;j++)
                                printf("%f ", xtx[i+j*n]);
                        printf("\n");
                }
}