/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2003-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 "igraph_memory.h"
#include "igraph_random.h"
#include "igraph_error.h"

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

/**
 * \section about_igraph_matrix_t_objects About \type igraph_matrix_t objects
 *
 * <para>This type is just an interface to \type igraph_vector_t.</para>
 *
 * <para>The \type igraph_matrix_t type usually stores n
 * elements in O(n) space, but not always. See the documentation of
 * the vector type.</para>
 */

/**
 * \section igraph_matrix_constructor_and_destructor Matrix constructors and
 * destructors
 */

/**
 * \ingroup matrix
 * \function igraph_matrix_init
 * \brief Initializes a matrix.
 *
 * </para><para>
 * Every 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_matrix_destroy().
 * \param m Pointer to a not yet initialized 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: usually O(n),
 * n is the
 * number of elements in the matrix.
 */

int FUNCTION(igraph_matrix, init)(TYPE(igraph_matrix) *m, long int nrow, long int ncol) {
    int ret1;
    ret1 = FUNCTION(igraph_vector, init)(&m->data, nrow * ncol);
    m->nrow = nrow;
    m->ncol = ncol;
    return ret1;
}

const TYPE(igraph_matrix) *FUNCTION(igraph_matrix, view)(const TYPE(igraph_matrix) *m,
        const BASE *data,
        long int nrow,
        long int ncol) {
    TYPE(igraph_matrix) *m2 = (TYPE(igraph_matrix)*)m;
    FUNCTION(igraph_vector, view)(&m2->data, data, nrow * ncol);
    m2->nrow = nrow;
    m2->ncol = ncol;
    return m;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_destroy
 * \brief Destroys a matrix object.
 *
 * </para><para>
 * This function frees all the memory allocated for a 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 FUNCTION(igraph_matrix, destroy)(TYPE(igraph_matrix) *m) {
    FUNCTION(igraph_vector, destroy)(&m->data);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_capacity
 * \brief Returns the number of elements allocated for a matrix.
 *
 * Note that this might be different from the size of the matrix (as
 * queried by \ref igraph_matrix_size(), and specifies how many elements
 * the matrix can hold, without reallocation.
 * \param v Pointer to the (previously initialized) matrix object
 *          to query.
 * \return The allocated capacity.
 *
 * \sa \ref igraph_matrix_size(), \ref igraph_matrix_nrow(),
 * \ref igraph_matrix_ncol().
 *
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_matrix, capacity)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, capacity)(&m->data);
}


/**
 * \section igraph_matrix_accessing_elements Accessing elements of a matrix
 */

/**
 * \ingroup matrix
 * \function igraph_matrix_resize
 * \brief Resizes a matrix.
 *
 * </para><para>
 * This function resizes a matrix by adding more elements to it.
 * The matrix contains arbitrary data after resizing it.
 * That is, after calling this function you cannot expect that element
 * (i,j) in the matrix remains the
 * same as before.
 * \param m Pointer to an already initialized 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(1) if the
 * matrix gets smaller, usually O(n)
 * if it gets larger, n is the
 * number of elements in the resized matrix.
 */

