/* -*- mode: C -*-  */
/* 
   IGraph library.
   Copyright (C) 2003, 2004, 2005  Gabor Csardi <csardi@rmki.kfki.hu>
   MTA RMKI, Konkoly-Thege Miklos st. 29-33, Budapest 1121, Hungary
   
   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 2 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 "types.h"
#include "memory.h"
#include "random.h"
#include "error.h"

#include <assert.h>
#include <string.h>     /* memcpy & co. */
#include <stdlib.h>

/**
 * \section igraph_spmatrix_constructor_and_destructor Sparse matrix constructors
 * and destructors.
 */

/**
 * \ingroup matrix
 * \function igraph_spmatrix_init
 * \brief Initializes a sparse matrix.
 * 
 * </para><para>
 * Every sparse matrix needs to be initialized before using it, this is done
 * by calling this function. A matrix has to be destroyed if it is not
 * needed any more, see \ref igraph_spmatrix_destroy().
 * \param m Pointer to a not yet initialized sparse matrix object to be
 *        initialized. 
 * \param nrow The number of rows in the matrix.
 * \param ncol The number of columns in the matrix.
 * \return Error code.
 *
 * Time complexity: operating system dependent.
 */

int igraph_spmatrix_init(igraph_spmatrix_t *m, long int nrow, long int ncol) {
  assert(m != NULL);
  IGRAPH_VECTOR_INIT_FINALLY(&m->ridx, 0);
  IGRAPH_VECTOR_INIT_FINALLY(&m->cidx, ncol+1);
  IGRAPH_VECTOR_INIT_FINALLY(&m->data, 0);
  IGRAPH_FINALLY_CLEAN(3);
  m->nrow=nrow;
  m->ncol=ncol;
  return 0;
}

/** 
 * \ingroup matrix
 * \function igraph_spmatrix_destroy
 * \brief Destroys a sparse matrix object.
 * 
 * </para><para>
 * This function frees all the memory allocated for a sparse matrix
 * object. The destroyed object needs to be reinitialized before using
 * it again.
 * \param m The matrix to destroy.
 * 
 * Time complexity: operating system dependent.
 */ 

