/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2009-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA

   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 "config.h"

#include "cs/cs.h"

#include "igraph_sparsemat.h"
#include "igraph_error.h"
#include "igraph_interface.h"
#include "igraph_constructors.h"
#include "igraph_memory.h"
#include "igraph_vector_ptr.h"
#include "igraph_attributes.h"

#include <string.h>

/**
 * \section about_sparsemat About sparse matrices
 *
 * <para>
 * The <code>igraph_sparsemat_t</code> data type stores sparse matrices,
 * i.e. matrices in which the majority of the elements are zero.
 * </para>
 *
 * <para>The data type is essentially a wrapper to some of the
 * functions in the CXSparse library, by Tim Davis, see
 * http://faculty.cse.tamu.edu/davis/suitesparse.html
 * </para>
 *
 * <para>
 * Matrices can be stored in two formats: triplet and
 * column-compressed. The triplet format is intended for sparse matrix
 * initialization, as it is easy to add new (non-zero) elements to
 * it. Most of the computations are done on sparse matrices in
 * column-compressed format, after the user has converted the triplet
 * matrix to column-compressed, via \ref igraph_sparsemat_compress().
 * </para>
 *
 * <para>
 * Both formats are dynamic, in the sense that new elements can be
 * added to them, possibly resulting the allocation of more memory.
 * </para>
 *
 * <para>
 * Row and column indices follow the C convention and are zero-based.
 * </para>
 *
 * <para>
 * \example examples/simple/igraph_sparsemat.c
 * \example examples/simple/igraph_sparsemat2.c
 * \example examples/simple/igraph_sparsemat3.c
 * \example examples/simple/igraph_sparsemat4.c
 * \example examples/simple/igraph_sparsemat5.c
 * \example examples/simple/igraph_sparsemat6.c
 * \example examples/simple/igraph_sparsemat7.c
 * \example examples/simple/igraph_sparsemat8.c
 * </para>
 */

/**
 * \function igraph_sparsemat_init
 * Initialize a sparse matrix, in triplet format
 *
 * This is the most common way to create a sparse matrix, together
 * with the \ref igraph_sparsemat_entry() function, which can be used to
 * add the non-zero elements one by one. Once done, the user can call
 * \ref igraph_sparsemat_compress() to convert the matrix to
 * column-compressed, to allow computations with it.
 *
 * </para><para>The user must call \ref igraph_sparsemat_destroy() on
 * the matrix to deallocate the memory, once the matrix is no more
 * needed.
 * \param A Pointer to a not yet initialized sparse matrix.
 * \param rows The number of rows in the matrix.
 * \param cols The number of columns.
 * \param nzmax The maximum number of non-zero elements in the
 *    matrix. It is not compulsory to get this right, but it is
 *    useful for the allocation of the proper amount of memory.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_init(igraph_sparsemat_t *A, int rows, int cols, int nzmax) {

    if (rows < 0) {
        IGRAPH_ERROR("Negative number of rows", IGRAPH_EINVAL);
    }
    if (cols < 0) {
        IGRAPH_ERROR("Negative number of columns", IGRAPH_EINVAL);
    }

    A->cs = cs_spalloc( rows, cols, nzmax, /*values=*/ 1,
                        /*triplet=*/ 1);
    if (!A->cs) {
        IGRAPH_ERROR("Cannot allocate memory for sparse matrix", IGRAPH_ENOMEM);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_copy
 * Copy a sparse matrix
 *
 * Create a sparse matrix object, by copying another one. The source
 * matrix can be either in triplet or column-compressed format.
 *
 * </para><para>
 * Exactly the same amount of memory will be allocated to the
 * copy matrix, as it is currently for the original one.
 * \param to Pointer to an uninitialized sparse matrix, the copy will
 *    be created here.
 * \param from The sparse matrix to copy.
 * \return Error code.
 *
 * Time complexity: O(n+nzmax), the number of columns plus the maximum
 * number of non-zero elements.
 */

int igraph_sparsemat_copy(igraph_sparsemat_t *to,
                          const igraph_sparsemat_t *from) {

    int ne = from->cs->nz == -1 ? from->cs->n + 1 : from->cs->nzmax;

    to->cs = cs_spalloc(from->cs->m, from->cs->n, from->cs->nzmax,
                        /*values=*/ 1,
                        /*triplet=*/ igraph_sparsemat_is_triplet(from));

    to->cs->nzmax = from->cs->nzmax;
    to->cs->m     = from->cs->m;
    to->cs->n     = from->cs->n;
    to->cs->nz    = from->cs->nz;

    memcpy(to->cs->p, from->cs->p, sizeof(int) * (size_t) ne);
    memcpy(to->cs->i, from->cs->i, sizeof(int) * (size_t) (from->cs->nzmax));
    memcpy(to->cs->x, from->cs->x, sizeof(double) * (size_t) (from->cs->nzmax));

    return 0;
}

/**
 * \function igraph_sparsemat_destroy
 * Deallocate memory used by a sparse matrix
 *
 * One destroyed, the sparse matrix must be initialized again, before
 * calling any other operation on it.
 * \param A The sparse matrix to destroy.
 *
 * Time complexity: O(1).
 */

void igraph_sparsemat_destroy(igraph_sparsemat_t *A) {
    cs_spfree(A->cs);
}

/**
 * \function igraph_sparsemat_realloc
 * Allocate more (or less) memory for a sparse matrix
 *
 * Sparse matrices automatically allocate more memory, as needed. To
 * control memory allocation, the user can call this function, to
 * allocate memory for a given number of non-zero elements.
 * \param A The sparse matrix, it can be in triplet or
 *    column-compressed format.
 * \param nzmax The new maximum number of non-zero elements.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_realloc(igraph_sparsemat_t *A, int nzmax) {
    return !cs_sprealloc(A->cs, nzmax);
}

/**
 * \function igraph_sparsemat_nrow
 * Number of rows
 *
 * \param A The input matrix, in triplet or column-compressed format.
 * \return The number of rows in the \p A matrix.
 *
 * Time complexity: O(1).
 */

long int igraph_sparsemat_nrow(const igraph_sparsemat_t *A) {
    return A->cs->m;
}

/**
 * \function igraph_sparsemat_ncol
 * Number of columns.
 *
 * \param A The input matrix, in triplet or column-compressed format.
 * \return The number of columns in the \p A matrix.
 *
 * Time complexity: O(1).
 */

long int igraph_sparsemat_ncol(const igraph_sparsemat_t *A) {
    return A->cs->n;
}

/**
 * \function igraph_sparsemat_type
 * Type of a sparse matrix (triplet or column-compressed)
 *
 * Gives whether a sparse matrix is stored in the triplet format or in
 * column-compressed format.
 * \param A The input matrix.
 * \return Either \c IGRAPH_SPARSEMAT_CC or \c
 * IGRAPH_SPARSEMAT_TRIPLET.
 *
 * Time complexity: O(1).
 */

igraph_sparsemat_type_t igraph_sparsemat_type(const igraph_sparsemat_t *A) {
    return A->cs->nz < 0 ? IGRAPH_SPARSEMAT_CC : IGRAPH_SPARSEMAT_TRIPLET;
}

/**
 * \function igraph_sparsemat_is_triplet
 * Is this sparse matrix in triplet format?
 *
 * Decides whether a sparse matrix is in triplet format.
 * \param A The input matrix.
 * \return One if the input matrix is in triplet format, zero
 * otherwise.
 *
 * Time complexity: O(1).
 */

igraph_bool_t igraph_sparsemat_is_triplet(const igraph_sparsemat_t *A) {
    return A->cs->nz >= 0;
}

/**
 * \function igraph_sparsemat_is_cc
 * Is this sparse matrix in column-compressed format?
 *
 * Decides whether a sparse matrix is in column-compressed format.
 * \param A The input matrix.
 * \return One if the input matrix is in column-compressed format, zero
 * otherwise.
 *
 * Time complexity: O(1).
 */

igraph_bool_t igraph_sparsemat_is_cc(const igraph_sparsemat_t *A) {
    return A->cs->nz < 0;
}

/**
 * \function igraph_sparsemat_permute
 * Permute the rows and columns of a sparse matrix
 *
 * \param A The input matrix, it must be in column-compressed format.
 * \param p Integer vector, giving the permutation of the rows.
 * \param q Integer vector, the permutation of the columns.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *   stored here.
 * \return Error code.
 *
 * Time complexity: O(m+n+nz), the number of rows plus the number of
 * columns plus the number of non-zero elements in the matrix.
 */

int igraph_sparsemat_permute(const igraph_sparsemat_t *A,
                             const igraph_vector_int_t *p,
                             const igraph_vector_int_t *q,
                             igraph_sparsemat_t *res) {

    long int nrow = A->cs->m, ncol = A->cs->n;
    igraph_vector_int_t pinv;
    long int i;

    if (nrow != igraph_vector_int_size(p)) {
        IGRAPH_ERROR("Invalid row permutation length", IGRAPH_FAILURE);
    }
    if (ncol != igraph_vector_int_size(q)) {
        IGRAPH_ERROR("Invalid column permutation length", IGRAPH_FAILURE);
    }

    /* We invert the permutation by hand */
    IGRAPH_CHECK(igraph_vector_int_init(&pinv, nrow));
    IGRAPH_FINALLY(igraph_vector_int_destroy, &pinv);
    for (i = 0; i < nrow; i++) {
        VECTOR(pinv)[ VECTOR(*p)[i] ] = (int) i;
    }

    /* And call the permutation routine */
    if (! (res->cs = cs_permute(A->cs, VECTOR(pinv), VECTOR(*q), /*values=*/ 1))) {
        IGRAPH_ERROR("Cannot index sparse matrix", IGRAPH_FAILURE);
    }

    igraph_vector_int_destroy(&pinv);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

int igraph_i_sparsemat_index_rows(const igraph_sparsemat_t *A,
                                  const igraph_vector_int_t *p,
                                  igraph_sparsemat_t *res,
                                  igraph_real_t *constres) {

    igraph_sparsemat_t II, II2;
    long int nrow = A->cs->m;
    long int idx_rows = igraph_vector_int_size(p);
    long int k;

    /* Create index matrix */
    IGRAPH_CHECK(igraph_sparsemat_init(&II2, (int) idx_rows, (int) nrow,
                                       (int) idx_rows));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &II2);
    for (k = 0; k < idx_rows; k++) {
        igraph_sparsemat_entry(&II2, (int) k, VECTOR(*p)[k], 1.0);
    }
    IGRAPH_CHECK(igraph_sparsemat_compress(&II2, &II));
    igraph_sparsemat_destroy(&II2);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &II);

    /* Multiply */
    IGRAPH_CHECK(igraph_sparsemat_multiply(&II, A, res));
    igraph_sparsemat_destroy(&II);
    IGRAPH_FINALLY_CLEAN(1);

    if (constres) {
        if (res->cs->p[1] != 0) {
            *constres = res->cs->x[0];
        } else {
            *constres = 0.0;
        }
    }

    return 0;
}