int FUNCTION(igraph_matrix, resize)(TYPE(igraph_matrix) *m, long int nrow, long int ncol) {
    FUNCTION(igraph_vector, resize)(&m->data, nrow * ncol);
    m->nrow = nrow;
    m->ncol = ncol;
    return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_resize_min
 * \brief Deallocates unused memory for a matrix.
 *
 * </para><para>
 * Note that this function might fail if there is not enough memory
 * available.
 *
 * </para><para>
 * Also note, that this function leaves the matrix intact, i.e.
 * it does not destroy any of the elements. However, usually it involves
 * copying the matrix in memory.
 * \param m Pointer to an initialized matrix.
 * \return Error code.
 *
 * \sa \ref igraph_matrix_resize().
 *
 * Time complexity: operating system dependent.
 */

int FUNCTION(igraph_matrix, resize_min)(TYPE(igraph_matrix) *m) {
    TYPE(igraph_vector) tmp;
    long int size = FUNCTION(igraph_matrix, size)(m);
    long int capacity = FUNCTION(igraph_matrix, capacity)(m);
    if (size == capacity) {
        return 0;
    }

    IGRAPH_CHECK(FUNCTION(igraph_vector, init)(&tmp, size));
    FUNCTION(igraph_vector, update)(&tmp, &m->data);
    FUNCTION(igraph_vector, destroy)(&m->data);
    m->data = tmp;

    return 0;
}


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

long int FUNCTION(igraph_matrix, size)(const TYPE(igraph_matrix) *m) {
    return (m->nrow) * (m->ncol);
}

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

long int FUNCTION(igraph_matrix, nrow)(const TYPE(igraph_matrix) *m) {
    return m->nrow;
}

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

long int FUNCTION(igraph_matrix, ncol)(const TYPE(igraph_matrix) *m) {
    return m->ncol;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_copy_to
 * \brief Copies a 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.
 * \param m Pointer to an initialized 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.
 */

void FUNCTION(igraph_matrix, copy_to)(const TYPE(igraph_matrix) *m, BASE *to) {
    FUNCTION(igraph_vector, copy_to)(&m->data, to);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_null
 * \brief Sets all elements in a matrix to zero.
 *
 * \param m Pointer to an initialized matrix object.
 *
 * Time complexity: O(n),
 * n is the number of  elements in
 * the matrix.
 */

void FUNCTION(igraph_matrix, null)(TYPE(igraph_matrix) *m) {
    FUNCTION(igraph_vector, null)(&m->data);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_add_cols
 * \brief Adds columns to a matrix.
 * \param m The matrix object.
 * \param n The number of columns to add.
 * \return Error code, \c IGRAPH_ENOMEM if there is
 *   not enough memory to perform the operation.
 *
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix, add_cols)(TYPE(igraph_matrix) *m, long int n) {
    FUNCTION(igraph_matrix, resize)(m, m->nrow, m->ncol + n);
    return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_add_rows
 * \brief Adds rows to a matrix.
 * \param m The matrix object.
 * \param n The number of rows to add.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *   isn't enough memory for the operation.
 *
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix, add_rows)(TYPE(igraph_matrix) *m, long int n) {
    long int i;
    FUNCTION(igraph_vector, resize)(&m->data, (m->ncol) * (m->nrow + n));
    for (i = m->ncol - 1; i >= 0; i--) {
        FUNCTION(igraph_vector, move_interval2)(&m->data, (m->nrow)*i, (m->nrow) * (i + 1),
                                                (m->nrow + n)*i);
    }
    m->nrow += n;
    return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_remove_col
 * \brief Removes a column from a matrix.
 *
 * \param m The matrix object.
 * \param col The column to remove.
 * \return Error code, always returns with success.
 *
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix, remove_col)(TYPE(igraph_matrix) *m, long int col) {
    FUNCTION(igraph_vector, remove_section)(&m->data, (m->nrow)*col, (m->nrow) * (col + 1));
    m->ncol--;
    return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_permdelete_rows
 * \brief Removes rows from a matrix (for internal use).
 *
 * Time complexity: linear with the number of elements of the original
 * matrix.
 */

int FUNCTION(igraph_matrix, permdelete_rows)(TYPE(igraph_matrix) *m, long int *index, long int nremove) {
    long int i, j;
    for (j = 0; j < m->nrow; j++) {
        if (index[j] != 0) {
            for (i = 0; i < m->ncol; i++) {
                MATRIX(*m, index[j] - 1, i) = MATRIX(*m, j, i);
            }
        }
    }
    /* Remove unnecessary elements from the end of each column */
    for (i = 0; i < m->ncol; i++)
        FUNCTION(igraph_vector, remove_section)(&m->data,
                                                (i + 1) * (m->nrow - nremove), (i + 1) * (m->nrow - nremove) + nremove);
    FUNCTION(igraph_matrix, resize)(m, m->nrow - nremove, m->ncol);

    return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_delete_rows_neg
 * \brief Removes columns from a matrix (for internal use).
 *
 * Time complexity: linear with the number of elements of the original
 * matrix.
 */

int FUNCTION(igraph_matrix, delete_rows_neg)(TYPE(igraph_matrix) *m,
        const igraph_vector_t *neg, long int nremove) {
    long int i, j, idx = 0;
    for (i = 0; i < m->ncol; i++) {
        for (j = 0; j < m->nrow; j++) {
            if (VECTOR(*neg)[j] >= 0) {
                MATRIX(*m, idx++, i) = MATRIX(*m, j, i);
            }
        }
        idx = 0;
    }
    FUNCTION(igraph_matrix, resize)(m, m->nrow - nremove, m->ncol);

    return 0;
}

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

int FUNCTION(igraph_matrix, copy)(TYPE(igraph_matrix) *to, const TYPE(igraph_matrix) *from) {
    to->nrow = from->nrow;
    to->ncol = from->ncol;
    return FUNCTION(igraph_vector, copy)(&to->data, &from->data);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_max
 *
 * Returns the maximal element of a matrix.
 * \param m The matrix object.
 * \return The maximum element. For empty matrix the returned value is
 * undefined.
 *
 * Added in version 0.2.</para><para>
 *
 * Time complexity: O(n), the number of elements in the matrix.
 */

igraph_real_t FUNCTION(igraph_matrix, max)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, max)(&m->data);
}

#endif

/**
 * \function igraph_matrix_scale
 *
 * Multiplies each element of the matrix by a constant.
 * \param m The matrix.
 * \param by The constant.
 *
 * Added in version 0.2.</para><para>
 *
 * Time complexity: O(n), the number of elements in the matrix.
 */

void FUNCTION(igraph_matrix, scale)(TYPE(igraph_matrix) *m, BASE by) {
    FUNCTION(igraph_vector, scale)(&m->data, by);
}

/**
 * \function igraph_matrix_select_rows
 * \brief Select some rows of a matrix.
 *
 * This function selects some rows of a matrix and returns them in a
 * new matrix. The result matrix should be initialized before calling
 * the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param rows Vector; it contains the row indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 *
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix, select_rows)(const TYPE(igraph_matrix) *m,
        TYPE(igraph_matrix) *res,
        const igraph_vector_t *rows) {
    long int norows = igraph_vector_size(rows);
    long int i, j, ncols = FUNCTION(igraph_matrix, ncol)(m);

    IGRAPH_CHECK(FUNCTION(igraph_matrix, resize)(res, norows, ncols));
    for (i = 0; i < norows; i++) {
        for (j = 0; j < ncols; j++) {
            MATRIX(*res, i, j) = MATRIX(*m, (long int)VECTOR(*rows)[i], j);
        }
    }

    return 0;
}

/**
 * \function igraph_matrix_select_rows_cols
 * \brief Select some rows and columns of a matrix.
 *
 * This function selects some rows and columns of a matrix and returns
 * them in a new matrix. The result matrix should be initialized before
 * calling the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param rows Vector; it contains the row indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \param cols Vector; it contains the column indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 *
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix, select_rows_cols)(const TYPE(igraph_matrix) *m,
        TYPE(igraph_matrix) *res,
        const igraph_vector_t *rows,
        const igraph_vector_t *cols) {
    long int nrows = igraph_vector_size(rows);
    long int ncols = igraph_vector_size(cols);
    long int i, j;

    IGRAPH_CHECK(FUNCTION(igraph_matrix, resize)(res, nrows, ncols));
    for (i = 0; i < nrows; i++) {
        for (j = 0; j < ncols; j++) {
            MATRIX(*res, i, j) = MATRIX(*m, (long int)VECTOR(*rows)[i],
                                        (long int)VECTOR(*cols)[j]);
        }
    }

    return 0;
}

/**
 * \function igraph_matrix_get_col
 * \brief Select a column.
 *
 * Extract a column of a matrix and return it as a vector.
 * \param m The input matrix.
 * \param res The result will we stored in this vector. It should be
 *   initialized and will be resized as needed.
 * \param index The index of the column to select.
 * \return Error code.
 *
 * Time complexity: O(n), the number of rows in the matrix.
 */

int FUNCTION(igraph_matrix, get_col)(const TYPE(igraph_matrix) *m,
                                     TYPE(igraph_vector) *res,
                                     long int index) {
    long int nrow = FUNCTION(igraph_matrix, nrow)(m);

    if (index >= m->ncol) {
        IGRAPH_ERROR("Index out of range for selecting matrix column", IGRAPH_EINVAL);
    }
    IGRAPH_CHECK(FUNCTION(igraph_vector, get_interval)(&m->data, res,
                 nrow * index, nrow * (index + 1)));
    return 0;
}

/**
 * \function igraph_matrix_sum
 * \brief Sum of elements.
 *
 * Returns the sum of the elements of a matrix.
 * \param m The input matrix.
 * \return The sum of the elements.
 *
 * Time complexity: O(mn), the number of elements in the matrix.
 */

BASE FUNCTION(igraph_matrix, sum)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, sum)(&m->data);
}

/**
 * \function igraph_matrix_all_e
 * \brief Are all elements equal?
 *
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all
 *    equal to the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 *
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix, all_e)(const TYPE(igraph_matrix) *lhs,
        const TYPE(igraph_matrix) *rhs) {
    return lhs->ncol == rhs->ncol && lhs->nrow == rhs->nrow &&
           FUNCTION(igraph_vector, all_e)(&lhs->data, &rhs->data);
}

igraph_bool_t
FUNCTION(igraph_matrix, is_equal)(const TYPE(igraph_matrix) *lhs,
                                  const TYPE(igraph_matrix) *rhs) {
    return FUNCTION(igraph_matrix, all_e)(lhs, rhs);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_all_l
 * \brief Are all elements less?
 *
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all
 *    less than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 *
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix, all_l)(const TYPE(igraph_matrix) *lhs,
        const TYPE(igraph_matrix) *rhs) {
    return lhs->ncol == rhs->ncol && lhs->nrow == rhs->nrow &&
           FUNCTION(igraph_vector, all_l)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_g
 * \brief Are all elements greater?
 *
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all
 *    greater than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 *
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix, all_g)(const TYPE(igraph_matrix) *lhs,
        const TYPE(igraph_matrix) *rhs) {
    return lhs->ncol == rhs->ncol && lhs->nrow == rhs->nrow &&
           FUNCTION(igraph_vector, all_g)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_le
 * \brief Are all elements less or equal?
 *
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all
 *    less than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the dimensions of the matrices
 *    don't match.
 *
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t
FUNCTION(igraph_matrix, all_le)(const TYPE(igraph_matrix) *lhs,
                                const TYPE(igraph_matrix) *rhs) {
    return lhs->ncol == rhs->ncol && lhs->nrow == rhs->nrow &&
           FUNCTION(igraph_vector, all_le)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_ge
 * \brief Are all elements greater or equal?
 *
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all
 *    greater than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the dimensions of the matrices
 *    don't match.
 *
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t
FUNCTION(igraph_matrix, all_ge)(const TYPE(igraph_matrix) *lhs,
                                const TYPE(igraph_matrix) *rhs) {
    return lhs->ncol == rhs->ncol && lhs->nrow == rhs->nrow &&
           FUNCTION(igraph_vector, all_ge)(&lhs->data, &rhs->data);
}

#endif

#ifndef NOTORDERED

/**
 * \function igraph_matrix_maxdifference
 * \brief Maximum absolute difference between two matrices.
 *
 * Calculate the maximum absolute difference of two matrices. Both matrices
 * must be non-empty. If their dimensions differ then a warning is given and
 * the comparison is performed by vectors columnwise from both matrices.
 * The remaining elements in the larger vector are ignored.
 * \param m1 The first matrix.
 * \param m2 The second matrix.
 * \return The element with the largest absolute value in \c m1 - \c m2.
 *
 * Time complexity: O(mn), the elements in the smaller matrix.
 */

igraph_real_t FUNCTION(igraph_matrix, maxdifference)(const TYPE(igraph_matrix) *m1,
        const TYPE(igraph_matrix) *m2) {
    long int col1 = FUNCTION(igraph_matrix, ncol)(m1);
    long int col2 = FUNCTION(igraph_matrix, ncol)(m2);
    long int row1 = FUNCTION(igraph_matrix, nrow)(m1);
    long int row2 = FUNCTION(igraph_matrix, nrow)(m2);
    if (col1 != col2 || row1 != row2) {
        IGRAPH_WARNING("Comparing non-conformant matrices");
    }
    return FUNCTION(igraph_vector, maxdifference)(&m1->data, &m2->data);
}

#endif

/**
 * \function igraph_matrix_transpose
 * \brief Transpose a matrix.
 *
 * Calculate the transpose of a matrix. Note that the function
 * reallocates the memory used for the matrix.
 * \param m The input (and output) matrix.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix, transpose)(TYPE(igraph_matrix) *m) {
    long int nrow = m->nrow;
    long int ncol = m->ncol;
    if (nrow > 1 && ncol > 1) {
        TYPE(igraph_vector) newdata;
        long int i, size = nrow * ncol, mod = size - 1;
        FUNCTION(igraph_vector, init)(&newdata, size);
        IGRAPH_FINALLY(FUNCTION(igraph_vector, destroy), &newdata);
        for (i = 0; i < size; i++) {
            VECTOR(newdata)[i] = VECTOR(m->data)[ (i * nrow) % mod ];
        }
        VECTOR(newdata)[size - 1] = VECTOR(m->data)[size - 1];
        FUNCTION(igraph_vector, destroy)(&m->data);
        IGRAPH_FINALLY_CLEAN(1);
        m->data = newdata;
    }
    m->nrow = ncol;
    m->ncol = nrow;

    return 0;
}

/**
 * \function igraph_matrix_e
 * Extract an element from a matrix.
 *
 * Use this if you need a function for some reason and cannot use the
 * \ref MATRIX macro. Note that no range checking is performed.
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \return The element in the given row and column.
 *
 * Time complexity: O(1).
 */

BASE FUNCTION(igraph_matrix, e)(const TYPE(igraph_matrix) *m,
                                long int row, long int col) {
    return MATRIX(*m, row, col);
}

/**
 * \function igraph_matrix_e_ptr
 * Pointer to an element of a matrix.
 *
 * The function returns a pointer to an element. No range checking is
 * performed.
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \return Pointer to the element in the given row and column.
 *
 * Time complexity: O(1).
 */

BASE* FUNCTION(igraph_matrix, e_ptr)(const TYPE(igraph_matrix) *m,
                                     long int row, long int col) {
    return &MATRIX(*m, row, col);
}

/**
 * \function igraph_matrix_set
 * Set an element.
 *
 * Set an element of a matrix. No range checking is performed.
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \param value The new value of the element.
 *
 * Time complexity: O(1).
 */

void FUNCTION(igraph_matrix, set)(TYPE(igraph_matrix)* m, long int row, long int col,
                                  BASE value) {
    MATRIX(*m, row, col) = value;
}

/**
 * \function igraph_matrix_fill
 * Fill with an element.
 *
 * Set the matrix to a constant matrix.
 * \param m The input matrix.
 * \param e The element to set.
 *
 * Time complexity: O(mn), the number of elements.
 */

void FUNCTION(igraph_matrix, fill)(TYPE(igraph_matrix) *m, BASE e) {
    FUNCTION(igraph_vector, fill)(&m->data, e);
}

/**
 * \function igraph_matrix_update
 * Update from another matrix.
 *
 * This function replicates \p from in the matrix \p to.
 * Note that \p to must be already initialized.
 * \param to The result matrix.
 * \param from The matrix to replicate; it is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, update)(TYPE(igraph_matrix) *to,
                                    const TYPE(igraph_matrix) *from) {

    IGRAPH_CHECK(FUNCTION(igraph_matrix, resize)(to, from->nrow, from->ncol));
    FUNCTION(igraph_vector, update)(&to->data, &from->data);
    return 0;
}

/**
 * \function igraph_matrix_rbind
 * Combine two matrices rowwise.
 *
 * This function places the rows of \p from below the rows of \c to
 * and stores the result in \p to. The number of columns in the two
 * matrices must match.
 * \param to The upper matrix; the result is also stored here.
 * \param from The lower matrix. It is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the newly created
 * matrix.
 */

int FUNCTION(igraph_matrix, rbind)(TYPE(igraph_matrix) *to,
                                   const TYPE(igraph_matrix) *from) {
    long int tocols = to->ncol, fromcols = from->ncol;
    long int torows = to->nrow, fromrows = from->nrow;
    long int offset, c, r, index, offset2;
    if (tocols != fromcols) {
        IGRAPH_ERROR("Cannot do rbind, number of columns do not match", IGRAPH_EINVAL);
    }

    IGRAPH_CHECK(FUNCTION(igraph_vector, resize)(&to->data,
                 tocols * (fromrows + torows)));
    to->nrow += fromrows;

    offset = (tocols - 1) * fromrows;
    index = tocols * torows - 1;
    for (c = tocols - 1; c > 0; c--) {
        for (r = 0; r < torows; r++, index--) {
            VECTOR(to->data)[index + offset] = VECTOR(to->data)[index];
        }
        offset -= fromrows;
    }

    offset = torows; offset2 = 0;
    for (c = 0; c < tocols; c++) {
        memcpy(VECTOR(to->data) + offset, VECTOR(from->data) + offset2,
               sizeof(BASE) * (size_t) fromrows);
        offset += fromrows + torows;
        offset2 += fromrows;
    }
    return 0;
}

/**
 * \function igraph_matrix_cbind
 * Combine matrices columnwise.
 *
 * This function places the columns of \p from on the right of \p to,
 * and stores the result in \p to.
 * \param to The left matrix; the result is stored here too.
 * \param from The right matrix. It is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements on the new matrix.
 */

int FUNCTION(igraph_matrix, cbind)(TYPE(igraph_matrix) *to,
                                   const TYPE(igraph_matrix) *from) {

    long int tocols = to->ncol, fromcols = from->ncol;
    long int torows = to->nrow, fromrows = from->nrow;
    if (torows != fromrows) {
        IGRAPH_ERROR("Cannot do rbind, number of rows do not match", IGRAPH_EINVAL);
    }
    IGRAPH_CHECK(FUNCTION(igraph_matrix, resize)(to, torows, tocols + fromcols));
    FUNCTION(igraph_vector, copy_to)(&from->data, VECTOR(to->data) + tocols * torows);
    return 0;
}

/**
 * \function igraph_matrix_swap
 * Swap two matrices.
 *
 * The contents of the two matrices will be swapped. They must have the
 * same dimensions.
 * \param m1 The first matrix.
 * \param m2 The second matrix.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the matrices.
 */

int FUNCTION(igraph_matrix, swap)(TYPE(igraph_matrix) *m1, TYPE(igraph_matrix) *m2) {
    if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
        IGRAPH_ERROR("Cannot swap non-conformant matrices", IGRAPH_EINVAL);
    }
    return FUNCTION(igraph_vector, swap)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_get_row
 * Extract a row.
 *
 * Extract a row from a matrix and return it as a vector.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; it will be resized if
 *   needed.
 * \param index The index of the row to select.
 * \return Error code.
 *
 * Time complexity: O(n), the number of columns in the matrix.
 */

int FUNCTION(igraph_matrix, get_row)(const TYPE(igraph_matrix) *m,
                                     TYPE(igraph_vector) *res, long int index) {
    long int rows = m->nrow, cols = m->ncol;
    long int i, j;

    if (index >= rows) {
        IGRAPH_ERROR("Index out of range for selecting matrix row", IGRAPH_EINVAL);
    }
    IGRAPH_CHECK(FUNCTION(igraph_vector, resize)(res, cols));

    for (i = index, j = 0; j < cols; i += rows, j++) {
        VECTOR(*res)[j] = VECTOR(m->data)[i];
    }
    return 0;
}

/**
 * \function igraph_matrix_set_row
 * Set a row from a vector.
 *
 * Sets the elements of a row with the given vector. This has the effect of
 * setting row \c index to have the elements in the vector \c v. The length of
 * the vector and the number of columns in the matrix must match,
 * otherwise an error is triggered.
 * \param m The input matrix.
 * \param v The vector containing the new elements of the row.
 * \param index Index of the row to set.
 * \return Error code.
 *
 * Time complexity: O(n), the number of columns in the matrix.
 */

int FUNCTION(igraph_matrix, set_row)(TYPE(igraph_matrix) *m,
                                     const TYPE(igraph_vector) *v, long int index) {
    long int rows = m->nrow, cols = m->ncol;
    long int i, j;

    if (index >= rows) {
        IGRAPH_ERROR("Index out of range for selecting matrix row", IGRAPH_EINVAL);
    }
    if (FUNCTION(igraph_vector, size)(v) != cols) {
        IGRAPH_ERROR("Cannot set matrix row, invalid vector length", IGRAPH_EINVAL);
    }
    for (i = index, j = 0; j < cols; i += rows, j++) {
        VECTOR(m->data)[i] = VECTOR(*v)[j];
    }
    return 0;
}

/**
 * \function igraph_matrix_set_col
 * Set a column from a vector.
 *
 * Sets the elements of a column with the given vector. In effect, column
 * \c index will be set with elements from the vector \c v. The length of
 * the vector and the number of rows in the matrix must match,
 * otherwise an error is triggered.
 * \param m The input matrix.
 * \param v The vector containing the new elements of the column.
 * \param index Index of the column to set.
 * \return Error code.
 *
 * Time complexity: O(m), the number of rows in the matrix.
 */

int FUNCTION(igraph_matrix, set_col)(TYPE(igraph_matrix) *m,
                                     const TYPE(igraph_vector) *v, long int index) {
    long int rows = m->nrow, cols = m->ncol;
    long int i, j;

    if (index >= cols) {
        IGRAPH_ERROR("Index out of range for setting matrix column", IGRAPH_EINVAL);
    }
    if (FUNCTION(igraph_vector, size)(v) != rows) {
        IGRAPH_ERROR("Cannot set matrix column, invalid vector length", IGRAPH_EINVAL);
    }
    for (i = index * rows, j = 0; j < rows; i++, j++) {
        VECTOR(m->data)[i] = VECTOR(*v)[j];
    }
    return 0;
}

/**
 * \function igraph_matrix_swap_rows
 * Swap two rows.
 *
 * Swap two rows in the matrix.
 * \param m The input matrix.
 * \param i The index of the first row.
 * \param j The index of the second row.
 * \return Error code.
 *
 * Time complexity: O(n), the number of columns.
 */

int FUNCTION(igraph_matrix, swap_rows)(TYPE(igraph_matrix) *m,
                                       long int i, long int j) {
    long int ncol = m->ncol, nrow = m->nrow;
    long int n = nrow * ncol;
    long int index1, index2;
    if (i >= nrow || j >= nrow) {
        IGRAPH_ERROR("Cannot swap rows, index out of range", IGRAPH_EINVAL);
    }
    if (i == j) {
        return 0;
    }
    for (index1 = i, index2 = j; index1 < n; index1 += nrow, index2 += nrow) {
        BASE tmp;
        tmp = VECTOR(m->data)[index1];
        VECTOR(m->data)[index1] = VECTOR(m->data)[index2];
        VECTOR(m->data)[index2] = tmp;
    }
    return 0;
}

/**
 * \function igraph_matrix_swap_cols
 * Swap two columns.
 *
 * Swap two columns in the matrix.
 * \param m The input matrix.
 * \param i The index of the first column.
 * \param j The index of the second column.
 * \return Error code.
 *
 * Time complexity: O(m), the number of rows.
 */

int FUNCTION(igraph_matrix, swap_cols)(TYPE(igraph_matrix) *m,
                                       long int i, long int j) {
    long int ncol = m->ncol, nrow = m->nrow;
    long int k, index1, index2;
    if (i >= ncol || j >= ncol) {
        IGRAPH_ERROR("Cannot swap columns, index out of range", IGRAPH_EINVAL);
    }
    if (i == j) {
        return 0;
    }
    for (index1 = i * nrow, index2 = j * nrow, k = 0; k < nrow; k++, index1++, index2++) {
        BASE tmp = VECTOR(m->data)[index1];
        VECTOR(m->data)[index1] = VECTOR(m->data)[index2];
        VECTOR(m->data)[index2] = tmp;
    }
    return 0;
}

/**
 * \function igraph_matrix_add_constant
 * Add a constant to every element.
 *
 * \param m The input matrix.
 * \param plud The constant to add.
 *
 * Time complexity: O(mn), the number of elements.
 */

void FUNCTION(igraph_matrix, add_constant)(TYPE(igraph_matrix) *m, BASE plus) {
    FUNCTION(igraph_vector, add_constant)(&m->data, plus);
}

/**
 * \function igraph_matrix_add
 * Add two matrices.
 *
 * Add \p m2 to \p m1, and store the result in \p m1. The dimensions of the
 * matrices must match.
 * \param m1 The first matrix; the result will be stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, add)(TYPE(igraph_matrix) *m1,
                                 const TYPE(igraph_matrix) *m2) {
    if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
        IGRAPH_ERROR("Cannot add non-conformant matrices", IGRAPH_EINVAL);
    }
    return FUNCTION(igraph_vector, add)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_sub
 * Difference of two matrices.
 *
 * Subtract \p m2 from \p m1 and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The first matrix; the result is stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, sub)(TYPE(igraph_matrix) *m1,
                                 const TYPE(igraph_matrix) *m2) {
    if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
        IGRAPH_ERROR("Cannot subtract non-conformant matrices", IGRAPH_EINVAL);
    }
    return FUNCTION(igraph_vector, sub)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_mul_elements
 * Elementwise multiplication.
 *
 * Multiply \p m1 by \p m2 elementwise and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The first matrix; the result is stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, mul_elements)(TYPE(igraph_matrix) *m1,
        const TYPE(igraph_matrix) *m2) {
    if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
        IGRAPH_ERROR("Cannot multiply non-conformant matrices", IGRAPH_EINVAL);
    }
    return FUNCTION(igraph_vector, mul)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_div_elements
 * Elementwise division.
 *
 * Divide \p m1 by \p m2 elementwise and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The dividend. The result is store here.
 * \param m2 The divisor. It is left unchanged.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, div_elements)(TYPE(igraph_matrix) *m1,
        const TYPE(igraph_matrix) *m2) {
    if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
        IGRAPH_ERROR("Cannot divide non-conformant matrices", IGRAPH_EINVAL);
    }
    return FUNCTION(igraph_vector, div)(&m1->data, &m2->data);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_min
 * Minimum element.
 *
 * Returns the smallest element of a non-empty matrix.
 * \param m The input matrix.
 * \return The smallest element.
 *
 * Time complexity: O(mn), the number of elements.
 */

igraph_real_t FUNCTION(igraph_matrix, min)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, min)(&m->data);
}