void igraph_spmatrix_destroy(igraph_spmatrix_t *m) {
  assert(m != NULL);
  igraph_vector_destroy(&m->ridx);
  igraph_vector_destroy(&m->cidx);
  igraph_vector_destroy(&m->data);
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_copy
 * \brief Copies a sparse matrix.
 *
 * </para><para>
 * Creates a sparse matrix object by copying another one.
 * \param to Pointer to an uninitialized sparse matrix object.
 * \param from The initialized sparse matrix object to copy.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *   isn't enough memory to allocate the new sparse matrix.
 * 
 * Time complexity: O(n), the number
 * of elements in the matrix.
 */

int igraph_spmatrix_copy(igraph_spmatrix_t *to, const igraph_spmatrix_t *from) {
  assert(from != NULL);
  assert(to != NULL);
  to->nrow = from->nrow;
  to->ncol = from->ncol;
  IGRAPH_CHECK(igraph_vector_copy(&to->ridx, &from->ridx));
  IGRAPH_CHECK(igraph_vector_copy(&to->cidx, &from->cidx));
  IGRAPH_CHECK(igraph_vector_copy(&to->data, &from->data));
  return 0;
}

/**
 * \section igraph_spmatrix_accessing_elements Accessing elements of a sparse matrix
 */

/**
 * \ingroup matrix
 * \function igraph_spmatrix_e
 * \brief Accessing an element of a sparse matrix.
 *
 * Note that there are no range checks right now. 
 * \param m The matrix object.
 * \param row The index of the row, starting with zero.
 * \param col The index of the column, starting with zero.
 *
 * Time complexity: O(log n), where n is the number of nonzero elements in
 * the requested column.
 */
igraph_real_t igraph_spmatrix_e(const igraph_spmatrix_t *m,
                                long int row, long int col) {
  long int start, end;
 
  assert(m != NULL);
  start=(long)igraph_vector_e(&m->cidx, col);
  end=(long)igraph_vector_e(&m->cidx, col+1)-1;

  if (end<start) return 0; 
  /* Elements residing in column col are between m->data[start] and
   * m->data[end], inclusive, ordered by row index */
  while (start < end-1) {
    long int mid=(start+end)/2;
    if (VECTOR(m->ridx)[mid] > row) {
      end=mid;
    } else if (VECTOR(m->ridx)[mid] < row) {
      start=mid;
    } else {
      start=mid;
      break;
    }
  }

  if (VECTOR(m->ridx)[start] == row)
    return VECTOR(m->data)[start];
  if (VECTOR(m->ridx)[start] != row && VECTOR(m->ridx)[end] == row)
    return VECTOR(m->data)[end];
  return 0;
}


/**
 * \ingroup matrix
 * \function igraph_spmatrix_set
 * \brief Setting an element of a sparse matrix.
 *
 * Note that there are no range checks right now. 
 * \param m The matrix object.
 * \param row The index of the row, starting with zero.
 * \param col The index of the column, starting with zero.
 * \param value The new value.
 *
 * Time complexity: O(log n), where n is the number of nonzero elements in
 * the requested column.
 */
int igraph_spmatrix_set(igraph_spmatrix_t *m, long int row, long int col,
            igraph_real_t value) {
  long int start, end;
 
  assert(m != NULL);
  start=(long)igraph_vector_e(&m->cidx, col);
  end=(long)igraph_vector_e(&m->cidx, col+1)-1;

  if (end<start) {
    /* First element in the column */
    if (value == 0.0) return 0;
    IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start, row));
    IGRAPH_CHECK(igraph_vector_insert(&m->data, start, value));
    for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]++;
    return 0;
  }

  /* Elements residing in column col are between m->data[start] and
   * m->data[end], inclusive, ordered by row index */
  while (start < end-1) {
    long int mid=(start+end)/2;
    if (VECTOR(m->ridx)[mid] > row) { 
      end=mid;
    } else if (VECTOR(m->ridx)[mid] < row) {
      start=mid;
    } else {
      start=mid;
      break;
    }
  }

  if (VECTOR(m->ridx)[start] == row) {
    /* Overwriting a value - or deleting it if it has been overwritten by zero */
    if (value == 0) {
      igraph_vector_remove(&m->ridx, start);
      igraph_vector_remove(&m->data, start);
      for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]--;
    } else {
      VECTOR(m->data)[start] = value;
    }
    return 0;
  } else if (VECTOR(m->ridx)[end] == row) {
    /* Overwriting a value - or deleting it if it has been overwritten by zero */
    if (value == 0) {
      igraph_vector_remove(&m->ridx, end);
      igraph_vector_remove(&m->data, end);
      for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]--;
    } else {
      VECTOR(m->data)[end] = value;
    }
    return 0;
  }
  
  /* New element has to be inserted, but only if not a zero is
   * being written into the matrix */
  if (value != 0.0) {
    if (VECTOR(m->ridx)[end] < row) {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, end+1, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, end+1, value));
    } else if (VECTOR(m->ridx)[start] < row) {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start+1, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, start+1, value));
    } else {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, start, value));
    }
    for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]++;
  }
  return 0;
}


/**
 * \ingroup matrix
 * \function igraph_spmatrix_add_e
 * \brief Adding a real value to an element of a sparse matrix.
 *
 * Note that there are no range checks right now. This is implemented to avoid
 * double lookup of a given element in the matrix by using \ref igraph_spmatrix_e()
 * and \ref igraph_spmatrix_set() consecutively.
 *
 * \param m The matrix object.
 * \param row The index of the row, starting with zero.
 * \param col The index of the column, starting with zero.
 * \param value The value to add.
 *
 * Time complexity: O(log n), where n is the number of nonzero elements in
 * the requested column.
 */