int igraph_i_sparsemat_index_cols(const igraph_sparsemat_t *A,
                                  const igraph_vector_int_t *q,
                                  igraph_sparsemat_t *res,
                                  igraph_real_t *constres) {

    igraph_sparsemat_t JJ, JJ2;
    long int ncol = A->cs->n;
    long int idx_cols = igraph_vector_int_size(q);
    long int k;

    /* Create index matrix */
    IGRAPH_CHECK(igraph_sparsemat_init(&JJ2, (int) ncol, (int) idx_cols,
                                       (int) idx_cols));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ2);
    for (k = 0; k < idx_cols; k++) {
        igraph_sparsemat_entry(&JJ2, VECTOR(*q)[k], (int) k, 1.0);
    }
    IGRAPH_CHECK(igraph_sparsemat_compress(&JJ2, &JJ));
    igraph_sparsemat_destroy(&JJ2);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ);

    /* Multiply */
    IGRAPH_CHECK(igraph_sparsemat_multiply(A, &JJ, res));
    igraph_sparsemat_destroy(&JJ);
    IGRAPH_FINALLY_CLEAN(1);

    if (constres) {
        if (res->cs->p [1] != 0) {
            *constres = res->cs->x [0];
        } else {
            *constres = 0.0;
        }
    }

    return 0;
}