/**
 * \function igraph_matrix_which_min
 * Indices of the minimum.
 *
 * Gives the indices of the (first) smallest element in a non-empty
 * matrix.
 * \param m The matrix.
 * \param i Pointer to a <type>long int</type>. The row index of the
 *   minimum is stored here.
 * \param j Pointer to a <type>long int</type>. The column index of
 *   the minimum is stored here.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, which_min)(const TYPE(igraph_matrix) *m,
                                       long int *i, long int *j) {
    long int vmin = FUNCTION(igraph_vector, which_min)(&m->data);
    *i = vmin % m->nrow;
    *j = vmin / m->nrow;
    return 0;
}

/**
 * \function igraph_matrix_which_max
 * Indices of the maximum.
 *
 * Gives the indices of the (first) largest element in a non-empty
 * matrix.
 * \param m The matrix.
 * \param i Pointer to a <type>long int</type>. The row index of the
 *   maximum is stored here.
 * \param j Pointer to a <type>long int</type>. The column index of
 *   the maximum is stored here.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, which_max)(const TYPE(igraph_matrix) *m,
                                       long int *i, long int *j) {
    long int vmax = FUNCTION(igraph_vector, which_max)(&m->data);
    *i = vmax % m->nrow;
    *j = vmax / m->nrow;
    return 0;
}

/**
 * \function igraph_matrix_minmax
 * Minimum and maximum
 *
 * The maximum and minimum elements of a non-empty matrix.
 * \param m The input matrix.
 * \param min Pointer to a base type. The minimum is stored here.
 * \param max Pointer to a base type. The maximum is stored here.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, minmax)(const TYPE(igraph_matrix) *m,
                                    BASE *min, BASE *max) {
    return FUNCTION(igraph_vector, minmax)(&m->data, min, max);
}

/**
 * \function igraph_matrix_which_minmax
 * Indices of the minimum and maximum
 *
 * Find the positions of the smallest and largest elements of a
 * non-empty matrix.
 * \param m The input matrix.
 * \param imin Pointer to a <type>long int</type>, the row index of
 *   the minimum is stored here.
 * \param jmin Pointer to a <type>long int</type>, the column index of
 *   the minimum is stored here.
 * \param imax Pointer to a <type>long int</type>, the row index of
 *   the maximum is stored here.
 * \param jmax Pointer to a <type>long int</type>, the column index of
 *   the maximum is stored here.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix, which_minmax)(const TYPE(igraph_matrix) *m,
        long int *imin, long int *jmin,
        long int *imax, long int *jmax) {
    long int vmin, vmax;
    FUNCTION(igraph_vector, which_minmax)(&m->data, &vmin, &vmax);
    *imin = vmin % m->nrow;
    *jmin = vmin / m->nrow;
    *imax = vmax % m->nrow;
    *jmax = vmax / m->nrow;
    return 0;
}

#endif

/**
 * \function igraph_matrix_isnull
 * Check for a null matrix.
 *
 * Checks whether all elements are zero.
 * \param m The input matrix.
 * \return Boolean, \c TRUE is \p m contains only zeros and \c FALSE
 *   otherwise.
 *
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix, isnull)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, isnull)(&m->data);
}

/**
 * \function igraph_matrix_empty
 * Check for an empty matrix.
 *
 * It is possible to have a matrix with zero rows or zero columns, or
 * even both. This functions checks for these.
 * \param m The input matrix.
 * \return Boolean, \c TRUE if the matrix contains zero elements, and
 *    \c FALSE otherwise.
 *
 * Time complexity: O(1).
 */

