nummatrix2d.pmc

#include "pla_matrix_types.h"
#include <math.h>
#include <cblas.h>

#define ALLOCATE_STORAGE(s) (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL))
#define PLATYPENAME "NumMatrix2D"
#define DECLATTRS(p, a) Parrot_NumMatrix2D_attributes * const (a) = \
    (Parrot_NumMatrix2D_attributes *)((p)->data)

/* Relative error between two large numbers */
#define FLOATVAL_EPSILON_REL 0.0001

/* Absolute difference between two small numbers */
#define FLOATVAL_EPSILON_ABS 0.00001

/* Above we do relative error, below we use abs */
#define FLOATVAL_REL_THRESHOLD 10

/* Resize the matrix internal storage to be able to hold a point at position
   (x, y). The matrix grows but does not shrink. New spaces in the matrix
   are initialized to 0.0. Parameters x and y are the indices that are trying
   to be accessed, so we must resize the matrix to be able to accomodate those
   indices. Notice that the matrix type is zero-indexed, so the size is one
   plus the highest index that we need to access.

   This function will not shrink the matrix, only grow it. After the call,
   the matrix will be at least large enough to hold an item at the given
   index. To shrink the matrix, use the get_block() method, or that
   algorithm. */
static void
resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col)
{
    DECLATTRS(self, attrs);
    /* Store the old values */
    const INTVAL old_rows = attrs->rows;
    const INTVAL old_cols = attrs->cols;
    FLOATVAL * old_s = attrs->storage;

    /* rows and cols are indices, not sizes. Resize the matrix to accomodate
       this new point without shrinking by taking the max. */
    const INTVAL new_rows = INDEX_MAX(old_rows, row + 1);
    const INTVAL new_cols = INDEX_MAX(old_cols, col + 1);
    const INTVAL newsize = new_rows * new_cols;
    FLOATVAL * new_s = ALLOCATE_STORAGE(newsize);

    INTVAL i, j;
    for (i = 0; i < old_rows; i++) {
        for (j = 0; j < old_cols; j++) {
            ITEM_XY_ROWMAJOR(new_s, new_rows, new_cols, i, j) =
                ITEM_XY_ROWMAJOR(old_s, old_rows, old_cols, i, j);
        }
    }
    attrs->storage = new_s;
    attrs->rows = new_rows;
    attrs->cols = new_cols;
    if (old_s)
        mem_sys_free(old_s);
}

/* If the matrix is lazily transposed, actually transpose the physical memory
   layout. This is necessary for calculations, especially BLAS calculations,
   which aren't lazy-transpose-aware. */
static void
normalize_lazy_transpose(PARROT_INTERP, PMC * self)
{
    DECLATTRS(self, attrs);
    if (IS_TRANSPOSED(attrs->flags)) {
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        const INTVAL size = rows_size * cols_size;
        FLOATVAL * const new_s = ALLOCATE_STORAGE(size);
        FLOATVAL * const old_s = attrs->storage;
        INTVAL i, j;

        for (i = 0; i < rows_size; ++i) {
            for (j = 0; j < cols_size; ++j) {
                ITEM_XY_COLMAJOR(new_s, rows_size, cols_size, i, j) =
                    ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, i, j);
            }
        }
        attrs->storage = new_s;
        mem_sys_free(old_s);
        FLAG_CLEAR(attrs->flags, FLAG_TRANSPOSED);
    }
}

/* Initialize the PMC from an array, filling the matrix row-at-a-time */
static void
init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size,
    INTVAL cols_size, PMC * values)
{
    DECLATTRS(self, attrs);
    FLOATVAL * s;
    INTVAL self_rows, self_cols, i, j, num = 0;
    const INTVAL init_elems = VTABLE_elements(interp, values);

    resize_matrix(interp, self, rows_size - 1, cols_size - 1);
    self_rows = attrs->rows;
    self_cols = attrs->cols;
    s = attrs->storage;

    for (i = 0; i < rows_size; i++) {
        for (j = 0; j < cols_size; j++) {
            const FLOATVAL value = VTABLE_get_number_keyed_int(interp, values, num);
            num++;
            ITEM_XY_ROWMAJOR(s, self_rows, self_cols, i, j) = value;
            if (num >= init_elems)
                return;
        }
    }
}

