complexmatrix2d.pmc

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

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

INTVAL __PLA_Type_ComplexMatrix2D;

static void
get_complex_value_from_pmc(PARROT_INTERP, PMC * value, FLOATVAL * real, FLOATVAL * imag)
{
    switch (value->vtable->base_type) {
        case enum_class_String:
            value = Parrot_pmc_new_init(interp, enum_class_Complex, value);
        case enum_class_Complex:
            *real = VTABLE_get_number_keyed_int(interp, value, 0);
            *imag = VTABLE_get_number_keyed_int(interp, value, 1);
            break;
        case enum_class_Float:
            *real = VTABLE_get_number(interp, value);
            *imag = 0.0;
            break;
        case enum_class_Integer:
            {
                const INTVAL _r = VTABLE_get_integer(interp, value);
                *real = (FLOATVAL)_r;
            }
            *imag = 0.0;
            break;
        default:
            /* TODO: We should support HLL-mapped types here. */
            Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": cannot set unknown PMC type");
    }
}

/* 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) {
                R_ITEM_XY_COLMAJOR(new_s, rows_size, cols_size, i, j) =
                    R_ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, i, j);
                I_ITEM_XY_COLMAJOR(new_s, rows_size, cols_size, i, j) =
                    I_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);
    }
}

static PMC *
convert_to_ComplexMatrix2D(PARROT_INTERP, PMC * p, INTVAL forcecopy)
{
    if (p->vtable->base_type == __PLA_Type_ComplexMatrix2D) {
        normalize_lazy_transpose(interp, p);
        if (forcecopy)
            return VTABLE_clone(interp, p);
        return p;
    }
    else if (!VTABLE_does(interp, p, CONST_STRING(interp, "matrix")))
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": cannot convert unknown PMC type");
    else {
        PMC * const meth = VTABLE_find_method(interp, p, CONST_STRING(interp, "convert_to_complex_matrix"));
        PMC * dest = PMCNULL;
        Parrot_ext_call(interp, meth, "Pi->P", p, &dest);
        return dest;
    }
}


/* 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 y 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++) {
            R_ITEM_XY_ROWMAJOR(new_s, new_rows, new_cols, i, j) =
                R_ITEM_XY_ROWMAJOR(old_s, old_rows, old_cols, i, j);
            I_ITEM_XY_ROWMAJOR(new_s, new_rows, new_cols, i, j) =
                I_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);
}

static void
init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, INTVAL cols_size, PMC * values) {
    DECLATTRS(self, attrs);
    INTVAL num = 0;
    const INTVAL init_elems = VTABLE_elements(interp, values);
    const INTVAL total_elems = rows_size * cols_size;
    resize_matrix(interp, self, rows_size - 1, cols_size - 1);

    for (; num < init_elems && num < total_elems ; num++) {
        PMC * const value = VTABLE_get_pmc_keyed_int(interp, values, num);
        VTABLE_set_pmc_keyed_int(interp, self, num, value);
    }
}

static PMC *
get_complex_pmc_at_xy(PARROT_INTERP, PMC *self, INTVAL rows, INTVAL cols)
{
    DECLATTRS(self, attrs);
    const INTVAL rows_size = attrs->rows;
    const INTVAL cols_size = attrs->cols;
    PMC * newcomplex;
    const INTVAL flags = attrs->flags;
    FLOATVAL real, imag;
    FLOATVAL * const s = attrs->storage;

    if (rows >= rows_size || cols >= cols_size || cols < 0 || rows < 0)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": indices out of bounds");
    real = R_ITEM_XY(s, flags, rows_size, cols_size, rows, cols);
    imag = I_ITEM_XY(s, flags, rows_size, cols_size, rows, cols);
    newcomplex = Parrot_pmc_new(interp, enum_class_Complex);
    VTABLE_set_number_keyed_int(interp, newcomplex, 0, real);
    VTABLE_set_number_keyed_int(interp, newcomplex, 1, imag);
    return newcomplex;
}

static void
set_complex_pmc_at_xy(PARROT_INTERP, PMC * self, PMC * value, INTVAL row, INTVAL col)
{
    DECLATTRS(self, attrs);
    INTVAL rows_size = attrs->rows;
    INTVAL cols_size = attrs->cols;
    FLOATVAL real, imag;
    const INTVAL flags = attrs->flags;
    FLOATVAL * s;

    if (row >= rows_size || col >= cols_size) {
        resize_matrix(interp, self, row, col);
        rows_size = attrs->rows;
        cols_size = attrs->cols;
    }
    s = attrs->storage;
    get_complex_value_from_pmc(interp, value, &real, &imag);
    R_ITEM_XY(s, flags, rows_size, cols_size, row, col) = real;
    I_ITEM_XY(s, flags, rows_size, cols_size, row, col) = imag;
}

static void
set_scalar_at_xy(PARROT_INTERP, PMC * self, FLOATVAL value, INTVAL row, INTVAL col)
{
    DECLATTRS(self, attrs);
    INTVAL rows_size = attrs->rows;
    INTVAL cols_size = attrs->cols;
    const INTVAL flags = attrs->flags;
    FLOATVAL * s;

    if (row >= rows_size || col >= cols_size) {
        resize_matrix(interp, self, row, col);
        rows_size = attrs->rows;
        cols_size = attrs->cols;
    }
    s = attrs->storage;
    R_ITEM_XY(s, flags, rows_size, cols_size, row, col) = value;
    I_ITEM_XY(s, flags, rows_size, cols_size, row, col) = 0;
}

/* Wrapper to call the dgemm function from BLAS with PMC arguments. Assumes
   A, B, and C are all ComplexMatrix2D. */
