Hi, while writing a program that uses exact diagonalization techniques I've encountered a strange problem with gsl_matrix_complex_fwrite() and its counterpart gsl_matrix_complex_fread(). Somehow the data is not stored correctly on the disk...
So, here's a short overview of the program: 1) Create a complex Hermitian matrix. 2) Diagonalize this matrix. 3) Store the eigensystem on the disk, using separate files for the eigenvectors and eigenvalues. In later steps I'll use the eigensystem to transform other complex matrices. Actually, there's quite a number of matrices, and the dimensions of the matrices can be of the order of 8000 by 8000. Because I cannot keep all of those matrices in memory, I want to store them on the disk temporarily. However, using the appropriate fwrite() and fread() methods is quite an issue, so it seems. After writing the matrix of eigenvectors to the disk and reloading them immediately afterwards, the matrix is useless, because rows are somehow missing then. I've written a small program that should demonstrate my problem. I kindly ask you to run this program and verify/falsify my observation that the matrix of eigenvectors before writing to and after reading from the disk differ. If you find any other GSL-related errors in my code, I'd be happy to know about them as well! Yours Matthias
run.sh
Description: application/shellscript
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_version.h>
// Sets matrix elements to predefined values depending on row and column
// indices. The resulting matrix is complex Hermitian.
int
set_matrix (gsl_matrix_complex * m, const int size)
{
int i, j;
for (i = 0; i < size; i++)
{
gsl_matrix_complex_set (m, i, i, (gsl_complex) {{ (i+1) * (i+1), 0}});
for (j = 0; j < i; j++)
{
gsl_matrix_complex_set (m, i, j, (gsl_complex) {{ (i+1) * (j+1), j-5 }});
gsl_matrix_complex_set (m, j, i, (gsl_complex) {{ (i+1) * (j+1), -j+5 }});
}
}
return 0;
}
// Prints the matrix line-by-line to stdout.
int
print_matrix (const gsl_matrix_complex * m, const int nrows, const int ncols,
const char * label)
{
int i, j;
gsl_complex z;
fprintf (stdout, "%s\n\n", label);
fprintf (stdout,
"tda = %i\nsize1 = %i\nsize2 = %i\n\n",
m->tda, m->size1, m->size2);
// gsl_matrix_complex_fprintf (stdout, m, "%f");
// fputc ('\n', stdout);
for (i = 0; i < nrows; i++)
{
for (j = 0; j < ncols; j++)
{
z = gsl_matrix_complex_get (m, i, j);
fprintf (stdout, "{% f, % f}", GSL_REAL(z), GSL_IMAG(z));
fputc ((j < ncols - 1) ? '\t' : '\n', stdout);
}
}
fputc ('\n', stdout);
return 0;
}
// Prints the vector on single line to stdout.
int
print_vector (const gsl_vector * v, const int nsize, const char * label)
{
int i;
fprintf (stdout, "%s\n\n", label);
// gsl_vector_fprintf (stdout, v, "%f");
// fputc ('\n', stdout);
for (i = 0; i < nsize; i++)
{
fprintf (stdout, "%f", gsl_vector_get (v, i));
fputc ((i < nsize - 1) ? '\t' : '\n', stdout);
}
fputc ('\n', stdout);
return 0;
}
// Prints the data block line-by-line assuming that matrix is complex-valued,
// but using double variables.
int
dump_matrix (const double * data, const int nrows, const int ncols,
const int tda, const char * label)
{
int i, j;
fprintf (stdout, "%s (dump)\n\n", label);
fprintf (stdout, "tda = %i\nsize1 = %i\nsize2 = %i\n\n", tda, nrows, ncols);
for (i = 0; i < 2 * nrows; i += 2)
{
for (j = 0; j < 2 * ncols; j += 2)
{
fprintf (stdout, "{% f, % f}",
data[i * tda + j], data[i * tda + j + 1]);
fputc ((j < 2 * (ncols - 1)) ? '\t' : '\n', stdout);
}
}
fputc ('\n', stdout);
return 0;
}
int
main (int argc, char ** argv)
{
int i, j;
gsl_complex z;
const int tda = 7;
const int size = 5;
// Print GSL version.
fprintf (stdout, "GSL v%s\n", gsl_version);
fputc ('\n', stdout);
gsl_matrix_complex * m = gsl_matrix_complex_calloc (tda, tda);
gsl_matrix_complex * evec = gsl_matrix_complex_calloc (tda, tda);
gsl_vector * eval = gsl_vector_calloc (tda);
m->size1 = m->size2 = size;
evec->size1 = evec->size2 = size;
eval->size = size;
// Fill matrix.
set_matrix (m, size);
// Print matrix.
print_matrix (m, size, size, "matrix");
// Diagonalize matrix.
gsl_eigen_hermv_workspace * ws = gsl_eigen_hermv_alloc (tda);
gsl_eigen_hermv (m, eval, evec, ws);
gsl_eigen_hermv_free (ws);
// Print diagonalized matrix.
print_matrix (m, size, size, "matrix after diagonalization routine");
// Print vector of eigenvalues.
print_vector (eval, size, "vector of eigenvalues");
// Print matrix of eigenvectors.
print_matrix (evec, size, size, "matrix of eigenvectors");
// Print full matrix of eigenvectors.
evec->size1 = evec->size2 = tda;
print_matrix (evec, tda, tda, "FULL matrix of eigenvectors");
evec->size1 = evec->size2 = size;
// Print matrix of eigenvector using doubles (just checking).
dump_matrix (evec->data, size, size, tda, "matrix of eigenvectors");
dump_matrix (evec->data, tda, tda, tda, "FULL matrix of eigenvectors");
// Restore matrix.
set_matrix (m, size);
// Similarity transformation: conjg(trans(evec)).m.evec
gsl_matrix_complex * tmp = gsl_matrix_complex_calloc (tda, tda);
tmp->size1 = tmp->size2 = size;
gsl_complex alpha = GSL_COMPLEX_ONE;
gsl_complex beta = GSL_COMPLEX_ZERO;
gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, alpha, m, evec, beta, tmp);
gsl_blas_zgemm (CblasConjTrans, CblasNoTrans, alpha, evec, tmp, beta, m);
gsl_matrix_complex_free (tmp);
// Print transformed matrix.
print_matrix (m, size, size, "matrix m after similarity transformation");
// Write eigenvectors to disk.
const char * filename = "evec.bin";
FILE * stream = fopen (filename, "w");
gsl_matrix_complex_fwrite (stream, evec);
fclose (stream);
// Check file size.
int filesize_est = size * size * sizeof (gsl_complex);
fprintf (stdout, "estimated size of evec file: %i byte\n", filesize_est);
stream = fopen (filename, "r");
fseek (stream, 0L, SEEK_END);
int filesize_real = ftell (stream);
fclose (stream);
fprintf (stdout, "real size of evec file: %i byte\n", filesize_real);
fputc ('\n', stdout);
// Read eigenvectors from disk.
gsl_matrix_complex_set_zero (evec);
stream = fopen (filename, "r");
gsl_matrix_complex_fread (stream, evec);
fclose (stream);
// Print matrix of eigenvectors.
print_matrix (evec, size, size, "matrix of eigenvectors after reading");
// Print full matrix of eigenvectors.
evec->size1 = evec->size2 = tda;
print_matrix (evec, tda, tda, "FULL matrix of eigenvectors after reading");
evec->size1 = evec->size2 = size;
// Brute-force check of the evec matrix on disk.
double * evec_bf = (double *) calloc (size * size, sizeof (gsl_complex));
stream = fopen (filename, "r");
int nread = fread (evec_bf, sizeof (gsl_complex), size * size, stream);
fclose (stream);
fprintf (stdout, "number of elements read: %i\n\n", nread);
dump_matrix (evec_bf, size, size, size, "matrix of eigenvectors after reading");
free (evec_bf);
gsl_matrix_complex_free (m);
gsl_matrix_complex_free (evec);
gsl_vector_free (eval);
return EXIT_SUCCESS;
}
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Help-gsl mailing list Help-gsl@gnu.org https://lists.gnu.org/mailman/listinfo/help-gsl