igraph_bool_t FUNCTION(igraph_matrix, empty)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, empty)(&m->data);
}

/**
 * \function igraph_matrix_is_symmetric
 * Check for symmetric matrix.
 *
 * A non-square matrix is not symmetric by definition.
 * \param m The input matrix.
 * \return Boolean, \c TRUE if the matrix is square and symmetric, \c
 *    FALSE otherwise.
 *
 * Time complexity: O(mn), the number of elements. O(1) for non-square
 * matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix, is_symmetric)(const TYPE(igraph_matrix) *m) {

    long int n = m->nrow;
    long int r, c;
    if (m->ncol != n) {
        return 0;
    }
    for (r = 1; r < n; r++) {
        for (c = 0; c < r; c++) {
            BASE a1 = MATRIX(*m, r, c);
            BASE a2 = MATRIX(*m, c, r);
#ifdef EQ
            if (!EQ(a1, a2)) {
                return 0;
            }
#else
            if (a1 != a2) {
                return 0;
            }
#endif
        }
    }
    return 1;
}

/**
 * \function igraph_matrix_prod
 * Product of the elements.
 *
 * Note this function can result in overflow easily, even for not too
 * big matrices.
 * \param m The input matrix.
 * \return The product of the elements.
 *
 * Time complexity: O(mn), the number of elements.
 */