/**
 * \function igraph_sparsemat_index
 * Index a sparse matrix, extract a submatrix, or a single element
 *
 * This function serves two purposes. First, it can extract
 * submatrices from a sparse matrix. Second, as a special case, it can
 * extract a single element from a sparse matrix.
 * \param A The input matrix, it must be in column-compressed format.
 * \param p An integer vector, or a null pointer. The selected row
 *    index or indices. A null pointer selects all rows.
 * \param q An integer vector, or a null pointer. The selected column
 *    index or indices. A null pointer selects all columns.
 * \param res Pointer to an uninitialized sparse matrix, or a null
 *    pointer. If not a null pointer, then the selected submatrix is
 *    stored here.
 * \param constres Pointer to a real variable or a null pointer. If
 *    not a null pointer, then the first non-zero element in the
 *    selected submatrix is stored here, if there is one. Otherwise
 *    zero is stored here. This behavior is handy if one
 *    wants to select a single entry from the matrix.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_index(const igraph_sparsemat_t *A,
                           const igraph_vector_int_t *p,
                           const igraph_vector_int_t *q,
                           igraph_sparsemat_t *res,
                           igraph_real_t *constres) {

    igraph_sparsemat_t II, JJ, II2, JJ2, tmp;
    long int nrow = A->cs->m;
    long int ncol = A->cs->n;
    long int idx_rows = p ? igraph_vector_int_size(p) : -1;
    long int idx_cols = q ? igraph_vector_int_size(q) : -1;
    long int k;

    igraph_sparsemat_t *myres = res, mres;

    if (!p && !q) {
        IGRAPH_ERROR("No index vectors", IGRAPH_EINVAL);
    }

    if (!res && (idx_rows != 1 || idx_cols != 1)) {
        IGRAPH_ERROR("Sparse matrix indexing: must give `res' if not a "
                     "single element is selected", IGRAPH_EINVAL);
    }

    if (!q) {
        return igraph_i_sparsemat_index_rows(A, p, res, constres);
    }
    if (!p) {
        return igraph_i_sparsemat_index_cols(A, q, res, constres);
    }

    if (!res) {
        myres = &mres;
    }

    /* Create first index matrix */
    IGRAPH_CHECK(igraph_sparsemat_init(&II2, (int) idx_rows, (int) nrow,
                                       (int) idx_rows));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &II2);
    for (k = 0; k < idx_rows; k++) {
        igraph_sparsemat_entry(&II2, (int) k, VECTOR(*p)[k], 1.0);
    }
    IGRAPH_CHECK(igraph_sparsemat_compress(&II2, &II));
    igraph_sparsemat_destroy(&II2);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &II);

    /* Create second index matrix */
    IGRAPH_CHECK(igraph_sparsemat_init(&JJ2, (int) ncol, (int) idx_cols,
                                       (int) idx_cols));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ2);
    for (k = 0; k < idx_cols; k++) {
        igraph_sparsemat_entry(&JJ2, VECTOR(*q)[k], (int) k, 1.0);
    }
    IGRAPH_CHECK(igraph_sparsemat_compress(&JJ2, &JJ));
    igraph_sparsemat_destroy(&JJ2);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &JJ);

    /* Multiply */
    IGRAPH_CHECK(igraph_sparsemat_multiply(&II, A, &tmp));
    igraph_sparsemat_destroy(&II);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
    IGRAPH_CHECK(igraph_sparsemat_multiply(&tmp, &JJ, myres));
    igraph_sparsemat_destroy(&tmp);
    igraph_sparsemat_destroy(&JJ);
    IGRAPH_FINALLY_CLEAN(2);

    if (constres) {
        if (myres->cs->p [1] != 0) {
            *constres = myres->cs->x [0];
        } else {
            *constres = 0.0;
        }
    }

    if (!res) {
        igraph_sparsemat_destroy(myres);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_entry
 * Add an element to a sparse matrix
 *
 * This function can be used to add the entries to a sparse matrix,
 * after initializing it with \ref igraph_sparsemat_init().
 * \param A The input matrix, it must be in triplet format.
 * \param row The row index of the entry to add.
 * \param col The column index of the entry to add.
 * \param elem The value of the entry.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_entry(igraph_sparsemat_t *A, int row, int col,
                           igraph_real_t elem) {

    if (!cs_entry(A->cs, row, col, elem)) {
        IGRAPH_ERROR("Cannot add entry to sparse matrix",
                     IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_compress
 * Compress a sparse matrix, i.e. convert it to column-compress format
 *
 * Almost all sparse matrix operations require that the matrix is in
 * column-compressed format.
 * \param A The input matrix, it must be in triplet format.
 * \param res Pointer to an uninitialized sparse matrix object, the
 *    compressed version of \p A is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_compress(const igraph_sparsemat_t *A,
                              igraph_sparsemat_t *res) {

    if (! (res->cs = cs_compress(A->cs)) ) {
        IGRAPH_ERROR("Cannot compress sparse matrix", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_transpose
 * Transpose a sparse matrix
 *
 * \param A The input matrix, column-compressed or triple format.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *    stored here.
 * \param values If this is non-zero, the matrix transpose is
 *    calculated the normal way. If it is zero, then only the pattern
 *    of the input matrix is stored in the result, the values are not.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_transpose(const igraph_sparsemat_t *A,
                               igraph_sparsemat_t *res,
                               int values) {

    if (A->cs->nz < 0) {
        /* column-compressed */
        if (! (res->cs = cs_transpose(A->cs, values)) ) {
            IGRAPH_ERROR("Cannot transpose sparse matrix", IGRAPH_FAILURE);
        }
    } else {
        /* triplets */
        int *tmp;
        IGRAPH_CHECK(igraph_sparsemat_copy(res, A));
        tmp = res->cs->p;
        res->cs->p = res->cs->i;
        res->cs->i = tmp;
    }
    return 0;
}

igraph_bool_t
igraph_i_sparsemat_is_symmetric_cc(const igraph_sparsemat_t *A) {
    igraph_sparsemat_t t, tt;
    igraph_bool_t res;
    int nz;

    IGRAPH_CHECK(igraph_sparsemat_transpose(A, &t, /*values=*/ 1));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &t);
    IGRAPH_CHECK(igraph_sparsemat_dupl(&t));
    IGRAPH_CHECK(igraph_sparsemat_transpose(&t, &tt, /*values=*/ 1));
    igraph_sparsemat_destroy(&t);
    IGRAPH_FINALLY_CLEAN(1);
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &tt);
    IGRAPH_CHECK(igraph_sparsemat_transpose(&tt, &t, /*values=*/ 1));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &t);

    nz = t.cs->p[t.cs->n];
    res = memcmp(t.cs->i, tt.cs->i, sizeof(int) * (size_t) nz) == 0;
    res = res && memcmp(t.cs->p, tt.cs->p, sizeof(int) *
                        (size_t)(t.cs->n + 1)) == 0;
    res = res && memcmp(t.cs->x, tt.cs->x, sizeof(igraph_real_t) * (size_t)nz) == 0;

    igraph_sparsemat_destroy(&t);
    igraph_sparsemat_destroy(&tt);
    IGRAPH_FINALLY_CLEAN(2);

    return res;
}

igraph_bool_t
igraph_i_sparsemat_is_symmetric_triplet(const igraph_sparsemat_t *A) {
    igraph_sparsemat_t tmp;
    igraph_bool_t res;
    IGRAPH_CHECK(igraph_sparsemat_compress(A, &tmp));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
    res = igraph_i_sparsemat_is_symmetric_cc(&tmp);
    igraph_sparsemat_destroy(&tmp);
    IGRAPH_FINALLY_CLEAN(1);
    return res;
}

igraph_bool_t igraph_sparsemat_is_symmetric(const igraph_sparsemat_t *A) {

    if (A->cs->m != A->cs->n) {
        return 0;
    }

    if (A->cs->nz < 0) {
        return igraph_i_sparsemat_is_symmetric_cc(A);
    } else {
        return igraph_i_sparsemat_is_symmetric_triplet(A);
    }
}

/**
 * \function igraph_sparsemat_dupl
 * Remove duplicate elements from a sparse matrix
 *
 * It is possible that a column-compressed sparse matrix stores a
 * single matrix entry in multiple pieces. The entry is then the sum
 * of all its pieces. (Some functions create matrices like this.) This
 * function eliminates the multiple pieces.
 * \param A The input matrix, in column-compressed format.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_dupl(igraph_sparsemat_t *A) {

    if (!cs_dupl(A->cs)) {
        IGRAPH_ERROR("Cannot remove duplicates from sparse matrix",
                     IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_fkeep
 * Filter the elements of a sparse matrix
 *
 * This function can be used to filter the (non-zero) elements of a
 * sparse matrix. For all entries, it calls the supplied function and
 * depending on the return values either keeps, or deleted the element
 * from the matrix.
 * \param A The input matrix, in column-compressed format.
 * \param fkeep The filter function. It must take four arguments: the
 *    first is an \c int, the row index of the entry, the second is
 *    another \c int, the column index. The third is \c igraph_real_t,
 *    the value of the entry. The fourth element is a \c void pointer,
 *    the \p other argument is passed here. The function must return
 *    an \c int. If this is zero, then the entry is deleted, otherwise
 *    it is kept.
 * \param other A \c void pointer that is passed to the filtering
 * function.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_fkeep(igraph_sparsemat_t *A,
                           int (*fkeep)(int, int, igraph_real_t, void*),
                           void *other) {

    if (!cs_fkeep(A->cs, fkeep, other)) {
        IGRAPH_ERROR("Cannot filter sparse matrix", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_dropzeros
 * Drop the zero elements from a sparse matrix
 *
 * As a result of matrix operations, some of the entries in a sparse
 * matrix might be zero. This function removes these entries.
 * \param A The input matrix, it must be in column-compressed format.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_dropzeros(igraph_sparsemat_t *A) {

    if (!cs_dropzeros(A->cs)) {
        IGRAPH_ERROR("Cannot drop zeros from sparse matrix", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_droptol
 * Drop the almost zero elements of a sparse matrix
 *
 * This function is similar to \ref igraph_sparsemat_dropzeros(), but it
 * also drops entries that are closer to zero than the given tolerance
 * threshold.
 * \param A The input matrix, it must be in column-compressed format.
 * \param tol Real number, giving the tolerance threshold.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_droptol(igraph_sparsemat_t *A, igraph_real_t tol) {

    if (!cs_droptol(A->cs, tol)) {
        IGRAPH_ERROR("Cannot drop (almost) zeros from sparse matrix",
                     IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_multiply
 * Matrix multiplication
 *
 * Multiplies two sparse matrices.
 * \param A The first input matrix (left hand side), in
 *   column-compressed format.
 * \param B The second input matrix (right hand side), in
 *   column-compressed format.
 * \param res Pointer to an uninitialized sparse matrix, the result is
 *   stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_multiply(const igraph_sparsemat_t *A,
                              const igraph_sparsemat_t *B,
                              igraph_sparsemat_t *res) {

    if (! (res->cs = cs_multiply(A->cs, B->cs))) {
        IGRAPH_ERROR("Cannot multiply matrices", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_add
 * Sum of two sparse matrices
 *
 * \param A The first input matrix, in column-compressed format.
 * \param B The second input matrix, in column-compressed format.
 * \param alpha Real scalar, \p A is multiplied by \p alpha before the
 *    addition.
 * \param beta Real scalar, \p B is multiplied by \p beta before the
 *    addition.
 * \param res Pointer to an uninitialized sparse matrix, the result
 *    is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_add(const igraph_sparsemat_t *A,
                         const igraph_sparsemat_t *B,
                         igraph_real_t alpha,
                         igraph_real_t beta,
                         igraph_sparsemat_t *res) {

    if (! (res->cs = cs_add(A->cs, B->cs, alpha, beta))) {
        IGRAPH_ERROR("Cannot add matrices", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_gaxpy
 * Matrix-vector product, added to another vector.
 *
 * \param A The input matrix, in column-compressed format.
 * \param x The input vector, its size must match the number of
 *    columns in \p A.
 * \param res This vector is added to the matrix-vector product
 *    and it is overwritten by the result.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_gaxpy(const igraph_sparsemat_t *A,
                           const igraph_vector_t *x,
                           igraph_vector_t *res) {

    if (A->cs->n != igraph_vector_size(x) ||
        A->cs->m != igraph_vector_size(res)) {
        IGRAPH_ERROR("Invalid matrix/vector size for multiplication",
                     IGRAPH_EINVAL);
    }

    if (! (cs_gaxpy(A->cs, VECTOR(*x), VECTOR(*res)))) {
        IGRAPH_ERROR("Cannot perform sparse matrix vector multiplication",
                     IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_lsolve
 * Solve a lower-triangular linear system
 *
 * Solve the Lx=b linear equation system, where the L coefficient
 * matrix is square and lower-triangular, with a zero-free diagonal.
 * \param L The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_lsolve(const igraph_sparsemat_t *L,
                            const igraph_vector_t *b,
                            igraph_vector_t *res) {

    if (L->cs->m != L->cs->n) {
        IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (! cs_lsolve(L->cs, VECTOR(*res))) {
        IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_ltsolve
 * Solve an upper-triangular linear system
 *
 * Solve the L'x=b linear equation system, where the L
 * matrix is square and lower-triangular, with a zero-free diagonal.
 * \param L The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_ltsolve(const igraph_sparsemat_t *L,
                             const igraph_vector_t *b,
                             igraph_vector_t *res) {

    if (L->cs->m != L->cs->n) {
        IGRAPH_ERROR("Cannot perform transposed lower triangular solve",
                     IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (!cs_ltsolve(L->cs, VECTOR(*res))) {
        IGRAPH_ERROR("Cannot perform lower triangular solve", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_usolve
 * Solve an upper-triangular linear system
 *
 * Solves the Ux=b upper triangular system.
 * \param U The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_usolve(const igraph_sparsemat_t *U,
                            const igraph_vector_t *b,
                            igraph_vector_t *res) {

    if (U->cs->m != U->cs->n) {
        IGRAPH_ERROR("Cannot perform upper triangular solve", IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (! cs_usolve(U->cs, VECTOR(*res))) {
        IGRAPH_ERROR("Cannot perform upper triangular solve", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_utsolve
 * Solve a lower-triangular linear system
 *
 * This is the same as \ref igraph_sparsemat_usolve(), but U'x=b is
 * solved, where the apostrophe denotes the transpose.
 * \param U The input matrix, in column-compressed format.
 * \param b The right hand side of the linear system.
 * \param res An initialized vector, the result is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_utsolve(const igraph_sparsemat_t *U,
                             const igraph_vector_t *b,
                             igraph_vector_t *res) {

    if (U->cs->m != U->cs->n) {
        IGRAPH_ERROR("Cannot perform transposed upper triangular solve",
                     IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (!cs_utsolve(U->cs, VECTOR(*res))) {
        IGRAPH_ERROR("Cannot perform transposed upper triangular solve",
                     IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_cholsol
 * Solve a symmetric linear system via Cholesky decomposition
 *
 * Solve Ax=b, where A is a symmetric positive definite matrix.
 * \param A The input matrix, in column-compressed format.
 * \param v The right hand side.
 * \param res An initialized vector, the result is stored here.
 * \param order An integer giving the ordering method to use for the
 *    factorization. Zero is the natural ordering; if it is one, then
 *    the fill-reducing minimum-degree ordering of A+A' is used.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_cholsol(const igraph_sparsemat_t *A,
                             const igraph_vector_t *b,
                             igraph_vector_t *res,
                             int order) {

    if (A->cs->m != A->cs->n) {
        IGRAPH_ERROR("Cannot perform sparse symmetric solve",
                     IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (! cs_cholsol(order, A->cs, VECTOR(*res))) {
        IGRAPH_ERROR("Cannot perform sparse symmetric solve", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_lusol
 * Solve a linear system via LU decomposition
 *
 * Solve Ax=b, via LU factorization of A.
 * \param A The input matrix, in column-compressed format.
 * \param b The right hand side of the equation.
 * \param res An initialized vector, the result is stored here.
 * \param order The ordering method to use, zero means the natural
 *    ordering, one means the fill-reducing minimum-degree ordering of
 *    A+A', two means the ordering of A'*A, after removing the dense
 *    rows from A. Three means the ordering of A'*A.
 * \param tol Real number, the tolerance limit to use for the numeric
 *    LU factorization.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_lusol(const igraph_sparsemat_t *A,
                           const igraph_vector_t *b,
                           igraph_vector_t *res,
                           int order,
                           igraph_real_t tol) {

    if (A->cs->m != A->cs->n) {
        IGRAPH_ERROR("Cannot perform LU solve",
                     IGRAPH_NONSQUARE);
    }

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    if (! cs_lusol(order, A->cs, VECTOR(*res), tol)) {
        IGRAPH_ERROR("Cannot perform LU solve", IGRAPH_FAILURE);
    }

    return 0;
}

int igraph_i_sparsemat_cc(igraph_t *graph, const igraph_sparsemat_t *A,
                          igraph_bool_t directed) {

    igraph_vector_t edges;
    long int no_of_nodes = A->cs->m;
    long int no_of_edges = A->cs->p[A->cs->n];
    int *p = A->cs->p;
    int *i = A->cs->i;
    long int from = 0;
    long int to = 0;
    long int e = 0;

    if (no_of_nodes != A->cs->n) {
        IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
    }

    IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges * 2);

    while (*p < no_of_edges) {
        while (to < * (p + 1)) {
            if (directed || from >= *i) {
                VECTOR(edges)[e++] = from;
                VECTOR(edges)[e++] = (*i);
            }
            to++;
            i++;
        }
        from++;
        p++;
    }
    igraph_vector_resize(&edges, e);

    IGRAPH_CHECK(igraph_create(graph, &edges, (igraph_integer_t) no_of_nodes,
                               directed));
    igraph_vector_destroy(&edges);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

int igraph_i_sparsemat_triplet(igraph_t *graph, const igraph_sparsemat_t *A,
                               igraph_bool_t directed) {

    igraph_vector_t edges;
    long int no_of_nodes = A->cs->m;
    long int no_of_edges = A->cs->nz;
    int *i = A->cs->p;
    int *j = A->cs->i;
    long int e;

    if (no_of_nodes != A->cs->n) {
        IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
    }

    IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges * 2);

    for (e = 0; e < 2 * no_of_edges; i++, j++) {
        if (directed || *i >= *j) {
            VECTOR(edges)[e++] = (*i);
            VECTOR(edges)[e++] = (*j);
        }
    }
    igraph_vector_resize(&edges, e);

    IGRAPH_CHECK(igraph_create(graph, &edges, (igraph_integer_t) no_of_nodes,
                               directed));
    igraph_vector_destroy(&edges);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

/**
 * \function igraph_sparsemat
 * Create an igraph graph from a sparse matrix
 *
 * One edge is created for each non-zero entry in the matrix. If you
 * have a symmetric matrix, and want to create an undirected graph,
 * then delete the entries in the upper diagonal first, or call \ref
 * igraph_simplify() on the result graph to eliminate the multiple
 * edges.
 * \param graph Pointer to an uninitialized igraph_t object, the
 *    graphs is stored here.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param directed Boolean scalar, whether to create a directed
 *    graph.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat(igraph_t *graph, const igraph_sparsemat_t *A,
                     igraph_bool_t directed) {

    if (A->cs->nz < 0) {
        return (igraph_i_sparsemat_cc(graph, A, directed));
    } else {
        return (igraph_i_sparsemat_triplet(graph, A, directed));
    }
}

int igraph_i_weighted_sparsemat_cc(const igraph_sparsemat_t *A,
                                   igraph_bool_t directed, const char *attr,
                                   igraph_bool_t loops,
                                   igraph_vector_t *edges,
                                   igraph_vector_t *weights) {

    long int no_of_edges = A->cs->p[A->cs->n];
    int *p = A->cs->p;
    int *i = A->cs->i;
    igraph_real_t *x = A->cs->x;
    long int from = 0;
    long int to = 0;
    long int e = 0, w = 0;

    IGRAPH_UNUSED(attr);

    igraph_vector_resize(edges, no_of_edges * 2);
    igraph_vector_resize(weights, no_of_edges);

    while (*p < no_of_edges) {
        while (to < * (p + 1)) {
            if ( (loops || from != *i) && (directed || from >= *i) && *x != 0) {
                VECTOR(*edges)[e++] = (*i);
                VECTOR(*edges)[e++] = from;
                VECTOR(*weights)[w++] = (*x);
            }
            to++;
            i++;
            x++;
        }
        from++;
        p++;
    }

    igraph_vector_resize(edges, e);
    igraph_vector_resize(weights, w);

    return 0;
}

int igraph_i_weighted_sparsemat_triplet(const igraph_sparsemat_t *A,
                                        igraph_bool_t directed,
                                        const char *attr,
                                        igraph_bool_t loops,
                                        igraph_vector_t *edges,
                                        igraph_vector_t *weights) {

    IGRAPH_UNUSED(A); IGRAPH_UNUSED(directed); IGRAPH_UNUSED(attr);
    IGRAPH_UNUSED(loops); IGRAPH_UNUSED(edges); IGRAPH_UNUSED(weights);

    /* TODO */
    IGRAPH_ERROR("Triplet matrices are not implemented",
                 IGRAPH_UNIMPLEMENTED);
    return 0;
}

int igraph_weighted_sparsemat(igraph_t *graph, const igraph_sparsemat_t *A,
                              igraph_bool_t directed, const char *attr,
                              igraph_bool_t loops) {

    igraph_vector_t edges, weights;
    int pot_edges = A->cs->nz < 0 ? A->cs->p[A->cs->n] : A->cs->nz;
    const char* default_attr = "weight";
    igraph_vector_ptr_t attr_vec;
    igraph_attribute_record_t attr_rec;
    long int no_of_nodes = A->cs->m;

    if (no_of_nodes != A->cs->n) {
        IGRAPH_ERROR("Cannot create graph object", IGRAPH_NONSQUARE);
    }

    IGRAPH_VECTOR_INIT_FINALLY(&edges, pot_edges * 2);
    IGRAPH_VECTOR_INIT_FINALLY(&weights, pot_edges);
    IGRAPH_VECTOR_PTR_INIT_FINALLY(&attr_vec, 1);

    if (A->cs->nz < 0) {
        IGRAPH_CHECK(igraph_i_weighted_sparsemat_cc(A, directed, attr, loops,
                     &edges, &weights));
    } else {
        IGRAPH_CHECK(igraph_i_weighted_sparsemat_triplet(A, directed, attr,
                     loops, &edges,
                     &weights));
    }

    /* Prepare attribute record */
    attr_rec.name = attr ? attr : default_attr;
    attr_rec.type = IGRAPH_ATTRIBUTE_NUMERIC;
    attr_rec.value = &weights;
    VECTOR(attr_vec)[0] = &attr_rec;

    /* Create graph */
    IGRAPH_CHECK(igraph_empty(graph, (igraph_integer_t) no_of_nodes, directed));
    IGRAPH_FINALLY(igraph_destroy, graph);
    if (igraph_vector_size(&edges) > 0) {
        IGRAPH_CHECK(igraph_add_edges(graph, &edges, &attr_vec));
    }
    IGRAPH_FINALLY_CLEAN(1);

    /* Cleanup */
    igraph_vector_destroy(&edges);
    igraph_vector_destroy(&weights);
    igraph_vector_ptr_destroy(&attr_vec);
    IGRAPH_FINALLY_CLEAN(3);

    return 0;
}

/**
 * \function igraph_get_sparsemat
 * Convert an igraph graph to a sparse matrix
 *
 * If the graph is undirected, then a symmetric matrix is created.
 * \param graph The input graph.
 * \param res Pointer to an uninitialized sparse matrix. The result
 *    will be stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_get_sparsemat(const igraph_t *graph, igraph_sparsemat_t *res) {

    long int no_of_nodes = igraph_vcount(graph);
    long int no_of_edges = igraph_ecount(graph);
    igraph_bool_t directed = igraph_is_directed(graph);
    long int nzmax = directed ? no_of_edges : no_of_edges * 2;
    long int i;

    IGRAPH_CHECK(igraph_sparsemat_init(res, (igraph_integer_t) no_of_nodes,
                                       (igraph_integer_t) no_of_nodes,
                                       (igraph_integer_t) nzmax));

    for (i = 0; i < no_of_edges; i++) {
        long int from = IGRAPH_FROM(graph, i);
        long int to = IGRAPH_TO(graph, i);
        IGRAPH_CHECK(igraph_sparsemat_entry(res, (int) from, (int) to, 1.0));
        if (!directed && from != to) {
            IGRAPH_CHECK(igraph_sparsemat_entry(res, (int) to, (int) from, 1.0));
        }
    }

    return 0;
}

#define CHECK(x) if ((x)<0) { IGRAPH_ERROR("Cannot write to file", IGRAPH_EFILE); }

/**
 * \function igraph_sparsemat_print
 * Print a sparse matrix to a file
 *
 * Only the non-zero entries are printed. This function serves more as
 * a debugging utility, as currently there is no function that could
 * read back the printed matrix from the file.
 * \param A The input matrix, triplet or column-compressed format.
 * \param outstream The stream to print it to.
 * \return Error code.
 *
 * Time complexity: O(nz) for triplet matrices, O(n+nz) for
 * column-compressed matrices. nz is the number of non-zero elements,
 * n is the number columns in the matrix.
 */

int igraph_sparsemat_print(const igraph_sparsemat_t *A,
                           FILE *outstream) {

    if (A->cs->nz < 0) {
        /* CC */
        int j, p;
        for (j = 0; j < A->cs->n; j++) {
            CHECK(fprintf(outstream, "col %i: locations %i to %i\n",
                          j, A->cs->p[j], A->cs->p[j + 1] - 1));
            for (p = A->cs->p[j]; p < A->cs->p[j + 1]; p++) {
                CHECK(fprintf(outstream, "%i : %g\n", A->cs->i[p], A->cs->x[p]));
            }
        }
    } else {
        /* Triplet */
        int p;
        for (p = 0; p < A->cs->nz; p++) {
            CHECK(fprintf(outstream, "%i %i : %g\n",
                          A->cs->i[p], A->cs->p[p], A->cs->x[p]));
        }
    }

    return 0;
}

#undef CHECK

int igraph_i_sparsemat_eye_triplet(igraph_sparsemat_t *A, int n, int nzmax,
                                   igraph_real_t value) {
    long int i;

    IGRAPH_CHECK(igraph_sparsemat_init(A, n, n, nzmax));

    for (i = 0; i < n; i++) {
        igraph_sparsemat_entry(A, (int) i, (int) i, value);
    }

    return 0;
}

int igraph_i_sparsemat_eye_cc(igraph_sparsemat_t *A, int n,
                              igraph_real_t value) {
    long int i;

    if (! (A->cs = cs_spalloc(n, n, n, /*values=*/ 1, /*triplet=*/ 0)) ) {
        IGRAPH_ERROR("Cannot create eye sparse matrix", IGRAPH_FAILURE);
    }

    for (i = 0; i < n; i++) {
        A->cs->p [i] = (int) i;
        A->cs->i [i] = (int) i;
        A->cs->x [i] = value;
    }
    A->cs->p [n] = n;

    return 0;
}

/**
 * \function igraph_sparsemat_eye
 * Create a sparse identity matrix
 *
 * \param A An uninitialized sparse matrix, the result is stored
 *   here.
 * \param n The number of rows and number of columns in the matrix.
 * \param nzmax The maximum number of non-zero elements, this
 *   essentially gives the amount of memory that will be allocated for
 *   matrix elements.
 * \param value The value to store in the diagonal.
 * \param compress Whether to create a column-compressed matrix. If
 *   false, then a triplet matrix is created.
 * \return Error code.
 *
 * Time complexity: O(n).
 */

int igraph_sparsemat_eye(igraph_sparsemat_t *A, int n, int nzmax,
                         igraph_real_t value,
                         igraph_bool_t compress) {
    if (compress) {
        return (igraph_i_sparsemat_eye_cc(A, n, value));
    } else {
        return (igraph_i_sparsemat_eye_triplet(A, n, nzmax, value));
    }
}

int igraph_i_sparsemat_diag_triplet(igraph_sparsemat_t *A, int nzmax,
                                    const igraph_vector_t *values) {

    int i, n = (int) igraph_vector_size(values);

    IGRAPH_CHECK(igraph_sparsemat_init(A, n, n, nzmax));

    for (i = 0; i < n; i++) {
        igraph_sparsemat_entry(A, i, i, VECTOR(*values)[i]);
    }

    return 0;

}

int igraph_i_sparsemat_diag_cc(igraph_sparsemat_t *A,
                               const igraph_vector_t *values) {

    int i, n = (int) igraph_vector_size(values);

    if (! (A->cs = cs_spalloc(n, n, n, /*values=*/ 1, /*triplet=*/ 0)) ) {
        IGRAPH_ERROR("Cannot create eye sparse matrix", IGRAPH_FAILURE);
    }

    for (i = 0; i < n; i++) {
        A->cs->p [i] = i;
        A->cs->i [i] = i;
        A->cs->x [i] = VECTOR(*values)[i];
    }
    A->cs->p [n] = n;

    return 0;

}

/**
 * \function igraph_sparsemat_diag
 * Create a sparse diagonal matrix
 *
 * \param A An uninitialized sparse matrix, the result is stored
 *    here.
 * \param nzmax The maximum number of non-zero elements, this
 *   essentially gives the amount of memory that will be allocated for
 *   matrix elements.
 * \param values The values to store in the diagonal, the size of the
 *    matrix defined by the length of this vector.
 * \param compress Whether to create a column-compressed matrix. If
 *   false, then a triplet matrix is created.
 * \return Error code.
 *
 * Time complexity: O(n), the length of the diagonal vector.
 */

int igraph_sparsemat_diag(igraph_sparsemat_t *A, int nzmax,
                          const igraph_vector_t *values,
                          igraph_bool_t compress) {

    if (compress) {
        return (igraph_i_sparsemat_diag_cc(A, values));
    } else {
        return (igraph_i_sparsemat_diag_triplet(A, nzmax, values));
    }
}

int igraph_i_sparsemat_arpack_multiply(igraph_real_t *to,
                                       const igraph_real_t *from,
                                       int n,
                                       void *extra) {
    igraph_sparsemat_t *A = extra;
    igraph_vector_t vto, vfrom;
    igraph_vector_view(&vto, to, n);
    igraph_vector_view(&vfrom, from, n);
    igraph_vector_null(&vto);
    IGRAPH_CHECK(igraph_sparsemat_gaxpy(A, &vfrom, &vto));
    return 0;
}

typedef struct igraph_i_sparsemat_arpack_rssolve_data_t {
    igraph_sparsemat_symbolic_t *dis;
    igraph_sparsemat_numeric_t *din;
    igraph_real_t tol;
    igraph_sparsemat_solve_t method;
} igraph_i_sparsemat_arpack_rssolve_data_t;

int igraph_i_sparsemat_arpack_solve(igraph_real_t *to,
                                    const igraph_real_t *from,
                                    int n,
                                    void *extra) {

    igraph_i_sparsemat_arpack_rssolve_data_t *data = extra;
    igraph_vector_t vfrom, vto;

    igraph_vector_view(&vfrom, from, n);
    igraph_vector_view(&vto, to, n);

    if (data->method == IGRAPH_SPARSEMAT_SOLVE_LU) {
        IGRAPH_CHECK(igraph_sparsemat_luresol(data->dis, data->din, &vfrom,
                                              &vto));
    } else if (data->method == IGRAPH_SPARSEMAT_SOLVE_QR) {
        IGRAPH_CHECK(igraph_sparsemat_qrresol(data->dis, data->din, &vfrom,
                                              &vto));

    }

    return 0;
}

/**
 * \function igraph_sparsemat_arpack_rssolve
 * Eigenvalues and eigenvectors of a symmetric sparse matrix via ARPACK
 *
 * \param The input matrix, must be column-compressed.
 * \param options It is passed to \ref igraph_arpack_rssolve(). See
 *    \ref igraph_arpack_options_t for the details. If \c mode is 1,
 *    then ARPACK uses regular mode, if \c mode is 3, then shift and
 *    invert mode is used and the \c sigma structure member defines
 *    the shift.
 * \param storage Storage for ARPACK. See \ref
 *    igraph_arpack_rssolve() and \ref igraph_arpack_storage_t for
 *    details.
 * \param values An initialized vector or a null pointer, the
 *    eigenvalues are stored here.
 * \param vectors An initialised matrix, or a null pointer, the
 *    eigenvectors are stored here, in the columns.
 * \param solvemethod The method to solve the linear system, if \c
 *    mode is 3, i.e. the shift and invert mode is used.
 *    Possible values:
 *    \clist
 *      \cli IGRAPH_SPARSEMAT_SOLVE_LU
 *           The linear system is solved using LU decomposition.
 *      \cli IGRAPH_SPARSEMAT_SOLVE_QR
 *           The linear system is solved using QR decomposition.
 *    \endclist
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_arpack_rssolve(const igraph_sparsemat_t *A,
                                    igraph_arpack_options_t *options,
                                    igraph_arpack_storage_t *storage,
                                    igraph_vector_t *values,
                                    igraph_matrix_t *vectors,
                                    igraph_sparsemat_solve_t solvemethod) {

    int n = (int) igraph_sparsemat_nrow(A);

    if (n != igraph_sparsemat_ncol(A)) {
        IGRAPH_ERROR("Non-square matrix for ARPACK", IGRAPH_NONSQUARE);
    }

    options->n = n;

    if (options->mode == 1) {
        IGRAPH_CHECK(igraph_arpack_rssolve(igraph_i_sparsemat_arpack_multiply,
                                           (void*) A, options, storage,
                                           values, vectors));
    } else if (options->mode == 3) {
        igraph_real_t sigma = options->sigma;
        igraph_sparsemat_t OP, eye;
        igraph_sparsemat_symbolic_t symb;
        igraph_sparsemat_numeric_t num;
        igraph_i_sparsemat_arpack_rssolve_data_t data;
        /*-----------------------------------*/
        /* We need to factor the (A-sigma*I) */
        /*-----------------------------------*/

        /* Create (A-sigma*I) */
        IGRAPH_CHECK(igraph_sparsemat_eye(&eye, /*n=*/ n, /*nzmax=*/ n,
                                          /*value=*/ -sigma, /*compress=*/ 1));
        IGRAPH_FINALLY(igraph_sparsemat_destroy, &eye);
        IGRAPH_CHECK(igraph_sparsemat_add(/*A=*/ A, /*B=*/ &eye, /*alpha=*/ 1.0,
                     /*beta=*/ 1.0, /*res=*/ &OP));
        igraph_sparsemat_destroy(&eye);
        IGRAPH_FINALLY_CLEAN(1);
        IGRAPH_FINALLY(igraph_sparsemat_destroy, &OP);

        if (solvemethod == IGRAPH_SPARSEMAT_SOLVE_LU) {
            /* Symbolic analysis */
            IGRAPH_CHECK(igraph_sparsemat_symblu(/*order=*/ 0, &OP, &symb));
            IGRAPH_FINALLY(igraph_sparsemat_symbolic_destroy, &symb);
            /* Numeric LU factorization */
            IGRAPH_CHECK(igraph_sparsemat_lu(&OP, &symb, &num, /*tol=*/ 0));
            IGRAPH_FINALLY(igraph_sparsemat_numeric_destroy, &num);
        } else if (solvemethod == IGRAPH_SPARSEMAT_SOLVE_QR) {
            /* Symbolic analysis */
            IGRAPH_CHECK(igraph_sparsemat_symbqr(/*order=*/ 0, &OP, &symb));
            IGRAPH_FINALLY(igraph_sparsemat_symbolic_destroy, &symb);
            /* Numeric QR factorization */
            IGRAPH_CHECK(igraph_sparsemat_qr(&OP, &symb, &num));
            IGRAPH_FINALLY(igraph_sparsemat_numeric_destroy, &num);
        }

        data.dis = &symb;
        data.din = &num;
        data.tol = options->tol;
        data.method = solvemethod;
        IGRAPH_CHECK(igraph_arpack_rssolve(igraph_i_sparsemat_arpack_solve,
                                           (void*) &data, options, storage,
                                           values, vectors));

        igraph_sparsemat_numeric_destroy(&num);
        igraph_sparsemat_symbolic_destroy(&symb);
        igraph_sparsemat_destroy(&OP);
        IGRAPH_FINALLY_CLEAN(3);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_arpack_rnsolve
 * Eigenvalues and eigenvectors of a nonsymmetric sparse matrix via ARPACK
 *
 * Eigenvalues and/or eigenvectors of a nonsymmetric sparse matrix.
 * \param A The input matrix, in column-compressed mode.
 * \param options ARPACK options, it is passed to \ref
 *    igraph_arpack_rnsolve(). See also \ref igraph_arpack_options_t
 *    for details.
 * \param storage Storage for ARPACK, this is passed to \ref
 *    igraph_arpack_rnsolve(). See \ref igraph_arpack_storage_t for
 *    details.
 * \param values An initialized matrix, or a null pointer. If not a
 *    null pointer, then the eigenvalues are stored here, the first
 *    column is the real part, the second column is the imaginary
 *    part.
 * \param vectors An initialized matrix, or a null pointer. If not a
 *    null pointer, then the eigenvectors are stored here, please see
 *    \ref igraph_arpack_rnsolve() for the format.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_arpack_rnsolve(const igraph_sparsemat_t *A,
                                    igraph_arpack_options_t *options,
                                    igraph_arpack_storage_t *storage,
                                    igraph_matrix_t *values,
                                    igraph_matrix_t *vectors) {

    int n = (int) igraph_sparsemat_nrow(A);

    if (n != igraph_sparsemat_ncol(A)) {
        IGRAPH_ERROR("Non-square matrix for ARPACK", IGRAPH_NONSQUARE);
    }

    options->n = n;

    return igraph_arpack_rnsolve(igraph_i_sparsemat_arpack_multiply,
                                 (void*) A, options, storage,
                                 values, vectors);
}

/**
 * \function igraph_sparsemat_symbqr
 * Symbolic QR decomposition
 *
 * QR decomposition of sparse matrices involves two steps, the first
 * is calling this function, and then \ref
 * igraph_sparsemat_qr().
 * \param order The ordering to use: 0 means natural ordering, 1 means
 *   minimum degree ordering of A+A', 2 is minimum degree ordering of
 *   A'A after removing the dense rows from A, and 3 is the minimum
 *   degree ordering of A'A.
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic analysis is stored here. Once
 *    not needed anymore, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy().
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_symbqr(long int order, const igraph_sparsemat_t *A,
                            igraph_sparsemat_symbolic_t *dis) {

    dis->symbolic = cs_sqr((int) order, A->cs, /*qr=*/ 1);
    if (!dis->symbolic) {
        IGRAPH_ERROR("Cannot do symbolic QR decomposition", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_symblu
 * Symbolic LU decomposition
 *
 * LU decomposition of sparse matrices involves two steps, the first
 * is calling this function, and then \ref igraph_sparsemat_lu().
 * \param order The ordering to use: 0 means natural ordering, 1 means
 *   minimum degree ordering of A+A', 2 is minimum degree ordering of
 *   A'A after removing the dense rows from A, and 3 is the minimum
 *   degree ordering of A'A.
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic analysis is stored here. Once
 *    not needed anymore, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy().
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_symblu(long int order, const igraph_sparsemat_t *A,
                            igraph_sparsemat_symbolic_t *dis) {

    dis->symbolic = cs_sqr((int) order, A->cs, /*qr=*/ 0);
    if (!dis->symbolic) {
        IGRAPH_ERROR("Cannot do symbolic LU decomposition", IGRAPH_FAILURE);
    }

    return 0;
}

/**
 * \function igraph_sparsemat_lu
 * LU decomposition of a sparse matrix
 *
 * Performs numeric sparse LU decomposition of a matrix.
 * \param A The input matrix, in column-compressed format.
 * \param dis The symbolic analysis for LU decomposition, coming from
 *    a call to the \ref igraph_sparsemat_symblu() function.
 * \param din The numeric decomposition, the result is stored here. It
 *    can be used to solve linear systems with changing right hand
 *    side vectors, by calling \ref igraph_sparsemat_luresol(). Once
 *    not needed any more, it must be destroyed by calling \ref
 *    igraph_sparsemat_symbolic_destroy() on it.
 * \param tol The tolerance for the numeric LU decomposition.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_lu(const igraph_sparsemat_t *A,
                        const igraph_sparsemat_symbolic_t *dis,
                        igraph_sparsemat_numeric_t *din, double tol) {
    din->numeric = cs_lu(A->cs, dis->symbolic, tol);
    if (!din->numeric) {
        IGRAPH_ERROR("Cannot do LU decomposition", IGRAPH_FAILURE);
    }
    return 0;
}

/**
 * \function igraph_sparsemat_qr
 * QR decomposition of a sparse matrix
 *
 * Numeric QR decomposition of a sparse matrix.
 * \param A The input matrix, in column-compressed format.
 * \param dis The result of the symbolic QR analysis, from the
 *    function \ref igraph_sparsemat_symbqr().
 * \param din The result of the decomposition is stored here, it can
 *    be used to solve many linear systems with the same coefficient
 *    matrix and changing right hand sides, using the \ref
 *    igraph_sparsemat_qrresol() function. Once not needed any more,
 *    one should call \ref igraph_sparsemat_numeric_destroy() on it to
 *    free the allocated memory.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_qr(const igraph_sparsemat_t *A,
                        const igraph_sparsemat_symbolic_t *dis,
                        igraph_sparsemat_numeric_t *din) {
    din->numeric = cs_qr(A->cs, dis->symbolic);
    if (!din->numeric) {
        IGRAPH_ERROR("Cannot do QR decomposition", IGRAPH_FAILURE);
    }
    return 0;
}

/**
 * \function igraph_sparsemat_luresol
 * Solve linear system using a precomputed LU decomposition
 *
 * Uses the LU decomposition of a matrix to solve linear systems.
 * \param dis The symbolic analysis of the coefficient matrix, the
 *    result of \ref igraph_sparsemat_symblu().
 * \param din The LU decomposition, the result of a call to \ref
 *    igraph_sparsemat_lu().
 * \param b A vector that defines the right hand side of the linear
 *    equation system.
 * \param res An initialized vector, the solution of the linear system
 *    is stored here.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_luresol(const igraph_sparsemat_symbolic_t *dis,
                             const igraph_sparsemat_numeric_t *din,
                             const igraph_vector_t *b,
                             igraph_vector_t *res) {
    int n = din->numeric->L->n;
    igraph_real_t *workspace;

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    workspace = igraph_Calloc(n, igraph_real_t);
    if (!workspace) {
        IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_ENOMEM);
    }
    IGRAPH_FINALLY(igraph_free, workspace);

    if (!cs_ipvec(din->numeric->pinv, VECTOR(*res), workspace, n)) {
        IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    if (!cs_lsolve(din->numeric->L, workspace)) {
        IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    if (!cs_usolve(din->numeric->U, workspace)) {
        IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    if (!cs_ipvec(dis->symbolic->q, workspace, VECTOR(*res), n)) {
        IGRAPH_ERROR("Cannot LU (re)solve sparse matrix", IGRAPH_FAILURE);
    }

    igraph_Free(workspace);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

/**
 * \function igraph_sparsemat_qrresol
 * Solve a linear system using a precomputed QR decomposition
 *
 * Solves a linear system using a QR decomposition of its coefficient
 * matrix.
 * \param dis Symbolic analysis of the coefficient matrix, the result
 *    of \ref igraph_sparsemat_symbqr().
 * \param din The QR decomposition of the coefficient matrix, the
 *    result of \ref igraph_sparsemat_qr().
 * \param b Vector, giving the right hand side of the linear equation
 *    system.
 * \param res An initialized vector, the solution is stored here. It
 *    is resized as needed.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_qrresol(const igraph_sparsemat_symbolic_t *dis,
                             const igraph_sparsemat_numeric_t *din,
                             const igraph_vector_t *b,
                             igraph_vector_t *res) {
    int n = din->numeric->L->n;
    igraph_real_t *workspace;
    int k;

    if (res != b) {
        IGRAPH_CHECK(igraph_vector_update(res, b));
    }

    workspace = igraph_Calloc(dis->symbolic ? dis->symbolic->m2 : 1,
                              igraph_real_t);
    if (!workspace) {
        IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    IGRAPH_FINALLY(igraph_free, workspace);

    if (!cs_ipvec(dis->symbolic->pinv, VECTOR(*res), workspace, n)) {
        IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    for (k = 0; k < n; k++) {
        if (!cs_happly(din->numeric->L, k, din->numeric->B[k], workspace)) {
            IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
        }
    }
    if (!cs_usolve(din->numeric->U, workspace)) {
        IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
    }
    if (!cs_ipvec(dis->symbolic->q, workspace, VECTOR(*res), n)) {
        IGRAPH_ERROR("Cannot QR (re)solve sparse matrix", IGRAPH_FAILURE);
    }

    igraph_Free(workspace);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

/**
 * \function igraph_sparsemat_symbolic_destroy
 * Deallocate memory for a symbolic decomposition
 *
 * Frees the memory allocated by \ref igraph_sparsemat_symbqr() or
 * \ref igraph_sparsemat_symblu().
 * \param dis The symbolic analysis.
 *
 * Time complexity: O(1).
 */

void igraph_sparsemat_symbolic_destroy(igraph_sparsemat_symbolic_t *dis) {
    cs_sfree(dis->symbolic);
    dis->symbolic = 0;
}

/**
 * \function igraph_sparsemat_numeric_destroy
 * Deallocate memory for a numeric decomposition
 *
 * Frees the memoty allocated by \ref igraph_sparsemat_qr() or \ref
 * igraph_sparsemat_lu().
 * \param din The LU or QR decomposition.
 *
 * Time complexity: O(1).
 */

void igraph_sparsemat_numeric_destroy(igraph_sparsemat_numeric_t *din) {
    cs_nfree(din->numeric);
    din->numeric = 0;
}

/**
 * \function igraph_matrix_as_sparsemat
 * Convert a dense matrix to a sparse matrix
 *
 * \param res An uninitialized sparse matrix, the result is stored
 *    here.
 * \param mat The dense input matrix.
 * \param tol Real scalar, the tolerance. Values closer than \p tol to
 *    zero are considered as zero, and will not be included in the
 *    sparse matrix.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the dense
 * matrix.
 */

int igraph_matrix_as_sparsemat(igraph_sparsemat_t *res,
                               const igraph_matrix_t *mat,
                               igraph_real_t tol) {
    int nrow = (int) igraph_matrix_nrow(mat);
    int ncol = (int) igraph_matrix_ncol(mat);
    int i, j, nzmax = 0;

    for (i = 0; i < nrow; i++) {
        for (j = 0; j < ncol; j++) {
            if (fabs(MATRIX(*mat, i, j)) > tol) {
                nzmax++;
            }
        }
    }

    IGRAPH_CHECK(igraph_sparsemat_init(res, nrow, ncol, nzmax));

    for (i = 0; i < nrow; i++) {
        for (j = 0; j < ncol; j++) {
            if (fabs(MATRIX(*mat, i, j)) > tol) {
                IGRAPH_CHECK(igraph_sparsemat_entry(res, i, j, MATRIX(*mat, i, j)));
            }
        }
    }

    return 0;
}

int igraph_i_sparsemat_as_matrix_cc(igraph_matrix_t *res,
                                    const igraph_sparsemat_t *spmat) {

    int nrow = (int) igraph_sparsemat_nrow(spmat);
    int ncol = (int) igraph_sparsemat_ncol(spmat);
    int *p = spmat->cs->p;
    int *i = spmat->cs->i;
    igraph_real_t *x = spmat->cs->x;
    int nzmax = spmat->cs->nzmax;
    int from = 0, to = 0;

    IGRAPH_CHECK(igraph_matrix_resize(res, nrow, ncol));
    igraph_matrix_null(res);

    while (*p < nzmax) {
        while (to < * (p + 1)) {
            MATRIX(*res, *i, from) += *x;
            to++;
            i++;
            x++;
        }
        from++;
        p++;
    }

    return 0;
}

int igraph_i_sparsemat_as_matrix_triplet(igraph_matrix_t *res,
        const igraph_sparsemat_t *spmat) {
    int nrow = (int) igraph_sparsemat_nrow(spmat);
    int ncol = (int) igraph_sparsemat_ncol(spmat);
    int *i = spmat->cs->p;
    int *j = spmat->cs->i;
    igraph_real_t *x = spmat->cs->x;
    int nz = spmat->cs->nz;
    int e;

    IGRAPH_CHECK(igraph_matrix_resize(res, nrow, ncol));
    igraph_matrix_null(res);

    for (e = 0; e < nz; e++, i++, j++, x++) {
        MATRIX(*res, *j, *i) += *x;
    }

    return 0;
}

/**
 * \function igraph_sparsemat_as_matrix
 * Convert a sparse matrix to a dense matrix
 *
 * \param res Pointer to an initialized matrix, the result is stored
 *    here. It will be resized to the required size.
 * \param spmat The input sparse matrix, in triplet or
 *    column-compressed format.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the dense
 * matrix.
 */

int igraph_sparsemat_as_matrix(igraph_matrix_t *res,
                               const igraph_sparsemat_t *spmat) {
    if (spmat->cs->nz < 0) {
        return (igraph_i_sparsemat_as_matrix_cc(res, spmat));
    } else {
        return (igraph_i_sparsemat_as_matrix_triplet(res, spmat));
    }
}

/**
 * \function igraph_sparsemat_max
 * Maximum of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \return The maximum in the input matrix, or \c IGRAPH_NEGINFINITY
 *    if the matrix has zero elements.
 *
 * Time complexity: TODO.
 */

igraph_real_t igraph_sparsemat_max(igraph_sparsemat_t *A) {
    int i, n;
    igraph_real_t *ptr;
    igraph_real_t res;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ptr = A->cs->x;
    n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    if (n == 0) {
        return IGRAPH_NEGINFINITY;
    }
    res = *ptr;
    for (i = 1; i < n; i++, ptr++) {
        if (*ptr > res) {
            res = *ptr;
        }
    }
    return res;
}

/* TODO: CC matrix don't actually need _dupl,
   because the elements are right beside each other.
   Same for max and minmax. */

/**
 * \function igraph_sparsemat_min
 * Minimum of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \return The minimum in the input matrix, or \c IGRAPH_POSINFINITY
 *    if the matrix has zero elements.
 *
 * Time complexity: TODO.
 */

igraph_real_t igraph_sparsemat_min(igraph_sparsemat_t *A) {
    int i, n;
    igraph_real_t *ptr;
    igraph_real_t res;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ptr = A->cs->x;
    n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    if (n == 0) {
        return IGRAPH_POSINFINITY;
    }
    res = *ptr;
    for (i = 1; i < n; i++, ptr++) {
        if (*ptr < res) {
            res = *ptr;
        }
    }
    return res;
}

/**
 * \function igraph_sparsemat_minmax
 * Minimum and maximum of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \param min The minimum in the input matrix is stored here, or \c
 *    IGRAPH_POSINFINITY if the matrix has zero elements.
 * \param max The maximum in the input matrix is stored here, or \c
 *    IGRAPH_NEGINFINITY if the matrix has zero elements.
 * \return Error code.
 *
 * Time complexity: TODO.
 */


int igraph_sparsemat_minmax(igraph_sparsemat_t *A,
                            igraph_real_t *min, igraph_real_t *max) {
    int i, n;
    igraph_real_t *ptr;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ptr = A->cs->x;
    n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    if (n == 0) {
        *min = IGRAPH_POSINFINITY;
        *max = IGRAPH_NEGINFINITY;
        return 0;
    }
    *min = *max = *ptr;
    for (i = 1; i < n; i++, ptr++) {
        if (*ptr > *max) {
            *max = *ptr;
        } else if (*ptr < *min) {
            *min = *ptr;
        }
    }
    return 0;
}

/**
 * \function igraph_sparsemat_count_nonzero
 * Count nonzero elements of a sparse matrix
 *
 * \param A The input matrix, column-compressed.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

long int igraph_sparsemat_count_nonzero(igraph_sparsemat_t *A) {
    int i, n;
    int res = 0;
    igraph_real_t *ptr;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ptr = A->cs->x;
    n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    if (n == 0) {
        return 0;
    }
    for (i = 0; i < n; i++, ptr++) {
        if (*ptr) {
            res++;
        }
    }
    return res;
}

/**
 * \function igraph_sparsemat_count_nonzerotol
 * Count nonzero elements of a sparse matrix, ignoring elements close to zero
 *
 * Count the number of matrix entries that are closer to zero than \p
 * tol.
 * \param The input matrix, column-compressed.
 * \param Real scalar, the tolerance.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

long int igraph_sparsemat_count_nonzerotol(igraph_sparsemat_t *A,
        igraph_real_t tol) {
    int i, n;
    int res = 0;
    igraph_real_t *ptr;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ptr = A->cs->x;
    n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    if (n == 0) {
        return 0;
    }
    for (i = 0; i < n; i++, ptr++) {
        if (*ptr < - tol || *ptr > tol) {
            res++;
        }
    }
    return res;
}

int igraph_i_sparsemat_rowsums_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pi = A->cs->i;
    double *px = A->cs->x;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_null(res);

    for (i = 0; i < A->cs->nz; i++, pi++, px++) {
        VECTOR(*res)[ *pi ] += *px;
    }

    return 0;
}

int igraph_i_sparsemat_rowsums_cc(const igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int ne = A->cs->p[A->cs->n];
    double *px = A->cs->x;
    int *pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_null(res);

    for (; pi < A->cs->i + ne; pi++, px++) {
        VECTOR(*res)[ *pi ] += *px;
    }

    return 0;
}

/**
 * \function igraph_sparsemat_rowsums
 * Row-wise sums.
 *
 * \param A The input matrix, in triplet or column-compressed format.
 * \param res An initialized vector, the result is stored here. It
 *    will be resized as needed.
 * \return Error code.
 *
 * Time complexity: O(nz), the number of non-zero elements.
 */

int igraph_sparsemat_rowsums(const igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_rowsums_triplet(A, res);
    } else {
        return igraph_i_sparsemat_rowsums_cc(A, res);
    }
}

int igraph_i_sparsemat_rowmins_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pi = A->cs->i;
    double *px = A->cs->x;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_fill(res, inf);

    for (i = 0; i < A->cs->nz; i++, pi++, px++) {
        if (*px < VECTOR(*res)[ *pi ]) {
            VECTOR(*res)[ *pi ] = *px;
        }
    }

    return 0;
}

int igraph_i_sparsemat_rowmins_cc(igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int ne;
    double *px;
    int *pi;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ne = A->cs->p[A->cs->n];
    px = A->cs->x;
    pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_fill(res, inf);

    for (; pi < A->cs->i + ne; pi++, px++) {
        if (*px < VECTOR(*res)[ *pi ]) {
            VECTOR(*res)[ *pi ] = *px;
        }
    }

    return 0;
}

int igraph_sparsemat_rowmins(igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_rowmins_triplet(A, res);
    } else {
        return igraph_i_sparsemat_rowmins_cc(A, res);
    }
}


int igraph_i_sparsemat_rowmaxs_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pi = A->cs->i;
    double *px = A->cs->x;
    double inf = IGRAPH_NEGINFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_fill(res, inf);

    for (i = 0; i < A->cs->nz; i++, pi++, px++) {
        if (*px > VECTOR(*res)[ *pi ]) {
            VECTOR(*res)[ *pi ] = *px;
        }
    }

    return 0;
}

