/* spio.c
*
* Copyright (C) 2016 Patrick Alken
* Copyright (C) 2016 Alexis Tantet
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>
/*
gsl_spmatrix_fprintf()
Print sparse matrix to file in MatrixMarket format:
M N NNZ
I1 J1 A(I1,J1)
...
Note that indices start at 1 and not 0
*/
int
gsl_spmatrix_fprintf(FILE *stream, const gsl_spmatrix *m,
const char *format)
{
int status;
/* print header */
status = fprintf(stream, "%%%%MatrixMarket matrix coordinate real general\n");
if (status < 0)
{
GSL_ERROR("fprintf failed for header", GSL_EFAILED);
}
/* print rows,columns,nnz */
status = fprintf(stream, "%u\t%u\t%u\n",
(unsigned int) m->size1,
(unsigned int) m->size2,
(unsigned int) m->nz);
if (status < 0)
{
GSL_ERROR("fprintf failed for dimension header", GSL_EFAILED);
}
if (GSL_SPMATRIX_ISTRIPLET(m))
{
size_t n;
for (n = 0; n < m->nz; ++n)
{
status = fprintf(stream, "%u\t%u\t", (unsigned int) m->i[n] + 1, (unsigned int) m->p[n] + 1);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = fprintf(stream, format, m->data[n]);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = putc('\n', stream);
if (status == EOF)
{
GSL_ERROR("putc failed", GSL_EFAILED);
}
}
}
else if (GSL_SPMATRIX_ISCCS(m))
{
size_t j, p;
for (j = 0; j < m->size2; ++j)
{
for (p = m->p[j]; p < m->p[j + 1]; ++p)
{
status = fprintf(stream, "%u\t%u\t", (unsigned int) m->i[p] + 1, (unsigned int) j + 1);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = fprintf(stream, format, m->data[p]);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = putc('\n', stream);
if (status == EOF)
{
GSL_ERROR("putc failed", GSL_EFAILED);
}
}
}
}
else if (GSL_SPMATRIX_ISCRS(m))
{
size_t i, p;
for (i = 0; i < m->size1; ++i)
{
for (p = m->p[i]; p < m->p[i + 1]; ++p)
{
status = fprintf(stream, "%u\t%u\t", (unsigned int) i + 1, (unsigned int) m->i[p] + 1);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = fprintf(stream, format, m->data[p]);
if (status < 0)
{
GSL_ERROR("fprintf failed", GSL_EFAILED);
}
status = putc('\n', stream);
if (status == EOF)
{
GSL_ERROR("putc failed", GSL_EFAILED);
}
}
}
}
else
{
GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
}
return GSL_SUCCESS;
}
gsl_spmatrix *
gsl_spmatrix_fscanf(FILE *stream)
{
gsl_spmatrix *m;
unsigned int size1, size2, nz;
char buf[1024];
int found_header = 0;
/* read file until we find rows,cols,nz header */
while (fgets(buf, 1024, stream) != NULL)
{
int c;
/* skip comments */
if (*buf == '%')
continue;
c = sscanf(buf, "%u %u %u",
&size1, &size2, &nz);
if (c == 3)
{
found_header = 1;
break;
}
}
if (!found_header)
{
GSL_ERROR_NULL ("fscanf failed reading header", GSL_EFAILED);
}
m = gsl_spmatrix_alloc_nzmax((size_t) size1, (size_t) size2, (size_t) nz, GSL_SPMATRIX_TRIPLET);
if (!m)
{
GSL_ERROR_NULL ("error allocating m", GSL_ENOMEM);
}
{
unsigned int i, j;
double val;
while (fgets(buf, 1024, stream) != NULL)
{
int c = sscanf(buf, "%u %u %lg", &i, &j, &val);
if (c < 3 || (i == 0) || (j == 0))
{
GSL_ERROR_NULL ("error in input file format", GSL_EFAILED);
}
else if ((i > size1) || (j > size2))
{
GSL_ERROR_NULL ("element exceeds matrix dimensions", GSL_EBADLEN);
}
else
{
/* subtract 1 from (i,j) since indexing starts at 1 */
gsl_spmatrix_set(m, i - 1, j - 1, val);
}
}
}
return m;
}
int
gsl_spmatrix_fwrite(FILE *stream, const gsl_spmatrix *m)
{
size_t items;
/* write header: size1, size2, nz */
items = fwrite(&(m->size1), sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fwrite failed on size1", GSL_EFAILED);
}
items = fwrite(&(m->size2), sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fwrite failed on size2", GSL_EFAILED);
}
items = fwrite(&(m->nz), sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fwrite failed on nz", GSL_EFAILED);
}
/* write m->i and m->data which are size nz in all storage formats */
items = fwrite(m->i, sizeof(size_t), m->nz, stream);
if (items != m->nz)
{
GSL_ERROR("fwrite failed on row indices", GSL_EFAILED);
}
items = fwrite(m->data, sizeof(double), m->nz, stream);
if (items != m->nz)
{
GSL_ERROR("fwrite failed on data", GSL_EFAILED);
}
if (GSL_SPMATRIX_ISTRIPLET(m))
{
items = fwrite(m->p, sizeof(size_t), m->nz, stream);
if (items != m->nz)
{
GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
}
}
else if (GSL_SPMATRIX_ISCCS(m))
{
items = fwrite(m->p, sizeof(size_t), m->size2 + 1, stream);
if (items != m->size2 + 1)
{
GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
}
}
else if (GSL_SPMATRIX_ISCRS(m))
{
items = fwrite(m->p, sizeof(size_t), m->size1 + 1, stream);
if (items != m->size1 + 1)
{
GSL_ERROR("fwrite failed on column indices", GSL_EFAILED);
}
}
return GSL_SUCCESS;
}
int
gsl_spmatrix_fread(FILE *stream, gsl_spmatrix *m)
{
size_t size1, size2, nz;
size_t items;
/* read header: size1, size2, nz */
items = fread(&size1, sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fread failed on size1", GSL_EFAILED);
}
items = fread(&size2, sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fread failed on size2", GSL_EFAILED);
}
items = fread(&nz, sizeof(size_t), 1, stream);
if (items != 1)
{
GSL_ERROR("fread failed on nz", GSL_EFAILED);
}
if (m->size1 != size1)
{
GSL_ERROR("matrix has wrong size1", GSL_EBADLEN);
}
else if (m->size2 != size2)
{
GSL_ERROR("matrix has wrong size2", GSL_EBADLEN);
}
else if (nz > m->nzmax)
{
GSL_ERROR("matrix nzmax is too small", GSL_EBADLEN);
}
else
{
/* read m->i and m->data arrays, which are size nz for all formats */
items = fread(m->i, sizeof(size_t), nz, stream);
if (items != nz)
{
GSL_ERROR("fread failed on row indices", GSL_EFAILED);
}
items = fread(m->data, sizeof(double), nz, stream);
if (items != nz)
{
GSL_ERROR("fread failed on data", GSL_EFAILED);
}
m->nz = nz;
if (GSL_SPMATRIX_ISTRIPLET(m))
{
items = fread(m->p, sizeof(size_t), nz, stream);
if (items != nz)
{
GSL_ERROR("fread failed on column indices", GSL_EFAILED);
}
/* build binary search tree for m */
gsl_spmatrix_tree_rebuild(m);
}
else if (GSL_SPMATRIX_ISCCS(m))
{
items = fread(m->p, sizeof(size_t), size2 + 1, stream);
if (items != size2 + 1)
{
GSL_ERROR("fread failed on row pointers", GSL_EFAILED);
}
}
else if (GSL_SPMATRIX_ISCRS(m))
{
items = fread(m->p, sizeof(size_t), size1 + 1, stream);
if (items != size1 + 1)
{
GSL_ERROR("fread failed on column pointers", GSL_EFAILED);
}
}
}
return GSL_SUCCESS;
}