static void
call_gemm(PARROT_INTERP, FLOATVAL alpha_r, FLOATVAL alpha_i, PMC * A, PMC *B,
    FLOATVAL beta_r, FLOATVAL beta_i, 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;
    FLOATVAL alpha_p[2];
    FLOATVAL beta_p[2];
    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");
    alpha_p[0] = alpha_r;
    alpha_p[1] = alpha_i;
    beta_p[0] = beta_r;
    beta_p[1] = beta_i;
    cblas_zgemm(CblasRowMajor,
        IS_TRANSPOSED_BLAS(attrs_a->flags),
        IS_TRANSPOSED_BLAS(attrs_b->flags),
        M,
        N,
        K,
        alpha_p,
        attrs_a->storage,
        M,
        attrs_b->storage,
        N,
        beta_p,
        attrs_c->storage,
        M
    );
}

/* 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++) {
            R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) += v;
        }
    }
}

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

    if (sub) {
        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) -= v_r;
                I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) -= v_i;
            }
        }
    }
    else {
        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) += v_r;
                I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) += v_i;
            }
        }
    }
}

/* 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++) {
            R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) *= v;
            I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) *= v;
        }
    }
}

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

    for (i = 0; i < rows_size; i++) {
        for (j = 0; j < cols_size; j++) {
            const FLOATVAL s_r = R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j);
            const FLOATVAL s_i = I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j);
            const FLOATVAL new_r = (s_r * v_r) - (s_i * v_i);
            const FLOATVAL new_i = (s_r * v_i) + (s_i * v_r);
            R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) = new_r;
            I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, i, j) = new_i;
        }
    }
}

/* 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++) {
                R_ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) -=
                    R_ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
                I_ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) -=
                    R_ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
            }
        }
    }
    else {
        for (i = 0; i < rows; i++) {
            for (j = 0; j < cols; j++) {
                R_ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) +=
                    R_ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
                I_ITEM_XY_ROWMAJOR(s_a, rows, cols, i, j) +=
                    R_ITEM_XY_ROWMAJOR(s_b, rows, cols, i, j);
            }
        }
    }
}


static void
intkey_to_coords(PARROT_INTERP, const INTVAL rows, const INTVAL cols,
    const INTVAL key, INTVAL * row, INTVAL * col)
{
    if (rows == 0) {
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
            PLATYPENAME ": Index out of bounds.");
    } else {
        const INTVAL r = key / rows;
        const INTVAL c = key % rows;
        if (key < 0 || r > rows || c > cols) {
            Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Index out of bounds.");
        }
        *row = r;
        *col = c;
    }
}

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

    void class_init() {
        __PLA_Type_ComplexMatrix2D = entry;
    }

/*

=head1 VTABLEs

=over 4

=item* init

=cut

*/

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

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