/* Add a scalar to every element of the matrix */
static void
add_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v)
{
    DECLATTRS(self, attrs);
    const INTVAL rows_size = attrs->rows;
    const INTVAL cols_size = attrs->cols;
    FLOATVAL * const s = attrs->storage;
    INTVAL i, j;

    for (i = 0; i < rows_size; i++) {
        for (j = 0; j < cols_size; j++) {
            ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) += v;
        }
    }
}

/* Multiply every entry by a scalar */
static void
multiply_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v)
{
    DECLATTRS(self, attrs);
    const INTVAL rows_size = attrs->rows;
    const INTVAL cols_size = attrs->cols;
    FLOATVAL * const s = attrs->storage;
    INTVAL i, j;

    for (i = 0; i < rows_size; i++) {
        for (j = 0; j < cols_size; j++) {
            ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) *= v;
        }
    }
}

/* item-by-item addition or subtraction
   A = A + B */
static void
add_matrices(PARROT_INTERP, PMC * A, PMC * B, INTVAL sub)
{
    DECLATTRS(A, attrs_a);
    DECLATTRS(B, attrs_b);
    const INTVAL rows = attrs_a->rows;
    const INTVAL cols = attrs_a->cols;
    FLOATVAL * const s_a = attrs_a->storage;
    FLOATVAL * const s_b = attrs_b->storage;
    INTVAL i, j;
    if (attrs_b->rows != rows || attrs_b->cols != cols)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": sizes do not match");
    if (sub) {
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) -=
                    ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
            }
        }
    }
    else {
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) +=
                    ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
            }
        }
    }
}

/* Wrapper to call the dgemm function from BLAS with PMC arguments. Assumes
   A, B, and C are all NumMatrix2D. */
static void
call_gemm(PARROT_INTERP, FLOATVAL alpha, PMC * A, PMC *B, FLOATVAL beta, PMC *C)
{
    DECLATTRS(A, attrs_a);
    DECLATTRS(B, attrs_b);
    DECLATTRS(C, attrs_c);
    const INTVAL M = attrs_a->rows;
    const INTVAL N = attrs_b->cols;
    const INTVAL K = attrs_a->cols;
    if (attrs_c->rows != M)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": A, C indices do not match in gemm");
    if (attrs_c->cols != N)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": B, C indices do not match in gemm");
    if (attrs_b->rows != K)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": A, B indeces do not match in gemm");
    cblas_dgemm(CblasRowMajor,
        IS_TRANSPOSED_BLAS(attrs_a->flags),
        IS_TRANSPOSED_BLAS(attrs_b->flags),
        M,
        N,
        K,
        alpha,
        attrs_a->storage,
        M,
        attrs_b->storage,
        N,
        beta,
        attrs_c->storage,
        M
    );
}

/* Check if two floats are equalish. Above a certain threshold we want to do a
   relative error calculation because precision between subsequent large
   numbers can be much larger than standard epsilon values. Below the threshold
   we want to do an absolute calculation because as the values approach zero the
   error value quickly begins to diverge, even for very similar values. The
   values for absolute difference, relative error, and threshold need to be
   tuned. */
static INTVAL
floats_are_equal(FLOATVAL a, FLOATVAL b)
{
    FLOATVAL diff = fabs(a - b);
    if (a == b)
        return 1;
    if (diff > FLOATVAL_EPSILON_ABS)
        return 0;
    if (a > FLOATVAL_REL_THRESHOLD) {
        FLOATVAL err = diff / b;
        if (err > FLOATVAL_EPSILON_REL)
            return 0;
    }
    return 1;
}