int igraph_i_sparsemat_rowmaxs_cc(igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int ne;
    double *px;
    int *pi;
    double inf = IGRAPH_NEGINFINITY;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    ne = A->cs->p[A->cs->n];
    px = A->cs->x;
    pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    igraph_vector_fill(res, inf);

    for (; pi < A->cs->i + ne; pi++, px++) {
        if (*px > VECTOR(*res)[ *pi ]) {
            VECTOR(*res)[ *pi ] = *px;
        }
    }

    return 0;
}

int igraph_sparsemat_rowmaxs(igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_rowmaxs_triplet(A, res);
    } else {
        return igraph_i_sparsemat_rowmaxs_cc(A, res);
    }
}

int igraph_i_sparsemat_colmins_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pp = A->cs->p;
    double *px = A->cs->x;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->n));
    igraph_vector_fill(res, inf);

    for (i = 0; i < A->cs->nz; i++, pp++, px++) {
        if (*px < VECTOR(*res)[ *pp ]) {
            VECTOR(*res)[ *pp ] = *px;
        }
    }

    return 0;
}

int igraph_i_sparsemat_colmins_cc(igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int n;
    double *px;
    int *pp;
    int *pi;
    double *pr;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    n = A->cs->n;
    px = A->cs->x;
    pp = A->cs->p;
    pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, n));
    igraph_vector_fill(res, inf);
    pr = VECTOR(*res);

    for (; pp < A->cs->p + n; pp++, pr++) {
        for (; pi < A->cs->i + * (pp + 1); pi++, px++) {
            if (*px < *pr) {
                *pr = *px;
            }
        }
    }
    return 0;
}

