/* 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;

}