pmclass NumMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix {
    ATTR FLOATVAL * storage;
    ATTR INTVAL rows;
    ATTR INTVAL cols;
    ATTR INTVAL flags;

/*

=head1 VTABLEs

=over 4

=item* init

Create a new NumMatrix2D

=item* destroy

Destrow the matrix and free it's allocated storage

=cut

*/

    VTABLE void init() {
        DECLATTRS(SELF, a);
        a->storage = NULL;
        a->rows = 0;
        a->cols = 0;
        a->flags = 0;
        PObj_custom_destroy_SET(SELF);
    }

    VTABLE void destroy() {
        DECLATTRS(SELF, a);
        FLOATVAL * const s = a->storage;
        if (s)
            mem_sys_free(s);
    }

/*

=item* get_number_keyed

Get the number at the location X, Y. The key must have two elements.

=item* get_integer_keyed

Get the integer at the location X, Y. The key must have two elements.

=item* get_string_keyed

Get a string representation of the number at the location X, Y. The key must
have two elements.

=item* get_pmc_keyed

Get a Float PMC of the number at the location X, Y. The key must have two
elements.

=cut

*/

    VTABLE FLOATVAL get_number_keyed(PMC * key) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        INTVAL rows, cols;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        if (rows >= rows_size || cols >= cols_size)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": indices out of bounds");
        return ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols);
    }

    VTABLE INTVAL get_integer_keyed(PMC * key) {
        const FLOATVAL f = VTABLE_get_number_keyed(INTERP, SELF, key);
        return (INTVAL)f;
    }

    VTABLE STRING * get_string_keyed(PMC * key) {
        const FLOATVAL f = VTABLE_get_number_keyed(INTERP, SELF, key);
        return Parrot_str_from_num(INTERP, f);
    }

    VTABLE PMC * get_pmc_keyed(PMC * key) {
        const FLOATVAL f = VTABLE_get_number_keyed(INTERP, SELF, key);
        PMC * const item = Parrot_pmc_new(INTERP, enum_class_Float);
        VTABLE_set_number_native(INTERP, item, f);
        return item;
    }

/*

=item* get_number_keyed_int

Treating the matrix memory block like a C array, get the Nth number. The
array is arranged in memory row by row.

=item* get_integer_keyed_int

Treating the matrix memory block like a C array, get the Nth integer. The
array is arranged in memory row by row.

=item* get_string_keyed_int

Treating the matrix memory block like a C array, get the string representation
of the Nth number. The array is arranged in memory row by row.

=item* get_pmc_keyed_int

Treating the matrix memory block like a C array, get a Float PMC for the
Nth number. The array is arranged in memory row by row.

=cut

*/

    VTABLE FLOATVAL get_number_keyed_int(INTVAL key) {
        DECLATTRS(SELF, attrs);
        const INTVAL total_size = attrs->rows * attrs->cols;
        if (key >= total_size) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Matrix dimensions must match in add.");
        }
        return attrs->storage[key];
    }

    VTABLE INTVAL get_integer_keyed_int(INTVAL key) {
        const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, key);
        return (INTVAL)f;
    }

    VTABLE STRING * get_string_keyed_int(INTVAL key) {
        const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, key);
        return Parrot_str_from_num(INTERP, f);
    }

    VTABLE PMC * get_pmc_keyed_int(INTVAL key) {
        const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, key);
        PMC * const item = Parrot_pmc_new(INTERP, enum_class_Float);
        VTABLE_set_number_native(INTERP, item, f);
        return item;
    }

/*

=item* set_number_keyed

Set the number at position (X, Y), growing the matrix if necessary and padding
empty spaces with 0.0. The key must have two elements.

=item* set_integer_keyed

Set the integer at position (X, Y), growing the matrix if necessary and
padding empty spaces with 0.0. The key must have two elements.

=item* set_pmc_keyed

Get the numeric value from the PMC argument and set it at position (X, Y),
growing the matrix if necessary and padding empty spaces with 0.0. The key
must have two elements.

=cut

*/

    VTABLE void set_number_keyed(PMC * key, FLOATVAL value) {
        DECLATTRS(SELF, attrs);
        INTVAL rows, cols, rows_size = attrs->rows, cols_size = attrs->cols;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        if (rows >= rows_size || cols >= cols_size) {
            resize_matrix(INTERP, SELF, rows, cols);
            rows_size = attrs->rows;
            cols_size = attrs->cols;
        }
        ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols) = value;
    }

    VTABLE void set_integer_keyed(PMC * key, INTVAL value) {
        VTABLE_set_number_keyed(INTERP, SELF, key, (FLOATVAL)value);
    }

    VTABLE void set_pmc_keyed(PMC * key, PMC * value) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        VTABLE_set_number_keyed(INTERP, SELF, key, v);
    }