BASE FUNCTION(igraph_matrix, prod)(const TYPE(igraph_matrix) *m) {
    return FUNCTION(igraph_vector, prod)(&m->data);
}

/**
 * \function igraph_matrix_rowsum
 * Rowwise sum.
 *
 * Calculate the sum of the elements in each row.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; the result is stored
 *   here. It will be resized if necessary.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix, rowsum)(const TYPE(igraph_matrix) *m,
                                    TYPE(igraph_vector) *res) {
    long int nrow = m->nrow, ncol = m->ncol;
    long int r, c;
    BASE sum;
    IGRAPH_CHECK(FUNCTION(igraph_vector, resize)(res, nrow));
    for (r = 0; r < nrow; r++) {
        sum = ZERO;
        for (c = 0; c < ncol; c++) {
#ifdef SUM
            SUM(sum, sum, MATRIX(*m, r, c));
#else
            sum += MATRIX(*m, r, c);
#endif
        }
        VECTOR(*res)[r] = sum;
    }
    return 0;
}

/**
 * \function igraph_matrix_colsum
 * Columnwise sum.
 *
 * Calculate the sum of the elements in each column.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; the result is stored
 *   here. It will be resized if necessary.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix, colsum)(const TYPE(igraph_matrix) *m,
                                    TYPE(igraph_vector) *res) {
    long int nrow = m->nrow, ncol = m->ncol;
    long int r, c;
    BASE sum;
    IGRAPH_CHECK(FUNCTION(igraph_vector, resize)(res, ncol));
    for (c = 0; c < ncol; c++) {
        sum = ZERO;
        for (r = 0; r < nrow; r++) {
#ifdef SUM
            SUM(sum, sum, MATRIX(*m, r, c));
#else
            sum += MATRIX(*m, r, c);
#endif
        }
        VECTOR(*res)[c] = sum;
    }
    return 0;
}

/**
 * \function igraph_matrix_contains
 * Search for an element.
 *
 * Search for the given element in the matrix.
 * \param m The input matrix.
 * \param e The element to search for.
 * \return Boolean, \c TRUE if the matrix contains \p e, \c FALSE
 * otherwise.
 *
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix, contains)(const TYPE(igraph_matrix) *m,
        BASE e) {
    return FUNCTION(igraph_vector, contains)(&m->data, e);
}

/**
 * \function igraph_matrix_search
 * Search from a given position.
 *
 * Search for an element in a matrix and start the search from the
 * given position. The search is performed columnwise.
 * \param m The input matrix.
 * \param from The position to search from, the positions are
 *    enumerated columnwise.
 * \param what The element to search for.
 * \param pos Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to the position of its first appearance.
 * \param row Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to its row index.
 * \param col Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to its column index.
 * \return Boolean, \c TRUE if the element is found, \c FALSE
 *    otherwise.
 *
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix, search)(const TYPE(igraph_matrix) *m,
        long int from, BASE what,
        long int *pos,
        long int *row, long int *col) {
    igraph_bool_t find = FUNCTION(igraph_vector, search)(&m->data, from, what, pos);
    if (find) {
        *row = *pos % m->nrow;
        *col = *pos / m->nrow;
    }
    return find;
}

/**
 * \function igraph_matrix_remove_row
 * Remove a row.
 *
 * A row is removed from the matrix.
 * \param m The input matrix.
 * \param row The index of the row to remove.
 * \return Error code.
 *
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix, remove_row)(TYPE(igraph_matrix) *m, long int row) {

    long int c, r, index = row + 1, leap = 1, n = m->nrow * m->ncol;
    if (row >= m->nrow) {
        IGRAPH_ERROR("Cannot remove row, index out of range", IGRAPH_EINVAL);
    }

    for (c = 0; c < m->ncol; c++) {
        for (r = 0; r < m->nrow - 1 && index < n; r++) {
            VECTOR(m->data)[index - leap] = VECTOR(m->data)[index];
            index++;
        }
        leap++;
        index++;
    }
    m->nrow--;
    FUNCTION(igraph_vector, resize)(&m->data, m->nrow * m->ncol);
    return 0;
}

/**
 * \function igraph_matrix_select_cols
 * \brief Select some columns of a matrix.
 *
 * This function selects some columns of a matrix and returns them in a
 * new matrix. The result matrix should be initialized before calling
 * the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param cols Vector; it contains the column indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 *
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix, select_cols)(const TYPE(igraph_matrix) *m,
        TYPE(igraph_matrix) *res,
        const igraph_vector_t *cols) {
    long int ncols = igraph_vector_size(cols);
    long int nrows = m->nrow;
    long int i, j;

    IGRAPH_CHECK(FUNCTION(igraph_matrix, resize)(res, nrows, ncols));
    for (i = 0; i < nrows; i++) {
        for (j = 0; j < ncols; j++) {
            MATRIX(*res, i, j) = MATRIX(*m, i, (long int)VECTOR(*cols)[j]);
        }
    }
    return 0;
}

#ifdef OUT_FORMAT

#ifndef USING_R
int FUNCTION(igraph_matrix, print)(const TYPE(igraph_matrix) *m) {

    long int nr = FUNCTION(igraph_matrix, nrow)(m);
    long int nc = FUNCTION(igraph_matrix, ncol)(m);
    long int i, j;
    for (i = 0; i < nr; i++) {
        for (j = 0; j < nc; j++) {
            if (j != 0) {
                putchar(' ');
            }
            printf(OUT_FORMAT, MATRIX(*m, i, j));
        }
        printf("\n");
    }

    return 0;
}

int FUNCTION(igraph_matrix, printf)(const TYPE(igraph_matrix) *m,
                                    const char *format) {
    long int nr = FUNCTION(igraph_matrix, nrow)(m);
    long int nc = FUNCTION(igraph_matrix, ncol)(m);
    long int i, j;
    for (i = 0; i < nr; i++) {
        for (j = 0; j < nc; j++) {
            if (j != 0) {
                putchar(' ');
            }
            printf(format, MATRIX(*m, i, j));
        }
        printf("\n");
    }

    return 0;
}

#endif

int FUNCTION(igraph_matrix, fprint)(const TYPE(igraph_matrix) *m,
                                    FILE *file) {

    long int nr = FUNCTION(igraph_matrix, nrow)(m);
    long int nc = FUNCTION(igraph_matrix, ncol)(m);
    long int i, j;
    for (i = 0; i < nr; i++) {
        for (j = 0; j < nc; j++) {
            if (j != 0) {
                fputc(' ', file);
            }
            fprintf(file, OUT_FORMAT, MATRIX(*m, i, j));
        }
        fprintf(file, "\n");
    }

    return 0;
}

#endif