int igraph_sparsemat_colmins(igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_colmins_triplet(A, res);
    } else {
        return igraph_i_sparsemat_colmins_cc(A, res);
    }
}

int igraph_i_sparsemat_colmaxs_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pp = A->cs->p;
    double *px = A->cs->x;
    double inf = IGRAPH_NEGINFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->n));
    igraph_vector_fill(res, inf);

    for (i = 0; i < A->cs->nz; i++, pp++, px++) {
        if (*px > VECTOR(*res)[ *pp ]) {
            VECTOR(*res)[ *pp ] = *px;
        }
    }

    return 0;
}

int igraph_i_sparsemat_colmaxs_cc(igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int n;
    double *px;
    int *pp;
    int *pi;
    double *pr;
    double inf = IGRAPH_NEGINFINITY;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    n = A->cs->n;
    px = A->cs->x;
    pp = A->cs->p;
    pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, n));
    igraph_vector_fill(res, inf);
    pr = VECTOR(*res);

    for (; pp < A->cs->p + n; pp++, pr++) {
        for (; pi < A->cs->i + * (pp + 1); pi++, px++) {
            if (*px > *pr) {
                *pr = *px;
            }
        }
    }
    return 0;
}

int igraph_sparsemat_colmaxs(igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_colmaxs_triplet(A, res);
    } else {
        return igraph_i_sparsemat_colmaxs_cc(A, res);
    }
}