/*

=item* get_string

Get a string representation of the matrix.

=cut

*/

    VTABLE STRING *get_string() {
        DECLATTRS(SELF, attrs);
        INTVAL i, j;
        PMC * const builder = Parrot_pmc_new(INTERP, enum_class_StringBuilder);
        STRING * const newline = Parrot_str_new(INTERP, "\n", 1);
        STRING * const tab = Parrot_str_new(INTERP, "\t", 1);
        FLOATVAL * const s = attrs->storage;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL flags = attrs->flags;

        for (i = 0; i < rows; ++i) {
            for (j = 0; j < cols; ++j) {
                const FLOATVAL f = ITEM_XY(s, flags, rows, cols, i, j);
                STRING * const item = Parrot_str_from_num(INTERP, f);
                VTABLE_push_string(INTERP, builder, tab);
                VTABLE_push_string(INTERP, builder, item);
            }
            VTABLE_push_string(INTERP, builder, newline);
        }

        return VTABLE_get_string(INTERP, builder);
    }

/*

=item* add(NumMatrix)

Add two matrices together, element-by-element.

=item* add(DEFAULT)

Take the number representation of the argument PMC and add it to every element
in the matrix.

=item* add_float

Add the float value to every element in the matrix.

=item* add_int

Add the integer value to every element in the matrix.

=cut

*/

    MULTI PMC *add(NumMatrix2D *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_matrices(INTERP, dest, value, 0);
        return dest;
    }

    MULTI PMC *add(DEFAULT *value, PMC *dest) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, v);
        return dest;
    }

    VTABLE PMC *add_float(FLOATVAL v, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, v);
        return dest;
    }

    VTABLE PMC *add_int(INTVAL v, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, (FLOATVAL)v);
        return dest;
    }

/*

=item* i_add(NumMatrix2D)

=item* i_add(DEFAULT)

=item* i_add_int

=item* i_add_float

=cut

*/

    MULTI void i_add(NumMatrix2D* value) {
        add_matrices(INTERP, SELF, value, 0);
    }

    MULTI void i_add(DEFAULT* value) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        add_scalar_float(INTERP, SELF, v);
    }

    void i_add_int(INTVAL v) {
        add_scalar_float(INTERP, SELF, (FLOATVAL)v);
    }

    void i_add_float(FLOATVAL v) {
        add_scalar_float(INTERP, SELF, v);
    }

/*

=item* subtract(NumMatrix)

Add two matrices together, element-by-element.

=item* subtract(DEFAULT)

Take the number representation of the argument PMC and add it to every element
in the matrix.

=item* subtract_float

Add the float value to every element in the matrix.

=item* subtract_int

Add the integer value to every element in the matrix.

=cut

*/

    MULTI PMC *subtract(NumMatrix2D *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_matrices(INTERP, dest, value, 1);
        return dest;
    }

    MULTI PMC *subtract(DEFAULT *value, PMC *dest) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, -v);
        return dest;
    }

    VTABLE PMC *subtract_float(FLOATVAL v, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, -v);
        return dest;
    }

    VTABLE PMC *subtract_int(INTVAL v, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_float(INTERP, dest, -(FLOATVAL)v);
        return dest;
    }

/*

=item* i_subtract(NumMatrix2D)

=item* i_subtract(DEFAULT)

=item* i_subtract_int

=item* i_subtract_float

=cut

*/

    MULTI void i_subtract(NumMatrix2D* value) {
        add_matrices(INTERP, SELF, value, 1);
    }

    MULTI void i_subtract(DEFAULT* value) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        add_scalar_float(INTERP, SELF, -v);
    }

    void i_subtract_int(INTVAL v) {
        add_scalar_float(INTERP, SELF, -(FLOATVAL)v);
    }

    void i_subtract_float(FLOATVAL v) {
        add_scalar_float(INTERP, SELF, -v);
    }

/*

=item* multiply(NumMatrix2D)

=item* multiply(DEFAULT)

=cut

*/

    MULTI PMC *multiply(NumMatrix2D *value, PMC *dest) {
        DECLATTRS(SELF, selfattr);
        DECLATTRS(value, valattr);
        const INTVAL new_rows = selfattr->rows;
        const INTVAL new_cols = valattr->cols;

        if (selfattr->cols != valattr->rows) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Dimensions do not match in matrix multiply");
        }

        dest = Parrot_pmc_new(INTERP, VTABLE_type(INTERP, SELF));
        resize_matrix(INTERP, dest, new_rows - 1, new_cols - 1);
        call_gemm(INTERP, 1.0, SELF, value, 0.0, dest);
        return dest;
    }

    MULTI PMC *multiply(DEFAULT *value, PMC *dest) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        dest = VTABLE_clone(INTERP, SELF);
        multiply_scalar_float(INTERP, dest, v);
        return dest;
    }