/*

=item* get_number_keyed

=item* get_integer_keyed

=item* get_string_keyed

=item* get_pmc_keyed

=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;
        FLOATVAL real, imag;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        if (rows < 0 || rows >= rows_size || cols >= cols_size || cols < 0)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": indices out of bounds");
        real = R_ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols);
        imag = I_ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols);
        return sqrt(real * real + imag + imag);
    }

    VTABLE FLOATVAL get_number_keyed_int(INTVAL key) {
        DECLATTRS(SELF, attrs);
        const INTVAL totalsize = attrs->rows * attrs->cols;
        FLOATVAL real, imag;
        if (key >= totalsize || key < 0) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Index out of bounds.");
        }
        real = attrs->storage[key * 2];
        imag = attrs->storage[key * 2 + 1];
        return sqrt(real * real + imag * imag);
    }

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

    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(PMC * key) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        INTVAL rows, cols;
        const INTVAL flags = attrs->flags;
        FLOATVAL real, imag;
        FLOATVAL * const s = attrs->storage;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        if (rows < 0 || rows >= rows_size || cols >= cols_size || cols < 0)
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": indices out of bounds");
        real = R_ITEM_XY(s, flags, rows_size, cols_size, rows, cols);
        imag = I_ITEM_XY(s, flags, rows_size, cols_size, rows, cols);
        return Parrot_sprintf_c(INTERP, "%vg%+vgi", real, imag);
    }

    VTABLE STRING * get_string_keyed_int(INTVAL key) {
        DECLATTRS(SELF, attrs);
        FLOATVAL real, imag;
        const INTVAL totalsize = attrs->rows * attrs->cols;
        if (key >= totalsize || totalsize < 0) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Index out of bounds.");
        }
        real = attrs->storage[key * 2];
        imag = attrs->storage[key * 2 + 1];
        return Parrot_sprintf_c(INTERP, "%vg%+vgi", real, imag);
    }

    VTABLE PMC * get_pmc_keyed(PMC * key) {
        INTVAL rows, cols;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        return get_complex_pmc_at_xy(INTERP, SELF, rows, cols);
    }

    VTABLE PMC * get_pmc_keyed_int(INTVAL key) {
        DECLATTRS(SELF, attrs);
        INTVAL row, col;
        intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col);
        return get_complex_pmc_at_xy(INTERP, SELF, row, col);
    }

/*

=item* set_pmc_keyed

=item* set_string_keyed

=cut

*/

    VTABLE void set_pmc_keyed(PMC * key, PMC * value) {
        INTVAL rows, cols;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        set_complex_pmc_at_xy(INTERP, SELF, value, rows, cols);
    }

    VTABLE void set_pmc_keyed_int(INTVAL key, PMC * value) {
        DECLATTRS(SELF, attrs);
        INTVAL row, col;
        intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col);
        set_complex_pmc_at_xy(INTERP, SELF, value, row, col);
    }

    VTABLE void set_number_keyed(PMC * key, FLOATVAL value) {
        INTVAL rows, cols;
        GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols);
        set_scalar_at_xy(INTERP, SELF, value, rows, cols);
    }

    VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) {
        DECLATTRS(SELF, attrs);
        INTVAL row, col;
        intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col);
        set_scalar_at_xy(INTERP, SELF, value, row, col);
    }

    VTABLE void set_string_keyed(PMC * key, STRING * value) {
        PMC * const item = Parrot_pmc_new(INTERP, enum_class_Complex);
        VTABLE_set_string_native(INTERP, item, value);
        VTABLE_set_pmc_keyed(INTERP, SELF, key, item);
    }

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

    VTABLE void set_integer_keyed_int(INTVAL key, INTVAL value) {
        VTABLE_set_number_keyed_int(INTERP, SELF, key, (FLOATVAL)value);
    }

    VTABLE void set_string_keyed_int(INTVAL key, STRING * value) {
        DECLATTRS(SELF, attrs);
        INTVAL row, col;
        PMC * const item = Parrot_pmc_new(INTERP, enum_class_Complex);
        VTABLE_set_string_native(INTERP, item, value);
        VTABLE_set_pmc_keyed_int(INTERP, SELF, key, item);
    }