int igraph_spmatrix_add_e(igraph_spmatrix_t *m, long int row, long int col,
            igraph_real_t value) {
  long int start, end;
 
  assert(m != NULL);
  start=(long)igraph_vector_e(&m->cidx, col);
  end=(long)igraph_vector_e(&m->cidx, col+1)-1;

  if (end<start) {
    /* First element in the column */
    if (value == 0.0) return 0;
    IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start, row));
    IGRAPH_CHECK(igraph_vector_insert(&m->data, start, value));
    for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]++;
    return 0;
  }

  /* Elements residing in column col are between m->data[start] and
   * m->data[end], inclusive, ordered by row index */
  while (start < end-1) {
    long int mid=(start+end)/2;
    if (VECTOR(m->ridx)[mid] > row) { 
      end=mid;
    } else if (VECTOR(m->ridx)[mid] < row) {
      start=mid;
    } else {
      start=mid;
      break;
    }
  }

  if (VECTOR(m->ridx)[start] == row) {
    /* Overwriting a value */
    if (VECTOR(m->data)[start] == -1) {
      igraph_vector_remove(&m->ridx, start);
      igraph_vector_remove(&m->data, start);
      for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]--;
    } else {
      VECTOR(m->data)[start] += value;
    }
    return 0;
  } else if (VECTOR(m->ridx)[end] == row) {
    /* Overwriting a value */
    if (VECTOR(m->data)[end] == -1) {
      igraph_vector_remove(&m->ridx, end);
      igraph_vector_remove(&m->data, end);
      for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]--;
    } else {
      VECTOR(m->data)[end] += value;
    }
    return 0;
  }
  
  /* New element has to be inserted, but only if not a zero is
   * being added to a zero element of the matrix */
  if (value != 0.0) {
    if (VECTOR(m->ridx)[end] < row) {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, end+1, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, end+1, value));
    } else if (VECTOR(m->ridx)[start] < row) {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start+1, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, start+1, value));
    } else {
      IGRAPH_CHECK(igraph_vector_insert(&m->ridx, start, row));
      IGRAPH_CHECK(igraph_vector_insert(&m->data, start, value));
    }
    for (start=col+1; start < m->ncol+1; start++) VECTOR(m->cidx)[start]++;
  }
  return 0;
}

/**
 * \function igraph_spmatrix_add_col_values
 *
 * Adds the values of a column to another column.
 *
 * \param to The index of the column to be added to
 * \param from The index of the column to be added
 * \return Error code.
 */
int igraph_spmatrix_add_col_values(igraph_spmatrix_t *m, long int to, long int from) {
  long int i;
  /* TODO: I think this implementation could be speeded up if I don't use
   * igraph_spmatrix_add_e directly -- but maybe it's not worth the fuss */
  for (i=VECTOR(m->cidx)[from]; i<VECTOR(m->cidx)[from+1]; i++) {
    IGRAPH_CHECK(igraph_spmatrix_add_e(m, VECTOR(m->ridx)[i], to, VECTOR(m->data)[i]));
  }

  return 0;
}


/**
 * \ingroup matrix
 * \function igraph_spmatrix_resize
 * \brief Resizes a sparse matrix.
 *
 * </para><para>
 * This function resizes a sparse matrix by adding more elements to it.
 * The matrix retains its data even after resizing it, except for the data
 * which lies outside the new boundaries (if the new size is smaller).
 * \param m Pointer to an already initialized sparse matrix object.
 * \param nrow The number of rows in the resized matrix.
 * \param ncol The number of columns in the resized matrix.
 * \return Error code.
 * 
 * Time complexity: O(n).
 * n is the number of elements in the old matrix.
 */