/*

=item* multiply_int

=item* multiply_float

=cut

*/

    VTABLE PMC *multiply_int(INTVAL v, PMC * dest) {
        dest = VTABLE_clone(INTERP, SELF);
        multiply_scalar_float(INTERP, dest, (FLOATVAL)v);
        return dest;
    }

    VTABLE PMC *multiply_float(FLOATVAL v, PMC * dest) {
        dest = VTABLE_clone(INTERP, SELF);
        multiply_scalar_float(INTERP, dest, v);
        return dest;
    }

/*

=item* i_multiply(NumMatrix2D)

=item* i_multiply(DEFAULT)

=cut

*/

    MULTI void i_multiply(NumMatrix2D* value) {
        PMC * const temp = VTABLE_clone(INTERP, SELF);
        call_gemm(INTERP, 1.0, temp, value, 0.0, SELF);
    }

    MULTI void i_multiply(DEFAULT* value) {
        const FLOATVAL v = VTABLE_get_number(INTERP, value);
        multiply_scalar_float(INTERP, SELF, v);
    }

/*

=item i_multiply_int

=item i_multiply_float

=cut

*/

    VTABLE void i_multiply_int(INTVAL v) {
        multiply_scalar_float(INTERP, SELF, (FLOATVAL)v);
    }

    VTABLE void i_multiply_float(FLOATVAL v) {
        multiply_scalar_float(INTERP, SELF, v);
    }

/*

=item* get_attr_str

=cut

*/

    VTABLE PMC * get_attr_str(STRING * idx) {
        DECLATTRS(SELF, attrs);
        if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "rows"))) {
            PMC * const rows = Parrot_pmc_new(INTERP, enum_class_Integer);
            VTABLE_set_integer_native(INTERP, rows, attrs->rows);
            return rows;
        }
        else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "cols"))) {
            PMC * const cols = Parrot_pmc_new(INTERP, enum_class_Integer);
            VTABLE_set_integer_native(INTERP, cols, attrs->cols);
            return cols;
        }
        else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "size"))) {
            PMC * const size = Parrot_pmc_new(INTERP, enum_class_Integer);
            VTABLE_set_integer_native(INTERP, size, attrs->cols * attrs->rows);
            return size;
        }
        return PMCNULL;
    }

/*

=item* clone

=cut

*/

    VTABLE PMC * clone() {
        PMC * const c = Parrot_pmc_new(INTERP, SELF->vtable->base_type);
        DECLATTRS(SELF, old_atts);
        DECLATTRS(c, new_atts);
        INTVAL const newsize = old_atts->rows * old_atts->cols;
        FLOATVAL * const old_s = old_atts->storage;
        FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize);
        memcpy(new_s, old_s, newsize * sizeof(FLOATVAL));
        memcpy(new_atts, old_atts, sizeof(Parrot_NumMatrix2D_attributes));
        new_atts->storage = new_s;
        return c;
    }

/*

=item* is_equal

=cut

*/

    VTABLE INTVAL is_equal(PMC * other) {
        if (other->vtable->base_type == SELF->vtable->base_type) {
            DECLATTRS(SELF, self_attrs);
            DECLATTRS(other, other_attrs);
            const INTVAL self_rows = self_attrs->rows;
            const INTVAL self_cols = self_attrs->cols;
            const INTVAL self_flags = self_attrs->flags;
            const INTVAL other_rows = other_attrs->rows;
            const INTVAL other_cols = other_attrs->cols;
            const INTVAL other_flags = other_attrs->flags;
            FLOATVAL * const self_s = self_attrs->storage;
            FLOATVAL * const other_s = other_attrs->storage;
            INTVAL i, j;

            if (self_rows != other_rows || self_cols != other_cols)
                return 0;

            for (i = 0; i < self_rows; i++) {
                for (j = 0; j < self_cols; j++) {
                    const FLOATVAL self_value =
                        ITEM_XY(self_s, self_flags, self_rows, self_cols, i, j);
                    const FLOATVAL other_value =
                        ITEM_XY(other_s, other_flags, other_rows, other_cols, i, j);
                    if (!floats_are_equal(self_value, other_value))
                        return 0;
                }
            }
            return 1;
        }
        return 0;
    }

