Hi,
I've discovered a bug in the I/O for complex matrices. If tda is larger
than the actual size of the matrix, the fwrite() function called
internally writes the wrong blocks of data to disk. In particular, the
code containing the address of the memory block ("m->data + i * tda")
does /not/ respect the MULTIPLICITY!I've written a small program to demonstrate the bug with some small complex-valued matrix that is stored and read immediately afterwards from disk. I've also included a patch proposal for "matrix/file_source.c". I'd be happy if you could verify/falsify my findings. Best, Matthias
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_version.h>
int
main (int argc, char ** argv)
{
fprintf (stdout, "GSL version: %s\n\n", gsl_version);
int i, j;
FILE * os;
gsl_complex z;
const int tda = 5;
const int size = 3;
gsl_matrix_complex * m = gsl_matrix_complex_alloc (tda, tda);
m->size1 = m->size2 = size;
for (i = 0; i < size; i++)
{
for (j = 0; j < size; j++)
{
z = (gsl_complex) {{ i+j+1, -(i+j+1) }};
gsl_matrix_complex_set (m, i, j, z);
}
}
fputs ("original complex-valued matrix:\n", stdout);
//gsl_matrix_complex_fprintf (stdout, m, "%g");
for (i = 0; i < tda; i++) // Use this version; it provides low-level access.
{
for (j = 0; j < tda; j++)
{
fprintf (stdout, "{%g, %g}",
m->data[2 * (i * tda + j)],
m->data[2 * (i * tda + j) + 1]);
fputc ((j < tda - 1) ? '\t' : '\n', stdout);
}
}
fputc ('\n', stdout);
os = fopen ("matrix.bin", "w");
gsl_matrix_complex_fwrite (os, m);
fclose (os);
gsl_matrix_complex_set_zero (m);
os = fopen ("matrix.bin", "r");
gsl_matrix_complex_fread (os, m);
fclose (os);
fputs ("matrix after writing to and reading from disk:\n", stdout);
//gsl_matrix_complex_fprintf (stdout, m, "%g");
for (i = 0; i < tda; i++) // Use this version; it provides low-level access.
{
for (j = 0; j < tda; j++)
{
fprintf (stdout, "{%g, %g}",
m->data[2 * (i * tda + j)],
m->data[2 * (i * tda + j) + 1]);
fputc ((j < tda - 1) ? '\t' : '\n', stdout);
}
}
fputc ('\n', stdout);
gsl_matrix_complex_free (m);
return EXIT_SUCCESS;
}
43c43 < m->data + i * tda, --- > m->data + i * MULTIPLICITY * tda, 76c76 < m->data + i * tda, --- > m->data + i * MULTIPLICITY * tda, 113c113 < m->data + i * tda, --- > m->data + i * MULTIPLICITY * tda, 147c147 < m->data + i * tda, --- > m->data + i * MULTIPLICITY * tda,
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Bug-gsl mailing list [email protected] https://lists.gnu.org/mailman/listinfo/bug-gsl