int igraph_spmatrix_resize(igraph_spmatrix_t *m, long int nrow, long int ncol) {
  long int i, j, ci, ei, mincol, nelem;
  assert(m != NULL);
  /* Iterating through the matrix data and deleting unnecessary data. */
  /* At the same time, we create the new indices as well */
  if (nrow < m->nrow) {
    nelem = igraph_vector_size(&m->data);
    ei = j = 0;
    mincol = (m->ncol < ncol) ? m->ncol : ncol;
    for (ci=0; ci < mincol; ci++) {
      for (; ei<VECTOR(m->cidx)[ci+1]; ei++) {
        if (VECTOR(m->ridx)[ei] < nrow) {
          VECTOR(m->ridx)[j]=VECTOR(m->ridx)[ei];
          VECTOR(m->data)[j]=VECTOR(m->data)[ei];
          j++;
        }
      }
      VECTOR(m->cidx)[ci]=j;
    }
    /* Contract the row index and the data vector */
    IGRAPH_CHECK(igraph_vector_resize(&m->ridx, j));
    IGRAPH_CHECK(igraph_vector_resize(&m->cidx, j));
  }
  /* Updating cidx */
  IGRAPH_CHECK(igraph_vector_resize(&m->cidx, ncol+1));
  for (i=m->ncol+1; i<ncol+1; i++) 
    VECTOR(m->cidx)[i] = VECTOR(m->cidx)[m->ncol];
  m->nrow=nrow;
  m->ncol=ncol;
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_count_nonzero
 * \brief The number of non-zero elements in a sparse matrix.
 * 
 * \param m Pointer to an initialized sparse matrix object.
 * \return The size of the matrix.
 *
 * Time complexity: O(1).
 */

long int igraph_spmatrix_count_nonzero(const igraph_spmatrix_t *m) {
  assert(m != NULL);
  return igraph_vector_size(&m->data);
}


/**
 * \ingroup matrix
 * \function igraph_spmatrix_size
 * \brief The number of elements in a sparse matrix.
 * 
 * \param m Pointer to an initialized sparse matrix object.
 * \return The size of the matrix.
 *
 * Time complexity: O(1).
 */

long int igraph_spmatrix_size(const igraph_spmatrix_t *m) {
  assert(m != NULL);
  return (m->nrow) * (m->ncol);
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_nrow
 * \brief The number of rows in a sparse matrix.
 * 
 * \param m Pointer to an initialized sparse matrix object.
 * \return The number of rows in the matrix.
 * 
 * Time complexity: O(1).
 */

long int igraph_spmatrix_nrow(const igraph_spmatrix_t *m) {
  assert(m != NULL);
  return m->nrow;
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_ncol
 * \brief The number of columns in a sparse matrix.
 * 
 * \param m Pointer to an initialized sparse matrix object.
 * \return The number of columns in the sparse matrix.
 * 
 * Time complexity: O(1).
 */

long int igraph_spmatrix_ncol(const igraph_spmatrix_t *m) {
  assert(m != NULL);
  return m->ncol;
}

/** 
 * \ingroup matrix
 * \brief Copies a sparse matrix to a regular C array.
 *
 * </para><para>
 * The matrix is copied columnwise, as this is the format most
 * programs and languages use.
 * The C array should be of sufficient size, there are (of course) no
 * range checks done.
 * \param m Pointer to an initialized sparse matrix object.
 * \param to Pointer to a C array, the place to copy the data to.
 * \return Error code.
 *
 * Time complexity: O(n),
 * n is the number of 
 * elements in the matrix.
 */

int igraph_spmatrix_copy_to(const igraph_spmatrix_t *m, igraph_real_t *to) {
  long int c, dest_idx, idx, no_of_elements;

  no_of_elements=igraph_spmatrix_count_nonzero(m);
  memset(to, 0, sizeof(igraph_real_t)*igraph_spmatrix_size(m));
  for (c=0, dest_idx=0; c < m->ncol; c++, dest_idx+=m->nrow) {
    for (idx=VECTOR(m->cidx)[c]; idx<VECTOR(m->cidx)[c+1]; idx++) {
      to[dest_idx+(long)VECTOR(m->ridx)[idx]]=VECTOR(m->data)[idx];
    }
  }
  return 0;
}

/** 
 * \ingroup matrix
 * \brief Sets all element in a sparse matrix to zero.
 * 
 * \param m Pointer to an initialized matrix object.
 * \return Error code, always returns with success.
 * 
 * Time complexity: O(n),
 * n is the number of columns in the matrix
 */

int igraph_spmatrix_null(igraph_spmatrix_t *m) {
  assert(m != NULL);
  igraph_vector_clear(&m->data);
  igraph_vector_clear(&m->ridx);
  igraph_vector_null(&m->cidx);
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_add_cols
 * \brief Adds columns to a sparse matrix.
 * \param m The sparse matrix object.
 * \param n The number of columns to add.
 * \return Error code.
 *
 * Time complexity: O(1).
 */

int igraph_spmatrix_add_cols(igraph_spmatrix_t *m, long int n) {
  igraph_spmatrix_resize(m, m->nrow, m->ncol+n);
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_spmatrix_add_rows
 * \brief Adds rows to a sparse matrix.
 * \param m The sparse matrix object.
 * \param n The number of rows to add.
 * \return Error code.
 * 
 * Time complexity: O(1).
 */

int igraph_spmatrix_add_rows(igraph_spmatrix_t *m, long int n) {
  igraph_spmatrix_resize(m, m->nrow+n, m->ncol);
  return 0;
}

/**
 * \function igraph_spmatrix_clear_row
 *
 * Clears a row in the matrix (sets all of its elements to zero)
 * \param m The matrix.
 * \param row The index of the row to be cleared.
 *
 * Time complexity: O(n), the number of nonzero elements in the matrix.
 */

int igraph_spmatrix_clear_row(igraph_spmatrix_t *m, long int row) {
  long int ci, ei, i, j, nremove=0, nremove_old=0, *permvec;

  assert(m != NULL);
  permvec=Calloc(igraph_vector_size(&m->data), long int);
  if (permvec == 0) {
    IGRAPH_ERROR("can't clear row in sparse matrix", IGRAPH_ENOMEM);
  }
  IGRAPH_FINALLY(free, permvec);
  for (ci=0, i=0, j=1; ci < m->ncol; ci++) {
    for (ei=VECTOR(m->cidx)[ci]; ei < VECTOR(m->cidx)[ci+1]; ei++) {
      if (VECTOR(m->ridx)[ei] == row) {
        /* this element will be deleted, so all elements in cidx from the
         * column index of this element will have to be decreased by one */
        nremove++;
      } else {
        /* this element will be kept */
        permvec[i] = j;
        j++;
      }
      i++;
    }
    if (ci > 0) {
      VECTOR(m->cidx)[ci] -= nremove_old;
    }
    nremove_old = nremove;
  }
  VECTOR(m->cidx)[m->ncol] -= nremove;
  igraph_vector_permdelete(&m->ridx, permvec, nremove);
  igraph_vector_permdelete(&m->data, permvec, nremove);
  free(permvec);
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

int igraph_i_spmatrix_clear_row_fast(igraph_spmatrix_t *m, long int row) {
  long int ei, n;

  assert(m != NULL);
  n = igraph_vector_size(&m->data);
  for (ei=0; ei<n; ei++) {
    if (VECTOR(m->ridx)[ei] == row) VECTOR(m->data)[ei]=0.0;
  }
  return 0;
}

int igraph_i_spmatrix_cleanup(igraph_spmatrix_t *m) {
  long int ci, ei, i, j, nremove=0, nremove_old=0, *permvec;

  assert(m != NULL);
  permvec=Calloc(igraph_vector_size(&m->data), long int);
  if (permvec == 0) {
    IGRAPH_ERROR("can't perform cleanup on sparse matrix", IGRAPH_ENOMEM);
  }
  IGRAPH_FINALLY(free, permvec);
  for (ci=0, i=0, j=1; ci < m->ncol; ci++) {
    for (ei=VECTOR(m->cidx)[ci]; ei < VECTOR(m->cidx)[ci+1]; ei++) {
      if (VECTOR(m->data)[ei] == 0.0) {
        /* this element will be deleted, so all elements in cidx from the
         * column index of this element will have to be decreased by one */
        nremove++;
      } else {
        /* this element will be kept */
        permvec[i] = j;
        j++;
      }
      i++;
    }
    if (ci > 0) {
      VECTOR(m->cidx)[ci] -= nremove_old;
    }
    nremove_old = nremove;
  }
  VECTOR(m->cidx)[m->ncol] -= nremove;
  igraph_vector_permdelete(&m->ridx, permvec, nremove);
  igraph_vector_permdelete(&m->data, permvec, nremove);
  free(permvec);
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

/**
 * \function igraph_spmatrix_clear_col
 *
 * Clears a column in the matrix (sets all of its elements to zero)
 * \param m The matrix.
 * \param col The index of the column to be cleared.
 * \return Error code. The current implementation always succeeds.
 *
 * Time complexity: TODO
 */

int igraph_spmatrix_clear_col(igraph_spmatrix_t *m, long int col) {
  long int i, n;
  assert(m != NULL);
  n = (long)VECTOR(m->cidx)[col+1] - (long)VECTOR(m->cidx)[col];
  if (n == 0) return 0;
  igraph_vector_remove_section(&m->ridx, VECTOR(m->cidx)[col], VECTOR(m->cidx)[col+1]);
  igraph_vector_remove_section(&m->data, VECTOR(m->cidx)[col], VECTOR(m->cidx)[col+1]);
  for (i=col+1; i <= m->ncol; i++) {
    VECTOR(m->cidx)[i] -= n;
  }
  return 0;
}

/**
 * \function igraph_spmatrix_multiply
 * 
 * Multiplies each element of the sparse matrix by a constant.
 * \param m The matrix.
 * \param by The constant.
 * 
 * Time complexity: O(n), the number of elements in the matrix.
 */

void igraph_spmatrix_multiply(igraph_spmatrix_t *m, igraph_real_t by) {
  assert(m != NULL);
  igraph_vector_multiply(&m->data, by);
}

/**
 * \function igraph_spmatrix_colsums
 *
 * Calculates the column sums of the matrix.
 * \param m The matrix.
 * \param res An initialized \c igraph_vector_t, the result will be stored here.
 *   The vector will be resized as needed.
 *
 * Time complexity: O(n), the number of nonzero elements in the matrix.
 */

int igraph_spmatrix_colsums(const igraph_spmatrix_t *m, igraph_vector_t *res) {
  long int i, c;
  assert(m != NULL);
  IGRAPH_CHECK(igraph_vector_resize(res, m->ncol));
  igraph_vector_null(res);
  for (c=0; c < m->ncol; c++) {
    for (i=VECTOR(m->cidx)[c]; i<VECTOR(m->cidx)[c+1]; i++) {
      VECTOR(*res)[c] += VECTOR(m->data)[i];
    }
  }
  return 0;
}

/**
 * \function igraph_spmatrix_max_nonzero
 *
 * Returns the maximum nonzero element of a matrix.
 * If the matrix is empty, zero is returned.
 *
 * \param m the matrix object.
 * \param ridx the row index of the maximum element if not \c NULL.
 * \param cidx the column index of the maximum element if not \c NULL.
 *
 * Time complexity: O(n), the number of nonzero elements in the matrix.
 */
igraph_real_t igraph_spmatrix_max_nonzero(const igraph_spmatrix_t *m,
    igraph_real_t *ridx, igraph_real_t *cidx) {
  igraph_real_t res;
  long int i, n, maxidx;

  assert(m != NULL);
  n=igraph_vector_size(&m->data);
  if (n == 0) return 0.0;
   
  maxidx = -1;
  for (i=0; i<n; i++)
    if (VECTOR(m->data)[i] != 0.0 &&
	(maxidx == -1 || VECTOR(m->data)[i] >= VECTOR(m->data)[maxidx]))
      maxidx = i;
	
  if (maxidx == -1) return 0.0;
	
  res=VECTOR(m->data)[maxidx];
  if (ridx != 0) *ridx = VECTOR(m->ridx)[maxidx];
  if (cidx != 0) {
    igraph_vector_binsearch(&m->cidx, maxidx, &i);
    while (VECTOR(m->cidx)[i+1] == VECTOR(m->cidx)[i]) i++;
    *cidx = (igraph_real_t)i;
  }
  return res;
}

/**
 * \function igraph_spmatrix_max
 *
 * Returns the maximum element of a matrix.
 * If the matrix is empty, zero is returned.
 *
 * \param m the matrix object.
 * \param ridx the row index of the maximum element if not \c NULL.
 * \param cidx the column index of the maximum element if not \c NULL.
 *
 * Time complexity: O(n), the number of nonzero elements in the matrix.
 */
igraph_real_t igraph_spmatrix_max(const igraph_spmatrix_t *m,
    igraph_real_t *ridx, igraph_real_t *cidx) {
  igraph_real_t res;
  long int i, j, k, maxidx;

  assert(m != NULL);
  i=igraph_vector_size(&m->data);
  if (i == 0) return 0.0;

  maxidx=(long)igraph_vector_which_max(&m->data);
  res=VECTOR(m->data)[maxidx];
  if (res>=0.0 || i == m->nrow * m->ncol) {
    if (ridx != 0) *ridx = VECTOR(m->ridx)[maxidx];
    if (cidx != 0) {
      igraph_vector_binsearch(&m->cidx, maxidx, &i);
      while (VECTOR(m->cidx)[i+1] == VECTOR(m->cidx)[i]) i++;
      *cidx = (igraph_real_t)i;
    }
    return res;
  }
  /* the maximal nonzero element is negative and there is at least a
   * single zero
   */
  res=0.0;
  if (cidx != 0 || ridx != 0) {
    for (i=0; i < m->ncol; i++) {
      if (VECTOR(m->cidx)[i+1] - VECTOR(m->cidx)[i] < m->nrow) {
        if (cidx != 0) *cidx = i;
        if (ridx != 0) {
          for (j=VECTOR(m->cidx)[i], k=0; j < VECTOR(m->cidx)[i+1]; j++, k++) {
            if (VECTOR(m->ridx)[j] != k) {
              *ridx = k;
              break;
            }
          }
        }
        break;
      }
    }
  }

  return res;
}

int igraph_i_spmatrix_get_col_nonzero_indices(const igraph_spmatrix_t *m,
  igraph_vector_t *res, long int col) {
  assert(m != NULL);
  long int i, n = VECTOR(m->cidx)[col+1]-VECTOR(m->cidx)[col];
  IGRAPH_CHECK(igraph_vector_resize(res, n));
  for (i=VECTOR(m->cidx)[col], n=0; i<VECTOR(m->cidx)[col+1]; i++, n++)
    if (VECTOR(m->data)[i] != 0.0) VECTOR(*res)[n] = VECTOR(m->ridx)[i];
  return 0;
}