/*

=item* freeze

=item* thaw

=cut

*/

    VTABLE void freeze(PMC *info) {
        DECLATTRS(SELF, attrs);
        INTVAL const rows = attrs->rows;
        INTVAL const cols = attrs->cols;
        INTVAL const flags = attrs->flags;
        INTVAL i, j;
        FLOATVAL * const s = attrs->storage;
        VTABLE_push_integer(INTERP, info, rows);
        VTABLE_push_integer(INTERP, info, cols);
        VTABLE_push_integer(INTERP, info, flags);
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                const FLOATVAL f = ITEM_XY(s, flags, rows, cols, i, j);
                VTABLE_push_float(INTERP, info, f);
            }
        }
    }

    VTABLE void thaw(PMC *info) {
        DECLATTRS(SELF, attrs);
        INTVAL const rows = VTABLE_shift_integer(INTERP, info);
        INTVAL const cols = VTABLE_shift_integer(INTERP, info);
        INTVAL const flags = VTABLE_shift_integer(INTERP, info);
        INTVAL i, j;
        FLOATVAL * s;
        attrs->rows = 0;
        attrs->cols = 0;
        attrs->storage = NULL;
        attrs->flags = 0;
        resize_matrix(INTERP, SELF, rows - 1, cols - 1);
        s = attrs->storage;
        attrs->flags = flags;
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                const FLOATVAL f = VTABLE_shift_float(INTERP, info);
                ITEM_XY(s, flags, rows, cols, i, j) = f;
            }
        }
    }

/*

=back

=head1 METHODS

=over 4

=item resize()

=cut

*/

    METHOD resize(INTVAL new_rows, INTVAL new_cols) {
        resize_matrix(INTERP, SELF, new_rows - 1, new_cols - 1);
    }

/*

=item get_proxy()

=cut

*/

    METHOD get_proxy(PMC * key) {
        PMC * const proxy = Parrot_pmc_new_init(INTERP, __PLA_MatrixProxy_typeid, SELF);
        VTABLE_assign_pmc(INTERP, proxy, key);
        RETURN(PMC * proxy);
    }

/*

=item fill()

Fill the matrix with a single value. if sizes are provided, fill to those
sizes, growing the matrix if needed.

=cut

*/

    METHOD fill(FLOATVAL value,
                INTVAL rows_size :optional, INTVAL has_rows_size :opt_flag,
                INTVAL cols_size :optional, INTVAL has_cols_size :opt_flag
               ) {
        DECLATTRS(SELF, attrs);
        FLOATVAL * s  = attrs->storage;
        INTVAL const curr_rows_size = attrs->rows;
        INTVAL const curr_cols_size = attrs->cols;
        INTVAL i, j;

        if (!has_rows_size)
            rows_size = curr_rows_size;
        if (!has_cols_size)
            cols_size = curr_cols_size;

        if ((has_rows_size && rows_size > curr_rows_size) ||
            (has_cols_size && cols_size > curr_cols_size)) {
            resize_matrix(INTERP, SELF, rows_size - 1, cols_size - 1);
            s = attrs->storage;
        }

        for (i = 0; i < rows_size; ++i) {
            for (j = 0; j < cols_size; ++j) {
                ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) = value;
            }
        }
    }

/*

=item item_at()

=cut

*/

    METHOD item_at(INTVAL row, INTVAL col,
        FLOATVAL value :optional,
        INTVAL has_value :opt_flag) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        if (row >= rows || col >= cols || row < 0 || col < 0) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": indices out of bounds in item_at");
        }
        if (has_value)
            ITEM_XY(attrs->storage, attrs->flags, rows, cols, row, col) = value;
        else
            value = ITEM_XY(attrs->storage, attrs->flags, rows, cols, row, col);
        RETURN(FLOATVAL value);
    }

/*

=item transpose()

Transposes the matrix lazily.

=cut

*/

    METHOD transpose() {
        DECLATTRS(SELF, attrs);
        INTVAL transposed = IS_TRANSPOSED(attrs->flags);

        SWAP_XY(attrs);

        if (transposed)
            FLAG_CLEAR(attrs->flags, FLAG_TRANSPOSED);
        else
            FLAG_SET(attrs->flags, FLAG_TRANSPOSED);
        transposed = !transposed;
        RETURN(INTVAL transposed);
    }