int igraph_i_sparsemat_which_min_rows_triplet(igraph_sparsemat_t *A,
        igraph_vector_t *res,
        igraph_vector_int_t *pos) {
    int i;
    int *pi = A->cs->i;
    int *pp = A->cs->p;
    double *px = A->cs->x;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    IGRAPH_CHECK(igraph_vector_int_resize(pos, A->cs->m));
    igraph_vector_fill(res, inf);
    igraph_vector_int_null(pos);

    for (i = 0; i < A->cs->nz; i++, pi++, px++, pp++) {
        if (*px < VECTOR(*res)[ *pi ]) {
            VECTOR(*res)[ *pi ] = *px;
            VECTOR(*pos)[ *pi ] = *pp;
        }
    }

    return 0;
}

int igraph_i_sparsemat_which_min_rows_cc(igraph_sparsemat_t *A,
        igraph_vector_t *res,
        igraph_vector_int_t *pos) {
    int n;
    double *px;
    int *pp;
    int *pi;
    double inf = IGRAPH_INFINITY;
    int j;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    n = A->cs->n;
    px = A->cs->x;
    pp = A->cs->p;
    pi = A->cs->i;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->m));
    IGRAPH_CHECK(igraph_vector_int_resize(pos, A->cs->m));
    igraph_vector_fill(res, inf);
    igraph_vector_int_null(pos);

    for (j = 0; pp < A->cs->p + n; pp++, j++) {
        for (; pi < A->cs->i + * (pp + 1); pi++, px++) {
            if (*px < VECTOR(*res)[ *pi ]) {
                VECTOR(*res)[ *pi ] = *px;
                VECTOR(*pos)[ *pi ] = j;
            }
        }
    }

    return 0;
}