/*

=item* get_string

=cut

*/

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

        if (IS_TRANSPOSED(attrs->flags)) {
            for (cols = 0; cols < cols_size; ++cols) {
                for (rows = 0; rows < rows_size; ++rows) {
                    const FLOATVAL real = R_ITEM_XY_COLMAJOR(s, rows_size, cols_size, rows, cols);
                    const FLOATVAL imag = I_ITEM_XY_COLMAJOR(s, rows_size, cols_size, rows, cols);
                    STRING * const item = Parrot_sprintf_c(INTERP, "\t%vg%+vgi", real, imag);
                    VTABLE_push_string(INTERP, builder, item);
                }
                VTABLE_push_string(INTERP, builder, newline);
            }
        }
        else {
            for (rows = 0; rows < rows_size; ++rows) {
                for (cols = 0; cols < cols_size; ++cols) {
                    const FLOATVAL real = R_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, rows, cols);
                    const FLOATVAL imag = I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, rows, cols);
                    STRING * const item = Parrot_sprintf_c(INTERP, "\t%vg%+vgi", real, imag);
                    VTABLE_push_string(INTERP, builder, item);
                }
                VTABLE_push_string(INTERP, builder, newline);
            }
        }
        return VTABLE_get_string(INTERP, builder);
    }

    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;
    }

    VTABLE PMC * clone() {
        PMC * const c = Parrot_pmc_new(INTERP, SELF->vtable->base_type);
        DECLATTRS(SELF, old_atts);
        DECLATTRS(c, new_atts);
        INTVAL rows, cols;
        INTVAL const rows_size = old_atts->rows;
        INTVAL const cols_size = old_atts->cols;
        INTVAL const newsize = rows_size * cols_size;
        FLOATVAL * const old_s = old_atts->storage;
        FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize);
        for (rows = 0; rows < rows_size; ++rows) {
            for (cols = 0; cols < cols_size; ++cols) {
                R_ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, rows, cols) =
                    R_ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, rows, cols);
                I_ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, rows, cols) =
                    I_ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, rows, cols);
            }
        }
        new_atts->storage = new_s;
        new_atts->flags = old_atts->flags;
        new_atts->rows = rows_size;
        new_atts->cols = cols_size;
        return c;
    }

    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 rows, cols;

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

            for (cols = 0; cols < self_cols; cols++) {
                for (rows = 0; rows < self_rows; rows++) {
                    const FLOATVAL self_r_value =
                        R_ITEM_XY(self_s, self_flags, self_rows, self_cols, rows, cols);
                    const FLOATVAL other_r_value =
                        R_ITEM_XY(other_s, other_flags, other_rows, other_cols, rows, cols);
                    const FLOATVAL self_i_value =
                        I_ITEM_XY(self_s, self_flags, self_rows, self_cols, rows, cols);
                    const FLOATVAL other_i_value =
                        I_ITEM_XY(other_s, other_flags, other_rows, other_cols, rows, cols);
                    if (self_r_value != other_r_value || self_i_value != other_i_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++) {
                FLOATVAL f = R_ITEM_XY(s, flags, rows, cols, i, j);
                VTABLE_push_float(INTERP, info, f);
                f = I_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++) {
                FLOATVAL f = VTABLE_shift_float(INTERP, info);
                R_ITEM_XY(s, flags, rows, cols, i, j) = f;
                f = VTABLE_shift_float(INTERP, info);
                I_ITEM_XY(s, flags, rows, cols, i, j) = f;
            }
        }
    }