/*

=item mem_transpose()

Transposes the actual data storage of the matrix. More expensive up-front
than the transpose() method.

=cut

*/

    METHOD mem_transpose() {
        DECLATTRS(SELF, attrs);
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL newsize = rows * cols;
        FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize);
        FLOATVAL * const old_s = attrs->storage;

        INTVAL i, j;
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                ITEM_XY_ROWMAJOR(new_s, cols, rows, j, i) =
                    ITEM_XY_ROWMAJOR(old_s, rows, cols, i, j);
            }
        }
        attrs->storage = new_s;
        SWAP_XY(attrs);
        if (old_s)
            free(old_s);
    }

/*

=item iterate_function_inplace()

Calls a function for every element in the matrix, replacing the current
value with the return value of the called function.

=item iterate_function_external()

Calls a function for every element in the matrix, adding the result of each
computation to a new matrix. Return the new matrix of results.

=cut

*/

    METHOD iterate_function_inplace(PMC * func, PMC * args :slurpy) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL newsize = rows * cols;
        const INTVAL flags = attrs->flags;
        FLOATVAL * const old_s = attrs->storage;
        FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize);
        INTVAL i, j;

        if (newsize == 0 || old_s == NULL)
            RETURN();

        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                const FLOATVAL value = ITEM_XY(old_s, flags, rows, cols, i, j);
                FLOATVAL result = 0.0;
                Parrot_ext_call(INTERP, func, "PNIIPf->N", SELF, value, i, j, args, &result);
                ITEM_XY(new_s, flags, rows, cols, i, j) = result;
            }
        }
        attrs->storage = new_s;
        if (old_s)
            free(old_s);
    }

    METHOD iterate_function_external(PMC * func, PMC * args :slurpy) {
        DECLATTRS(SELF, attrs);
        PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type);
        Parrot_NumMatrix2D_attributes * new_attrs;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL newsize = rows * cols;
        const INTVAL flags = attrs->flags;
        FLOATVAL * const self_s = attrs->storage;
        FLOATVAL * new_s;
        INTVAL i, j;

        if (newsize == 0 || self_s == NULL)
            RETURN(PMC * new_matrix);

        resize_matrix(INTERP, new_matrix, rows - 1, cols - 1);
        new_attrs = PARROT_NUMMATRIX2D(new_matrix);
        new_s = new_attrs->storage;
        new_attrs->flags = flags;

        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                FLOATVAL value = ITEM_XY(self_s, flags, rows, cols, i, j);
                FLOATVAL result = 0.0;
                Parrot_ext_call(INTERP, func, "PNIIPf->N", SELF, value, i, j, args, &result);
                ITEM_XY(new_s, flags, rows, cols, i, j) = result;
            }
        }
        RETURN(PMC * new_matrix);
    }

/*

=item initialize_from_array()

Initialize matrix values from a linear array, filling each row with data
in order.

=item initialize_from_args()

Initialize matrix values from an array of function arguments, filling each row
with data in order.

=cut

*/

    METHOD initialize_from_array(INTVAL rows_size, INTVAL cols_size, PMC *values) {
        init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values);
    }

    METHOD initialize_from_args(INTVAL rows_size, INTVAL cols_size, PMC *values :slurpy) {
        init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values);
    }