int igraph_sparsemat_which_min_rows(igraph_sparsemat_t *A,
                                    igraph_vector_t *res,
                                    igraph_vector_int_t *pos) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_which_min_rows_triplet(A, res, pos);
    } else {
        return igraph_i_sparsemat_which_min_rows_cc(A, res, pos);
    }
}

int igraph_i_sparsemat_which_min_cols_triplet(igraph_sparsemat_t *A,
        igraph_vector_t *res,
        igraph_vector_int_t *pos) {

    int i;
    int *pi = A->cs->i;
    int *pp = A->cs->p;
    double *px = A->cs->x;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->n));
    IGRAPH_CHECK(igraph_vector_int_resize(pos, A->cs->n));
    igraph_vector_fill(res, inf);
    igraph_vector_int_null(pos);

    for (i = 0; i < A->cs->nz; i++, pi++, pp++, px++) {
        if (*px < VECTOR(*res)[ *pp ]) {
            VECTOR(*res)[ *pp ] = *px;
            VECTOR(*pos)[ *pp ] = *pi;
        }
    }

    return 0;
}

int igraph_i_sparsemat_which_min_cols_cc(igraph_sparsemat_t *A,
        igraph_vector_t *res,
        igraph_vector_int_t *pos) {
    int n, j, p;
    double *px;
    double *pr;
    int *ppos;
    double inf = IGRAPH_INFINITY;

    IGRAPH_CHECK(igraph_sparsemat_dupl(A));

    n = A->cs->n;
    px = A->cs->x;

    IGRAPH_CHECK(igraph_vector_resize(res, n));
    igraph_vector_fill(res, inf);
    pr = VECTOR(*res);
    IGRAPH_CHECK(igraph_vector_int_resize(pos, n));
    igraph_vector_int_null(pos);
    ppos = VECTOR(*pos);

    for (j = 0; j < A->cs->n; j++, pr++, ppos++) {
        for (p = A->cs->p[j]; p < A->cs->p[j + 1]; p++, px++) {
            if (*px < *pr) {
                *pr = *px;
                *ppos = A->cs->i[p];
            }
        }
    }
    return 0;
}

int igraph_sparsemat_which_min_cols(igraph_sparsemat_t *A,
                                    igraph_vector_t *res,
                                    igraph_vector_int_t *pos) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_which_min_cols_triplet(A, res, pos);
    } else {
        return igraph_i_sparsemat_which_min_cols_cc(A, res, pos);
    }
}

int igraph_i_sparsemat_colsums_triplet(const igraph_sparsemat_t *A,
                                       igraph_vector_t *res) {
    int i;
    int *pp = A->cs->p;
    double *px = A->cs->x;

    IGRAPH_CHECK(igraph_vector_resize(res, A->cs->n));
    igraph_vector_null(res);

    for (i = 0; i < A->cs->nz; i++, pp++, px++) {
        VECTOR(*res)[ *pp ] += *px;
    }

    return 0;
}

int igraph_i_sparsemat_colsums_cc(const igraph_sparsemat_t *A,
                                  igraph_vector_t *res) {
    int n = A->cs->n;
    double *px = A->cs->x;
    int *pp = A->cs->p;
    int *pi = A->cs->i;
    double *pr;

    IGRAPH_CHECK(igraph_vector_resize(res, n));
    igraph_vector_null(res);
    pr = VECTOR(*res);

    for (; pp < A->cs->p + n; pp++, pr++) {
        for (; pi < A->cs->i + * (pp + 1); pi++, px++) {
            *pr += *px;
        }
    }
    return 0;
}

/**
 * \function igraph_sparsemat_colsums
 * Column-wise sums
 *
 * \param A The input matrix, in triplet or column-compressed format.
 * \param res An initialized vector, the result is stored here. It
 *    will be resized as needed.
 * \return Error code.
 *
 * Time complexity: O(nz) for triplet matrices, O(nz+n) for
 * column-compressed ones, nz is the number of non-zero elements, n is
 * the number of columns.
 */

int igraph_sparsemat_colsums(const igraph_sparsemat_t *A,
                             igraph_vector_t *res) {
    if (igraph_sparsemat_is_triplet(A)) {
        return igraph_i_sparsemat_colsums_triplet(A, res);
    } else {
        return igraph_i_sparsemat_colsums_cc(A, res);
    }
}

/**
 * \function igraph_sparsemat_scale
 * Scale a sparse matrix
 *
 * Multiplies all elements of a sparse matrix, by the given scalar.
 * \param A The input matrix.
 * \param by The scaling factor.
 * \return Error code.
 *
 * Time complexity: O(nz), the number of non-zero elements in the
 * matrix.
 */

int igraph_sparsemat_scale(igraph_sparsemat_t *A, igraph_real_t by) {

    double *px = A->cs->x;
    int n = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    double *stop = px + n;

    for (; px < stop; px++) {
        *px *= by;
    }

    return 0;
}

/**
 * \function igraph_sparsemat_add_rows
 * Add rows to a sparse matrix
 *
 * The current matrix elements are retained and all elements in the
 * new rows are zero.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param n The number of rows to add.
 * \return Error code.
 *
 * Time complexity: O(1).
 */

int igraph_sparsemat_add_rows(igraph_sparsemat_t *A, long int n) {
    A->cs->m += n;
    return 0;
}

/**
 * \function igraph_sparsemat_add_cols
 * Add columns to a sparse matrix
 *
 * The current matrix elements are retained, and all elements in the
 * new columns are zero.
 * \param A The input matrix, in triplet or column-compressed format.
 * \param n The number of columns to add.
 * \return Error code.
 *
 * Time complexity: TODO.
 */