/*

=item* add(ComplexMatrix2D)

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(ComplexMatrix2D *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_matrices(INTERP, dest, value, 0);
        return dest;
    }

    MULTI PMC *add(Complex *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_complex(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(ComplexMatrix2D)

=item* i_add(DEFAULT)

=item* i_add_int

=item* i_add_float

=cut

*/

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

    MULTI void i_add(Complex* value) {
        add_scalar_complex(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(ComplexMatrix2D)

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(ComplexMatrix2D *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_matrices(INTERP, dest, value, 1);
        return dest;
    }

    MULTI PMC *subtract(Complex* value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        add_scalar_complex(INTERP, SELF, 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(ComplexMatrix2D)

=item* i_subtract(DEFAULT)

=item* i_subtract_int

=item* i_subtract_float

=cut

*/

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

    MULTI void i_sub(Complex* value) {
        add_scalar_complex(INTERP, SELF, value, 1);
    }

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

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

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

/*

=item* multiply(ComplexMatrix2D)

=item* multiply(DEFAULT)

=cut

*/

    MULTI PMC *multiply(ComplexMatrix2D *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, 0.0, SELF, value, 0.0, 0.0, dest);
        return dest;
    }

    MULTI PMC *mutiply(Complex *value, PMC *dest) {
        dest = VTABLE_clone(INTERP, SELF);
        multiply_scalar_complex(INTERP, SELF, 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(ComplexMatrix2D* value) {
        PMC * const temp = VTABLE_clone(INTERP, SELF);
        call_gemm(INTERP, 1.0, 0.0, temp, value, 0.0, 0.0, SELF);
    }

    MULTI void i_multiply(Complex *value) {
        multiply_scalar_complex(INTERP, SELF, value);
    }

    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);
    }


/*

=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 fill()

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

=cut

*/

    METHOD fill(PMC * value,
                INTVAL rows_size :optional, INTVAL has_rows_size :opt_flag,
                INTVAL cols_size :optional, INTVAL has_cols_size :opt_flag
               ) {
        /* TODO: Value here is going to be a Complex PMC. Handle that */
        DECLATTRS(SELF, attrs);
        FLOATVAL * s  = attrs->storage;
        INTVAL const curr_rows_size = attrs->rows;
        INTVAL const curr_cols_size = attrs->cols;
        INTVAL rows, cols;

        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 (rows = 0; rows < rows_size; ++rows) {
            for (cols = 0; cols < cols_size; ++cols) {
                set_complex_pmc_at_xy(INTERP, SELF, value, rows, cols);
            }
        }
    }

/*

=item item_at()

=cut

*/

    METHOD item_at(INTVAL row, INTVAL col,
        PMC * value :optional, INTVAL has_value :opt_flag) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows = attrs->rows;
        const INTVAL cols = attrs->cols;
        if (row < 0 || col < 0 || row >= rows || col >= cols) {
            Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS,
                PLATYPENAME ": Index out of bounds.");
        }
        if (has_value) {
            set_complex_pmc_at_xy(INTERP, SELF, value, row, col);
            RETURN(PMC * value);
        }
        else {
            PMC * const p = get_complex_pmc_at_xy(INTERP, SELF, row, col);
            RETURN(PMC * p);
        }
    }

/*

=item transpose()

Transposes the matrix.

=cut

*/

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

        SWAP_XY(attrs);

        if (transposed)
            attrs->flags -= FLAG_TRANSPOSED;
        else
            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_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        const INTVAL newsize = rows_size * cols_size;
        FLOATVAL * new_s = ALLOCATE_STORAGE(newsize);
        FLOATVAL * old_s = attrs->storage;

        INTVAL i, j;
        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                R_ITEM_XY_ROWMAJOR(new_s, cols_size, rows_size, j, i) =
                    R_ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, i, j);
                I_ITEM_XY_ROWMAJOR(new_s, cols_size, rows_size, j, i) =
                    I_ITEM_XY_ROWMAJOR(old_s, rows_size, cols_size, i, j);
            }
        }
        attrs->storage = new_s;
        SWAP_XY(attrs);
        if (old_s)
            free(old_s);
    }

/*

=item conjugate

=cut

*/

    METHOD conjugate() {
        DECLATTRS(SELF, attrs);
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        const INTVAL newsize = rows_size * cols_size;
        FLOATVAL * const s = attrs->storage;

        INTVAL i, j;
        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                I_ITEM_XY_ROWMAJOR(s, cols_size, rows_size, j, i) =
                    -I_ITEM_XY_ROWMAJOR(s, rows_size, cols_size, j, i);
            }
        }
    }

/*

=item iterate_function_inplace()

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

=cut

*/

    METHOD iterate_function_inplace(PMC * func, PMC * args :slurpy) {
        DECLATTRS(SELF, attrs);
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        INTVAL i, j;

        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                PMC * const item_xy = get_complex_pmc_at_xy(INTERP, SELF, i, j);
                PMC * result = PMCNULL;
                Parrot_ext_call(INTERP, func, "PPIIPf->P", SELF, item_xy, i, j, args, &result);
                set_complex_pmc_at_xy(INTERP, SELF, result, i, j);
            }
        }
    }

    METHOD iterate_function_external(PMC * func, PMC * args :slurpy) {
        DECLATTRS(SELF, attrs);
        PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type);
        const INTVAL rows_size = attrs->rows;
        const INTVAL cols_size = attrs->cols;
        INTVAL i, j;

        resize_matrix(INTERP, new_matrix, rows_size - 1, cols_size - 1);

        for (i = 0; i < rows_size; i++) {
            for (j = 0; j < cols_size; j++) {
                PMC * const item_xy = get_complex_pmc_at_xy(INTERP, SELF, i, j);
                PMC * result = PMCNULL;
                Parrot_ext_call(INTERP, func, "PPIIPf->P", SELF, item_xy, i, j, args, &result);
                VTABLE_set_pmc_keyed_int(INTERP, new_matrix,
                    INDEX_XY_ROWMAJOR(rows_size, cols_size, i, j), result);
            }
        }
        RETURN(PMC * new_matrix);
    }

    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);
    }

    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);
            Parrot_ComplexMatrix2D_attributes * const new_attrs = PARROT_COMPLEXMATRIX2D(new_matrix);
            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++) {
                    R_ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, j, i) =
                        R_ITEM_XY_ROWMAJOR(s, rows, cols, j + rows_idx, i + cols_idx);
                    I_ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, j, i) =
                        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++) {
                R_ITEM_XY_ROWMAJOR(self_s, self_rows, self_cols, j + rows_idx, i + cols_idx) =
                    R_ITEM_XY_ROWMAJOR(blck_s, blck_rows, blck_cols, j, i);
                I_ITEM_XY_ROWMAJOR(self_s, self_rows, self_cols, j + rows_idx, i + cols_idx) =
                    I_ITEM_XY_ROWMAJOR(blck_s, blck_rows, blck_cols, j, i);
            }
        }
    }

