/* spoper.c

 * 

 * Copyright (C) 2012 Patrick Alken

 * 

 * 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_matrix.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_spmatrix.h>

#include <gsl/gsl_spblas.h>

int

gsl_spmatrix_scale(gsl_spmatrix *m, const double x)

{

  size_t i;

  for (i = 0; i < m->nz; ++i)

    m->data[i] *= x;

  return GSL_SUCCESS;

} /* gsl_spmatrix_scale() */

int

gsl_spmatrix_minmax(const gsl_spmatrix *m, double *min_out, double *max_out)

{

  double min, max;

  size_t n;

  if (m->nz == 0)

    {

      GSL_ERROR("matrix is empty", GSL_EINVAL);

    }

  min = m->data[0];

  max = m->data[0];

  for (n = 1; n < m->nz; ++n)

    {

      double x = m->data[n];

      if (x < min)

        min = x;

      if (x > max)

        max = x;

    }

  *min_out = min;

  *max_out = max;

  return GSL_SUCCESS;

} /* gsl_spmatrix_minmax() */

/*

gsl_spmatrix_add()

  Add two sparse matrices

Inputs: c - (output) a + b

        a - (input) sparse matrix

        b - (input) sparse matrix

Return: success or error

*/

int

gsl_spmatrix_add(gsl_spmatrix *c, const gsl_spmatrix *a,

                 const gsl_spmatrix *b)

{

  const size_t M = a->size1;

  const size_t N = a->size2;

  if (b->size1 != M || b->size2 != N || c->size1 != M || c->size2 != N)

    {

      GSL_ERROR("matrices must have same dimensions", GSL_EBADLEN);

    }

  else if (a->sptype != b->sptype || a->sptype != c->sptype)

    {

      GSL_ERROR("matrices must have same sparse storage format",

                GSL_EINVAL);

    }

  else if (GSL_SPMATRIX_ISTRIPLET(a))

    {

      GSL_ERROR("triplet format not yet supported", GSL_EINVAL);

    }

  else

    {

      int status = GSL_SUCCESS;

      size_t *w = a->work_sze;

      double *x = c->work_dbl;

      size_t *Cp, *Ci;

      double *Cd;

      size_t j, p;

      size_t nz = 0; /* number of non-zeros in c */

      size_t inner_size, outer_size;

      if (GSL_SPMATRIX_ISCCS(a))

        {

          inner_size = M;

          outer_size = N;

        }

      else if (GSL_SPMATRIX_ISCRS(a))

        {

          inner_size = N;

          outer_size = M;

        }

      else

        {

          GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);

        }

      if (c->nzmax < a->nz + b->nz)

        {

          status = gsl_spmatrix_realloc(a->nz + b->nz, c);

          if (status)

            return status;

        }

      /* initialize w = 0 */

      for (j = 0; j < inner_size; ++j)

        w[j] = 0;

      Ci = c->i;

      Cp = c->p;

      Cd = c->data;

      for (j = 0; j < outer_size; ++j)

        {

          Cp[j] = nz;

          /* CCS: x += A(:,j); CRS: x += A(j,:) */

          nz = gsl_spblas_scatter(a, j, 1.0, w, x, j + 1, c, nz);

          /* CCS: x += B(:,j); CRS: x += B(j,:) */

          nz = gsl_spblas_scatter(b, j, 1.0, w, x, j + 1, c, nz);

          for (p = Cp[j]; p < nz; ++p)

            Cd[p] = x[Ci[p]];

        }

      /* finalize last column of c */

      Cp[j] = nz;

      c->nz = nz;

      return status;

    }

} /* gsl_spmatrix_add() */

/*

gsl_spmatrix_d2sp()

  Convert a dense gsl_matrix to sparse (triplet) format

Inputs: S - (output) sparse matrix in triplet format

        A - (input) dense matrix to convert

*/

int

gsl_spmatrix_d2sp(gsl_spmatrix *S, const gsl_matrix *A)

{

  int s = GSL_SUCCESS;

  size_t i, j;

  gsl_spmatrix_set_zero(S);

  S->size1 = A->size1;

  S->size2 = A->size2;

  for (i = 0; i < A->size1; ++i)

    {

      for (j = 0; j < A->size2; ++j)

        {

          double x = gsl_matrix_get(A, i, j);

          if (x != 0.0)

            gsl_spmatrix_set(S, i, j, x);

        }

    }

  return s;

} /* gsl_spmatrix_d2sp() */

/*

gsl_spmatrix_sp2d()

  Convert a sparse matrix to dense format

*/

int

gsl_spmatrix_sp2d(gsl_matrix *A, const gsl_spmatrix *S)

{

  if (A->size1 != S->size1 || A->size2 != S->size2)

    {

      GSL_ERROR("matrix sizes do not match", GSL_EBADLEN);

    }

  else

    {

      gsl_matrix_set_zero(A);

      if (GSL_SPMATRIX_ISTRIPLET(S))

        {

          size_t n;

          for (n = 0; n < S->nz; ++n)

            {

              size_t i = S->i[n];

              size_t j = S->p[n];

              double x = S->data[n];

              gsl_matrix_set(A, i, j, x);

            }

        }

      else

        {

          GSL_ERROR("non-triplet formats not yet supported", GSL_EINVAL);

        }

      return GSL_SUCCESS;

    }

} /* gsl_spmatrix_sp2d() */