/*

=item get_block

=item set_block

=cut

*/

    METHOD get_block(INTVAL rows_idx, INTVAL cols_idx, INTVAL rows_size, INTVAL cols_size) {
        DECLATTRS(SELF, attrs);
        FLOATVAL * const s = attrs->storage;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        INTVAL i, j;

        if ((rows < rows_idx + rows_size) || (cols < cols_idx + cols_size))
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Block boundaries exceed matrix size");
        if (rows_idx < 0 || cols_idx < 0)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Can not get block at negative location");
        if (rows_size < 0 || cols_size < 0)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Can not get block with negative size");
        else {
            PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type);
            DECLATTRS(new_matrix, new_attrs);
            FLOATVAL * new_s;
            resize_matrix(INTERP, new_matrix, rows_size - 1, cols_size - 1);
            new_s = new_attrs->storage;
            for (i = 0; i < cols_size; i++) {
                for (j = 0; j < rows_size; j++) {
                    ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, j, i) =
                        ITEM_XY_ROWMAJOR(s, rows, cols, j + rows_idx, i + cols_idx);
                }
            }
            RETURN(PMC * new_matrix);
        }
    }

    METHOD set_block(INTVAL rows_idx, INTVAL cols_idx, PMC * blck) {
        DECLATTRS(SELF, self_attrs);
        DECLATTRS(blck, blck_attrs);
        FLOATVAL * self_s = self_attrs->storage;
        FLOATVAL * const blck_s = blck_attrs->storage;
        INTVAL self_rows = self_attrs->rows;
        INTVAL self_cols = self_attrs->cols;
        const INTVAL blck_rows = blck_attrs->rows;
        const INTVAL blck_cols = blck_attrs->cols;
        INTVAL i, j;

        if (rows_idx < 0 || cols_idx < 0)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Can not set_block with negative indices");

        if (!VTABLE_does(INTERP, blck, CONST_STRING(interp, "matrix")))
            /* TODO: In this case, we should be able to do a .fill over the block
                     dimensions with VTABLE_get_number on the scalar value */
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Can not set_block on a non-matrix value");

        if (blck_rows + rows_idx > self_rows || blck_cols + cols_idx > self_cols) {
            resize_matrix(INTERP, SELF, blck_rows + rows_idx - 1, blck_cols + cols_idx - 1);
            self_s = self_attrs->storage;
            self_rows = self_attrs->rows;
            self_cols = self_attrs->cols;
        }
        for (i = 0; i < blck_cols; i++) {
            for (j = 0; j < blck_rows; j++) {
                ITEM_XY_ROWMAJOR(self_s, self_rows, self_cols, j + rows_idx, i + cols_idx) =
                    ITEM_XY_ROWMAJOR(blck_s, blck_rows, blck_cols, j, i);
            }
        }
    }

/*

=item* gemm

Calculates the matrix equation:

    Z = aAB + bC

=cut

*/

    METHOD gemm(FLOATVAL alpha, PMC * A, PMC *B, FLOATVAL beta, PMC *C) {
        PMC * const c_out = VTABLE_clone(INTERP, C);
        if (A->vtable->base_type != B->vtable->base_type &&
            C->vtable->base_type != A->vtable->base_type &&
            A->vtable->base_type != SELF->vtable->base_type)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": gemm only accepts NumMatrix2D arguments");
        call_gemm(INTERP, alpha, A, B, beta, c_out);
        RETURN(PMC* c_out);
    }

/*

=item row_combine(srcidx, destidx, gain)

add a multiple of the source row to the destination row.

=item row_scale(idx, gain)

Multiply all elements in the row by a gain factor.

=item row_swap(idx_a, idx_b)

Swap two rows

=cut

*/

    METHOD row_combine(INTVAL srcidx, INTVAL destidx, FLOATVAL gain) {
        DECLATTRS(SELF, attrs);
        FLOATVAL * const s = attrs->storage;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL flags = attrs->flags;
        INTVAL i;
        if (srcidx < 0 || srcidx >= rows || destidx < 0 || destidx >= rows)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Row index out of bounds");
        for (i = 0; i < cols; i++) {
            ITEM_XY(s, flags, rows, cols, destidx, i) +=
                ITEM_XY(s, flags, rows, cols, srcidx, i) * gain;
        }
    }

    METHOD row_scale(INTVAL idx, FLOATVAL gain) {
        DECLATTRS(SELF, attrs);
        FLOATVAL * const s = attrs->storage;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL flags = attrs->flags;
        INTVAL i;
        if (idx < 0 || idx >= rows)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Row index out of bounds");
        for (i = 0; i < cols; i++) {
            ITEM_XY(s, flags, rows, cols, idx, i) *= gain;
        }
    }

    METHOD row_swap(INTVAL idx_a, INTVAL idx_b) {
        DECLATTRS(SELF, attrs);
        FLOATVAL * const s = attrs->storage;
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        const INTVAL flags = attrs->flags;
        INTVAL i;
        if (idx_a < 0 || idx_a >= rows || idx_b < 0 || idx_b >= rows)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Row index out of bounds");
        for (i = 0; i < cols; i++) {
            const FLOATVAL t = ITEM_XY(s, flags, rows, cols, idx_b, i);
            ITEM_XY(s, flags, rows, cols, idx_b, i) =
                ITEM_XY(s, flags, rows, cols, idx_a, i);
            ITEM_XY(s, flags, rows, cols, idx_a, i) = t;
        }
    }

/*

=back

=end

*/
}