bug-gsl
[Top][All Lists]

## [Bug-gsl] gsl_blas_dsyrk in gsl-1.9

 From: Xiaogang Dong Subject: [Bug-gsl] gsl_blas_dsyrk in gsl-1.9 Date: Tue, 22 May 2007 15:57:55 -0700 User-agent: Thunderbird 1.5.0.10 (X11/20070304)

```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 .

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;
}
```

reply via email to