/*

=item* gemm

Calculates the matrix equation:

    Z = aAB + bC

=cut

*/

    METHOD gemm(PMC *alpha, PMC * A, PMC *B, PMC *beta, PMC *C) {
        FLOATVAL alpha_r, alpha_i, beta_r, beta_i;
        A = convert_to_ComplexMatrix2D(interp, A, 0);
        B = convert_to_ComplexMatrix2D(interp, B, 0);
        C = convert_to_ComplexMatrix2D(interp, C, 1);
        get_complex_value_from_pmc(interp, alpha, &alpha_r, &alpha_i);
        get_complex_value_from_pmc(interp, beta, &beta_r, &beta_i);
        call_gemm(INTERP, alpha_r, alpha_i, A, B, beta_r, beta_i, C);
        RETURN(PMC* C);
    }

/*

=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++) {
            R_ITEM_XY(s, flags, rows, cols, destidx, i) +=
                R_ITEM_XY(s, flags, rows, cols, srcidx, i) * gain;
            I_ITEM_XY(s, flags, rows, cols, destidx, i) +=
                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++) {
            R_ITEM_XY(s, flags, rows, cols, idx, i) *= gain;
            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_r = R_ITEM_XY(s, flags, rows, cols, idx_b, i);
            const FLOATVAL t_i = I_ITEM_XY(s, flags, rows, cols, idx_b, i);
            R_ITEM_XY(s, flags, rows, cols, idx_b, i) =
                R_ITEM_XY(s, flags, rows, cols, idx_a, i);
            I_ITEM_XY(s, flags, rows, cols, idx_b, i) =
                I_ITEM_XY(s, flags, rows, cols, idx_a, i);
            R_ITEM_XY(s, flags, rows, cols, idx_a, i) = t_r;
            I_ITEM_XY(s, flags, rows, cols, idx_a, i) = t_i;
        }
    }

/*

=item convert_to_number_matrix

Get a NumMatrix2D from the current matrix

=item convert_to_complex_matrix

Get a ComplexMatrix2D from the current matrix

=item convert_to_pmc_matrix

Get a PMCMatrix2D from the current matrix

=cut

*/

    METHOD convert_to_number_matrix() {
        DECLATTRS(SELF, attrs);
        PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_NumMatrix2D);
        const INTVAL totalsize = attrs->rows * attrs->cols;
        PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize"));
        INTVAL i;
        Parrot_ext_call(INTERP, meth, "PiII->", d, attrs->rows, attrs->cols);
        normalize_lazy_transpose(INTERP, SELF);
        for (i = 0; i < totalsize; i++) {
            const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, i);
            VTABLE_set_number_keyed_int(INTERP, d, i, f);
        }
        RETURN(PMC * d);
    }

    METHOD convert_to_complex_matrix() {
        PMC * const d = VTABLE_clone(INTERP, SELF);
        normalize_lazy_transpose(INTERP, d);
        RETURN(PMC * d);
    }

    METHOD convert_to_pmc_matrix() {
        DECLATTRS(SELF, attrs);
        PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_PMCMatrix2D);
        const INTVAL totalsize = attrs->rows * attrs->cols;
        PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize"));
        INTVAL i;
        Parrot_ext_call(INTERP, meth, "PiII->", d, attrs->rows, attrs->cols);
        normalize_lazy_transpose(INTERP, SELF);
        for (i = 0; i < totalsize; i++) {
            PMC * const p = VTABLE_get_pmc_keyed_int(INTERP, SELF, i);
            VTABLE_set_pmc_keyed_int(INTERP, d, i, p);
        }
        RETURN(PMC * d);
    }


/*

=back

=end

*/
}