[Bug-gsl] gsl_blas_dsyrk in gsl-1.9

[Xiaogang Dong]

To whom may concern,

I tried to use gsl_blas_dsyrk to calculate the matrix product
A^TA with A =
[ 1 2
 1 3
 1 4
 1 5 ].
The results produced by gsl_blas_dsyrk seems incorrect.
Bying trying different matrices, I figured out it only used
the first n rows of an m x n matrix ( assume m > n ).
For example, instead of calculating A^TA, it actually
calculate B^TB with B =
[ 1 2
 1 3 ].

My computer is a Dell XPS 410 with Cent OS 4.0 .
Please let me know if you need more info.

Thanks,
Xiaogang

#include <stdio.h>
#include "../gsl-1.9/gsl/gsl_blas.h"

int main()
{
  int m = 4, n = 2;
  gsl_matrix *a = gsl_matrix_alloc( m, n );

  int i, j;
  for ( i = 0; i < m; i++ ) 
    for ( j = 0; j < n; j++ )
      gsl_matrix_set( a, i, j, ( double )( i * j + j + 1 ) );

  for ( i = 0; i < m; i++ ) {
    for ( j = 0; j < n; j++ )
      printf( "%.5f ", gsl_matrix_get( a, i, j ) );
    printf( "\n" );
  }
  printf( "\n" );

  gsl_matrix *v = gsl_matrix_calloc( n, n ); 
  gsl_blas_dsyrk( CblasLower, CblasTrans, 1.0, a, 0.0, v );

  for ( i = 0; i < n; i++ ) {
    for ( j = 0; j < n; j++ )
      printf( "%.5f ", gsl_matrix_get( v, i, j ) );
    printf( "\n" );
  }
  printf( "\n" );

  gsl_matrix_free( v );
  gsl_matrix_free( a );

  return 0;
}