int igraph_sparsemat_add_cols(igraph_sparsemat_t *A, long int n) {
    if (igraph_sparsemat_is_triplet(A)) {
        A->cs->n += n;
    } else {
        int *newp = realloc(A->cs->p, sizeof(int) * (size_t) (A->cs->n + n + 1));
        int i;
        if (!newp) {
            IGRAPH_ERROR("Cannot add columns to sparse matrix", IGRAPH_ENOMEM);
        }
        if (newp != A->cs->p) {
            A->cs->p = newp;
        }
        for (i = A->cs->n + 1; i < A->cs->n + n + 1; i++) {
            A->cs->p[i] = A->cs->p[i - 1];
        }
        A->cs->n += n;
    }
    return 0;
}

/**
 * \function igraph_sparsemat_resize
 * Resize a sparse matrix
 *
 * This function resizes a sparse matrix. The resized sparse matrix
 * will be empty.
 *
 * \param A The initialized sparse matrix to resize.
 * \param nrow The new number of rows.
 * \param ncol The new number of columns.
 * \param nzmax The new maximum number of elements.
 * \return Error code.
 *
 * Time complexity: O(nzmax), the maximum number of non-zero elements.
 */

int igraph_sparsemat_resize(igraph_sparsemat_t *A, long int nrow,
                            long int ncol, int nzmax) {

    if (A->cs->nz < 0) {
        igraph_sparsemat_t tmp;
        IGRAPH_CHECK(igraph_sparsemat_init(&tmp, (int) nrow, (int) ncol, nzmax));
        igraph_sparsemat_destroy(A);
        *A = tmp;
    } else {
        IGRAPH_CHECK(igraph_sparsemat_realloc(A, nzmax));
        A->cs->m = (int) nrow;
        A->cs->n = (int) ncol;
        A->cs->nz = 0;
    }
    return 0;
}

int igraph_sparsemat_nonzero_storage(const igraph_sparsemat_t *A) {
    if (A->cs->nz < 0) {
        return A->cs->p[A->cs->n];
    } else {
        return A->cs->nz;
    }
}

int igraph_sparsemat_getelements(const igraph_sparsemat_t *A,
                                 igraph_vector_int_t *i,
                                 igraph_vector_int_t *j,
                                 igraph_vector_t *x) {
    int nz = A->cs->nz;
    if (nz < 0) {
        nz = A->cs->p[A->cs->n];
        IGRAPH_CHECK(igraph_vector_int_resize(i, nz));
        IGRAPH_CHECK(igraph_vector_int_resize(j, A->cs->n + 1));
        IGRAPH_CHECK(igraph_vector_resize(x, nz));
        memcpy(VECTOR(*i), A->cs->i, (size_t) nz * sizeof(int));
        memcpy(VECTOR(*j), A->cs->p, (size_t) (A->cs->n + 1) * sizeof(int));
        memcpy(VECTOR(*x), A->cs->x, (size_t) nz * sizeof(igraph_real_t));
    } else {
        IGRAPH_CHECK(igraph_vector_int_resize(i, nz));
        IGRAPH_CHECK(igraph_vector_int_resize(j, nz));
        IGRAPH_CHECK(igraph_vector_resize(x, nz));
        memcpy(VECTOR(*i), A->cs->i, (size_t) nz * sizeof(int));
        memcpy(VECTOR(*j), A->cs->p, (size_t) nz * sizeof(int));
        memcpy(VECTOR(*x), A->cs->x, (size_t) nz * sizeof(igraph_real_t));
    }
    return 0;
}

int igraph_sparsemat_scale_rows(igraph_sparsemat_t *A,
                                const igraph_vector_t *fact) {
    int *i = A->cs->i;
    igraph_real_t *x = A->cs->x;
    int no_of_edges = A->cs->nz < 0 ? A->cs->p[A->cs->n] : A->cs->nz;
    int e;

    for (e = 0; e < no_of_edges; e++, x++, i++) {
        igraph_real_t f = VECTOR(*fact)[*i];
        (*x) *= f;
    }

    return 0;
}

int igraph_i_sparsemat_scale_cols_cc(igraph_sparsemat_t *A,
                                     const igraph_vector_t *fact) {
    int *i = A->cs->i;
    igraph_real_t *x = A->cs->x;
    int no_of_edges = A->cs->p[A->cs->n];
    int e;
    int c = 0;        /* actual column */

    for (e = 0; e < no_of_edges; e++, x++, i++) {
        igraph_real_t f;
        while (c < A->cs->n && A->cs->p[c + 1] == e) {
            c++;
        }
        f = VECTOR(*fact)[c];
        (*x) *= f;
    }

    return 0;
}

int igraph_i_sparsemat_scale_cols_triplet(igraph_sparsemat_t *A,
        const igraph_vector_t *fact) {
    int *j = A->cs->p;
    igraph_real_t *x = A->cs->x;
    int no_of_edges = A->cs->nz;
    int e;

    for (e = 0; e < no_of_edges; e++, x++, j++) {
        igraph_real_t f = VECTOR(*fact)[*j];
        (*x) *= f;
    }

    return 0;
}

int igraph_sparsemat_scale_cols(igraph_sparsemat_t *A,
                                const igraph_vector_t *fact) {
    if (A->cs->nz < 0) {
        return igraph_i_sparsemat_scale_cols_cc(A, fact);
    } else {
        return igraph_i_sparsemat_scale_cols_triplet(A, fact);
    }
}

int igraph_sparsemat_multiply_by_dense(const igraph_sparsemat_t *A,
                                       const igraph_matrix_t *B,
                                       igraph_matrix_t *res) {

    int m = (int) igraph_sparsemat_nrow(A);
    int n = (int) igraph_sparsemat_ncol(A);
    int p = (int) igraph_matrix_ncol(B);
    int i;

    if (igraph_matrix_nrow(B) != n) {
        IGRAPH_ERROR("Invalid dimensions in sparse-dense matrix product",
                     IGRAPH_EINVAL);
    }

    IGRAPH_CHECK(igraph_matrix_resize(res, m, p));
    igraph_matrix_null(res);

    for (i = 0; i < p; i++) {
        if (!(cs_gaxpy(A->cs, &MATRIX(*B, 0, i), &MATRIX(*res, 0, i)))) {
            IGRAPH_ERROR("Cannot perform sparse-dense matrix multiplication",
                         IGRAPH_FAILURE);
        }
    }

    return 0;
}

int igraph_sparsemat_dense_multiply(const igraph_matrix_t *A,
                                    const igraph_sparsemat_t *B,
                                    igraph_matrix_t *res) {
    int m = (int) igraph_matrix_nrow(A);
    int n = (int) igraph_matrix_ncol(A);
    int p = (int) igraph_sparsemat_ncol(B);
    int r, c;
    int *Bp = B->cs->p;

    if (igraph_sparsemat_nrow(B) != n) {
        IGRAPH_ERROR("Invalid dimensions in dense-sparse matrix product",
                     IGRAPH_EINVAL);
    }

    if (!igraph_sparsemat_is_cc(B)) {
        IGRAPH_ERROR("Dense-sparse product is only implemented for "
                     "column-compressed sparse matrices", IGRAPH_EINVAL);
    }

    IGRAPH_CHECK(igraph_matrix_resize(res, m, p));
    igraph_matrix_null(res);

    for (c = 0; c < p; c++) {
        for (r = 0; r < m; r++) {
            int idx = *Bp;
            while (idx < * (Bp + 1)) {
                MATRIX(*res, r, c) += MATRIX(*A, r, B->cs->i[idx]) * B->cs->x[idx];
                idx++;
            }
        }
        Bp++;
    }

    return 0;
}

int igraph_i_sparsemat_view(igraph_sparsemat_t *A, int nzmax, int m, int n,
                            int *p, int *i, double *x, int nz) {

    A->cs = cs_calloc(1, sizeof(cs_di));
    A->cs->nzmax = nzmax;
    A->cs->m = m;
    A->cs->n = n;
    A->cs->p = p;
    A->cs->i = i;
    A->cs->x = x;
    A->cs->nz = nz;

    return 0;
}

int igraph_sparsemat_sort(const igraph_sparsemat_t *A,
                          igraph_sparsemat_t *sorted) {

    igraph_sparsemat_t tmp;

    IGRAPH_CHECK(igraph_sparsemat_transpose(A, &tmp, /*values=*/ 1));
    IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
    IGRAPH_CHECK(igraph_sparsemat_transpose(&tmp, sorted, /*values=*/ 1));
    igraph_sparsemat_destroy(&tmp);
    IGRAPH_FINALLY_CLEAN(1);

    return 0;
}

int igraph_sparsemat_getelements_sorted(const igraph_sparsemat_t *A,
                                        igraph_vector_int_t *i,
                                        igraph_vector_int_t *j,
                                        igraph_vector_t *x) {
    if (A->cs->nz < 0) {
        igraph_sparsemat_t tmp;
        IGRAPH_CHECK(igraph_sparsemat_sort(A, &tmp));
        IGRAPH_FINALLY(igraph_sparsemat_destroy, &tmp);
        IGRAPH_CHECK(igraph_sparsemat_getelements(&tmp, i, j, x));
        igraph_sparsemat_destroy(&tmp);
        IGRAPH_FINALLY_CLEAN(1);
    } else {
        IGRAPH_CHECK(igraph_sparsemat_getelements(A, i, j, x));
    }

    return 0;
}

int igraph_sparsemat_nzmax(const igraph_sparsemat_t *A) {
    return A->cs->nzmax;
}

int igraph_sparsemat_neg(igraph_sparsemat_t *A) {
    int i, nz = A->cs->nz == -1 ? A->cs->p[A->cs->n] : A->cs->nz;
    igraph_real_t *px = A->cs->x;

    for (i = 0; i < nz; i++, px++) {
        *px = - (*px);
    }

    return 0;
}

int igraph_sparsemat_iterator_init(igraph_sparsemat_iterator_t *it,
                                   igraph_sparsemat_t *sparsemat) {

    it->mat = sparsemat;
    igraph_sparsemat_iterator_reset(it);
    return 0;
}

int igraph_sparsemat_iterator_reset(igraph_sparsemat_iterator_t *it) {
    it->pos = 0;
    if (!igraph_sparsemat_is_triplet(it->mat)) {
        it->col = 0;
        while (it->col < it->mat->cs->n &&
               it->mat->cs->p[it->col + 1] == it->pos) {
            it->col ++;
        }
    }
    return 0;
}

igraph_bool_t
igraph_sparsemat_iterator_end(const igraph_sparsemat_iterator_t *it) {
    int nz = it->mat->cs->nz == -1 ? it->mat->cs->p[it->mat->cs->n] :
             it->mat->cs->nz;
    return it->pos >= nz;
}

int igraph_sparsemat_iterator_row(const igraph_sparsemat_iterator_t *it) {
    return it->mat->cs->i[it->pos];
}

int igraph_sparsemat_iterator_col(const igraph_sparsemat_iterator_t *it) {
    if (igraph_sparsemat_is_triplet(it->mat)) {
        return it->mat->cs->p[it->pos];
    } else {
        return it->col;
    }
}

igraph_real_t
igraph_sparsemat_iterator_get(const igraph_sparsemat_iterator_t *it) {
    return it->mat->cs->x[it->pos];
}

int igraph_sparsemat_iterator_next(igraph_sparsemat_iterator_t *it) {
    it->pos += 1;
    while (it->col < it->mat->cs->n &&
           it->mat->cs->p[it->col + 1] == it->pos) {
        it->col++;
    }
    return it->pos;
}

int igraph_sparsemat_iterator_idx(const igraph_sparsemat_iterator_t *it) {
    return it->pos;
}