matblock.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <inttypes.h>

#include "jwsc/base/error.h"
#include "jwsc/math/complex.h"
#include "jwsc/math/blas.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matblock/matblock.h"

#ifdef JWSC_HAVE_DMALLOC
#include <dmalloc.h>
#endif


void create_matblock_u8
(
    Matblock_u8** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    uint32_t     mat;
    uint32_t     row;
    uint8_t*     elts;
    uint8_t**    row_ptrs;
    Matblock_u8* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_u8)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(uint8_t));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(uint8_t*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(uint8_t**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_u32
(
    Matblock_u32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    uint32_t      mat;
    uint32_t      row;
    uint32_t*     elts;
    uint32_t**    row_ptrs;
    Matblock_u32* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_u32)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(uint32_t));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(uint32_t*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(uint32_t**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_i32
(
    Matblock_i32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    uint32_t      mat;
    uint32_t      row;
    int32_t*      elts;
    int32_t**     row_ptrs;
    Matblock_i32* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_i32)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(int32_t));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(int32_t*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(int32_t**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_i64
(
    Matblock_i64** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    uint32_t      mat;
    uint32_t      row;
    int64_t*      elts;
    int64_t**     row_ptrs;
    Matblock_i64* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_i64)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(int64_t));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(int64_t*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(int64_t**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_f
(
    Matblock_f** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    uint32_t    mat;
    uint32_t    row;
    float*      elts;
    float**     row_ptrs;
    Matblock_f* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_f)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(float));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(float*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(float**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_d
(
    Matblock_d** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    uint32_t    mat;
    uint32_t    row;
    double*     elts;
    double**    row_ptrs;
    Matblock_d* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_d)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(double));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(double*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(double**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_cf
(
    Matblock_cf** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    uint32_t     mat;
    uint32_t     row;
    Complex_f*   elts;
    Complex_f**  row_ptrs;
    Matblock_cf* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_cf)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(Complex_f));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(Complex_f*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(Complex_f**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_matblock_cd
(
    Matblock_cd** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    uint32_t     mat;
    uint32_t     row;
    Complex_d*   elts;
    Complex_d**  row_ptrs;
    Matblock_cd* m;

    assert(num_mats > 0 && num_rows > 0 && num_cols > 0);

    if (*m_out == NULL)
    {
        assert(*m_out = malloc(sizeof(Matblock_cd)));
        (*m_out)->num_mats = 0;
        (*m_out)->num_rows = 0;
        (*m_out)->num_cols = 0;
        (*m_out)->elts     = NULL;
    }
    m = *m_out;

    if (m->num_mats != num_mats || 
        m->num_rows != num_rows || 
        m->num_cols != num_cols)
    {
        elts = (m->num_mats > 0 && m->num_rows > 0) ? **(m->elts) : NULL;
        elts = realloc(elts, num_mats*num_rows*num_cols*sizeof(Complex_d));
        assert(elts);

        row_ptrs = (m->num_mats > 0) ? *(m->elts) : NULL;
        row_ptrs = realloc(row_ptrs, num_mats*num_rows*sizeof(Complex_d*));
        assert(row_ptrs);

        m->elts = realloc(m->elts, num_mats*sizeof(Complex_d**));
        assert(m->elts);

        for (mat = 0; mat < num_mats; mat++)
        {
            m->elts[ mat ] = &(row_ptrs[ mat*num_rows ]);

            for (row = 0; row < num_rows; row++)
            {
                m->elts[ mat ][ row ] = &(elts[ (row+mat*num_rows)*num_cols ]);
            }
        }
        m->num_mats = num_mats;
        m->num_rows = num_rows;
        m->num_cols = num_cols;
    }
}

void create_init_matblock_u8
(
    Matblock_u8** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols,
    uint8_t       val
)
{
    uint8_t  num_elts;
    uint8_t  elt;
    uint8_t* m_elts;

    create_matblock_u8(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_u32
(
    Matblock_u32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    uint32_t       val
)
{
    uint32_t  num_elts;
    uint32_t  elt;
    uint32_t* m_elts;

    create_matblock_u32(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_i32
(
    Matblock_i32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    int32_t        val
)
{
    uint32_t num_elts;
    uint32_t elt;
    int32_t* m_elts;

    create_matblock_i32(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_i64
(
    Matblock_i64** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    int64_t        val
)
{
    uint32_t num_elts;
    uint32_t elt;
    int64_t* m_elts;

    create_matblock_i64(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_f
(
    Matblock_f** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    float        val
)
{
    uint32_t num_elts;
    uint32_t elt;
    float*   m_elts;

    create_matblock_f(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_d
(
    Matblock_d** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    double       val
)
{
    uint32_t num_elts;
    uint32_t elt;
    double*  m_elts;

    create_matblock_d(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_cf
(
    Matblock_cf** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols,
    Complex_f     val
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_f* m_elts;

    create_matblock_cf(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_init_matblock_cd
(
    Matblock_cd** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols,
    Complex_d     val
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_d* m_elts;

    create_matblock_cd(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        m_elts[ elt ] = val;
    }
}

void create_zero_matblock_u8
(
    Matblock_u8** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    create_init_matblock_u8(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_u32
(
    Matblock_u32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    create_init_matblock_u32(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_i32
(
    Matblock_i32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    create_init_matblock_i32(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_i64
(
    Matblock_i64** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols
)
{
    create_init_matblock_i64(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_f
(
    Matblock_f** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    create_init_matblock_f(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_d
(
    Matblock_d** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    create_init_matblock_d(m_out, num_mats, num_rows, num_cols, 0);
}

void create_zero_matblock_cf
(
    Matblock_cf** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    Complex_f z = {0};
    create_init_matblock_cf(m_out, num_mats, num_rows, num_cols, z);
}

void create_zero_matblock_cd
(
    Matblock_cd** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols
)
{
    Complex_d z = {0};
    create_init_matblock_cd(m_out, num_mats, num_rows, num_cols, z);
}

void create_random_matblock_u32
(
    Matblock_u32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    uint32_t       min,
    uint32_t       max
)
{
    uint32_t  num_elts;
    uint32_t  elt;
    uint32_t* m_elts;
    double    r;

    create_matblock_u32(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ] = min + (uint32_t)floor((max - min)*r + 0.5);
    }
}

void create_random_matblock_i32
(
    Matblock_i32** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    int32_t        min,
    int32_t        max
)
{
    uint32_t num_elts;
    uint32_t elt;
    int32_t* m_elts;
    double   r;

    create_matblock_i32(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ] = min + (int32_t)floor((max - min)*r + 0.5);
    }
}

void create_random_matblock_i64
(
    Matblock_i64** m_out, 
    uint32_t       num_mats, 
    uint32_t       num_rows, 
    uint32_t       num_cols,
    int64_t        min,
    int64_t        max
)
{
    uint32_t num_elts;
    uint32_t elt;
    int64_t* m_elts;
    double   r;

    create_matblock_i64(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ] = min + (int64_t)floor((max - min)*r + 0.5);
    }
}

void create_random_matblock_f
(
    Matblock_f** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    float        min,
    float        max
)
{
    uint32_t num_elts;
    uint32_t elt;
    float*   m_elts;
    float    r;

    create_matblock_f(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (float)rand() / (float)RAND_MAX;
        m_elts[ elt ] = min + (max - min)*r;
    }
}

void create_random_matblock_d
(
    Matblock_d** m_out, 
    uint32_t     num_mats, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    double       min,
    double       max
)
{
    uint32_t num_elts;
    uint32_t elt;
    double*  m_elts;
    double   r;

    create_matblock_d(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ] = min + (max - min)*r;
    }
}

void create_random_matblock_cf
(
    Matblock_cf** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols,
    Complex_f     min,
    Complex_f     max
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_f* m_elts;
    float      r;

    create_matblock_cf(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (float)rand() / (float)RAND_MAX;
        m_elts[ elt ].r = min.r + (max.r - min.r)*r;
        r = (float)rand() / (float)RAND_MAX;
        m_elts[ elt ].i = min.i + (max.i - min.i)*r;
    }
}

void create_random_matblock_cd
(
    Matblock_cd** m_out, 
    uint32_t      num_mats, 
    uint32_t      num_rows, 
    uint32_t      num_cols,
    Complex_d     min,
    Complex_d     max
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_d* m_elts;
    double     r;

    create_matblock_cd(m_out, num_mats, num_rows, num_cols);
    m_elts = **((*m_out)->elts);
    num_elts = num_mats*num_rows*num_cols;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ].r = min.r + (max.r - min.r)*r;
        r = (double)rand() / (double)RAND_MAX;
        m_elts[ elt ].i = min.i + (max.i - min.i)*r;
    }
}

Error* copy_matblock_u32(Matblock_u32** m_out, const Matblock_u32* m_in)
{
    uint32_t  num_elts;
    uint32_t  elt;
    uint32_t* m_out_elts;
    uint32_t* m_in_elts;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_u32(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }

    return NULL;
}

Error* copy_matblock_i32(Matblock_i32** m_out, const Matblock_i32* m_in)
{
    uint32_t  num_elts;
    uint32_t  elt;
    int32_t*  m_out_elts;
    int32_t*  m_in_elts;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_i32(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }

    return NULL;
}

Error* copy_matblock_i64(Matblock_i64** m_out, const Matblock_i64* m_in)
{
    uint32_t  num_elts;
    uint32_t  elt;
    int64_t*  m_out_elts;
    int64_t*  m_in_elts;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_i64(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }

    return NULL;
}

Error* copy_matblock_f(Matblock_f** m_out, const Matblock_f* m_in)
{
    uint32_t num_elts;
    float*   m_out_elts;
    float*   m_in_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_f(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_scopy(num_elts, m_in_elts, 1, m_out_elts, 1);
#else
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }
#endif

    return NULL;
}

Error* copy_matblock_d(Matblock_d** m_out, const Matblock_d* m_in)
{
    uint32_t num_elts;
    double*  m_out_elts;
    double*  m_in_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_d(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_dcopy(num_elts, m_in_elts, 1, m_out_elts, 1);
#else
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }
#endif

    return NULL;
}

Error* copy_matblock_cf(Matblock_cf** m_out, const Matblock_cf* m_in)
{
    uint32_t   num_elts;
    Complex_f* m_out_elts;
    Complex_f* m_in_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_cf(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_ccopy(num_elts, m_in_elts, 1, m_out_elts, 1);
#else
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }
#endif

    return NULL;
}

Error* copy_matblock_cd(Matblock_cd** m_out, const Matblock_cd* m_in)
{
    uint32_t   num_elts;
    Complex_d* m_out_elts;
    Complex_d* m_in_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    create_matblock_cd(m_out, m_in->num_mats, m_in->num_rows, m_in->num_cols);
    num_elts = m_in->num_mats*m_in->num_rows*m_in->num_cols;

    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_zcopy(num_elts, m_in_elts, 1, m_out_elts, 1);
#else
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_in_elts[ elt ];
    }
#endif

    return NULL;
}

Error* copy_matblock_block_u32
(
    Matblock_u32**      m_out, 
    const Matblock_u32* m_in,
    uint32_t            mat_offset,
    uint32_t            row_offset,
    uint32_t            col_offset,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t      mat, row, col;
    Matblock_u32* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_u32(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_u32(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_i32
(
    Matblock_i32**      m_out, 
    const Matblock_i32* m_in,
    uint32_t            mat_offset,
    uint32_t            row_offset,
    uint32_t            col_offset,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t      mat, row, col;
    Matblock_i32* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_i32(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_i32(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_i64
(
    Matblock_i64**      m_out, 
    const Matblock_i64* m_in,
    uint32_t            mat_offset,
    uint32_t            row_offset,
    uint32_t            col_offset,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t      mat, row, col;
    Matblock_i64* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_i64(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_i64(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_in,
    uint32_t          mat_offset,
    uint32_t          row_offset,
    uint32_t          col_offset,
    uint32_t          num_mats,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t    mat, row, col;
    Matblock_f* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_f(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_f(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_in,
    uint32_t          mat_offset,
    uint32_t          row_offset,
    uint32_t          col_offset,
    uint32_t          num_mats,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t    mat, row, col;
    Matblock_d* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_d(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_d(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_cf
(
    Matblock_cf**      m_out, 
    const Matblock_cf* m_in,
    uint32_t           mat_offset,
    uint32_t           row_offset,
    uint32_t           col_offset,
    uint32_t           num_mats,
    uint32_t           num_rows,
    uint32_t           num_cols
)
{
    uint32_t     mat, row, col;
    Matblock_cf* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_cf(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_cf(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_cd
(
    Matblock_cd**      m_out, 
    const Matblock_cd* m_in,
    uint32_t           mat_offset,
    uint32_t           row_offset,
    uint32_t           col_offset,
    uint32_t           num_mats,
    uint32_t           num_rows,
    uint32_t           num_cols
)
{
    uint32_t     mat, row, col;
    Matblock_cd* m;

    if (*m_out == m_in)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    if (m_in->num_mats < (mat_offset+num_mats) ||
        m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matblock_cd(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input matblock");
    }

    create_matblock_cd(m_out, num_mats, num_rows, num_cols);
    m = *m_out;

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ mat ][ row ][ col ] = m_in->elts[ mat+mat_offset ]
                                                         [ row+row_offset ]
                                                         [ col+col_offset ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_u32
(
    Matblock_u32*       m_1, 
    uint32_t            mat_offset_1,
    uint32_t            row_offset_1,
    uint32_t            col_offset_1,
    const Matblock_u32* m_2,
    uint32_t            mat_offset_2,
    uint32_t            row_offset_2,
    uint32_t            col_offset_2,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_i32
(
    Matblock_i32*       m_1, 
    uint32_t            mat_offset_1,
    uint32_t            row_offset_1,
    uint32_t            col_offset_1,
    const Matblock_i32* m_2,
    uint32_t            mat_offset_2,
    uint32_t            row_offset_2,
    uint32_t            col_offset_2,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_i64
(
    Matblock_i64*       m_1, 
    uint32_t            mat_offset_1,
    uint32_t            row_offset_1,
    uint32_t            col_offset_1,
    const Matblock_i64* m_2,
    uint32_t            mat_offset_2,
    uint32_t            row_offset_2,
    uint32_t            col_offset_2,
    uint32_t            num_mats,
    uint32_t            num_rows,
    uint32_t            num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_f
(
    Matblock_f*       m_1, 
    uint32_t          mat_offset_1,
    uint32_t          row_offset_1,
    uint32_t          col_offset_1,
    const Matblock_f* m_2,
    uint32_t          mat_offset_2,
    uint32_t          row_offset_2,
    uint32_t          col_offset_2,
    uint32_t          num_mats,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_d
(
    Matblock_d*       m_1, 
    uint32_t          mat_offset_1,
    uint32_t          row_offset_1,
    uint32_t          col_offset_1,
    const Matblock_d* m_2,
    uint32_t          mat_offset_2,
    uint32_t          row_offset_2,
    uint32_t          col_offset_2,
    uint32_t          num_mats,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_cf
(
    Matblock_cf*       m_1, 
    uint32_t           mat_offset_1,
    uint32_t           row_offset_1,
    uint32_t           col_offset_1,
    const Matblock_cf* m_2,
    uint32_t           mat_offset_2,
    uint32_t           row_offset_2,
    uint32_t           col_offset_2,
    uint32_t           num_mats,
    uint32_t           num_rows,
    uint32_t           num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matblock_block_into_matblock_cd
(
    Matblock_cd*       m_1, 
    uint32_t           mat_offset_1,
    uint32_t           row_offset_1,
    uint32_t           col_offset_1,
    const Matblock_cd* m_2,
    uint32_t           mat_offset_2,
    uint32_t           row_offset_2,
    uint32_t           col_offset_2,
    uint32_t           num_mats,
    uint32_t           num_rows,
    uint32_t           num_cols
)
{
    uint32_t mat, row, col;

    if (m_1 == m_2)
    {
        return JWSC_EARG("Copying a matblock into itself is not supported");
    }

    /* Test if the size of the block to copy is in the bounds of the input
     * matblock. */
    if (m_2->num_mats < (mat_offset_2+num_mats) ||
        m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        m_1->num_mats < (mat_offset_1+num_mats) ||
        m_1->num_rows < (row_offset_1+num_rows) ||
        m_1->num_cols < (col_offset_1+num_cols))
    {
        return JWSC_EARG("Block to copy outside bounds of input matblock");
    }

    for (mat = 0; mat < num_mats; mat++)
    {
        for (row = 0; row < num_rows; row++)
        {
            for (col = 0; col < num_cols; col++)
            {
                m_1->elts[ mat+mat_offset_1 ]
                         [ row+row_offset_1 ]
                         [ col+col_offset_1 ] = m_2->elts[ mat+mat_offset_2 ]
                                                         [ row+row_offset_2 ]
                                                         [ col+col_offset_2 ];
            }
        }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_u32
(
    Matblock_u32**      mb_out,
    const Matblock_u32* mb_in,
    const Matrix_u32*   m,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_u32* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_u32(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_i32
(
    Matblock_i32**      mb_out,
    const Matblock_i32* mb_in,
    const Matrix_i32*   m,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_i32* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_i32(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_i64
(
    Matblock_i64**      mb_out,
    const Matblock_i64* mb_in,
    const Matrix_i64*   m,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_i64* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_i64(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_f
(
    Matblock_f**      mb_out,
    const Matblock_f* mb_in,
    const Matrix_f*   m,
    Matblock_matrix   orient,
    uint32_t          index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_f* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_f(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_d
(
    Matblock_d**      mb_out,
    const Matblock_d* mb_in,
    const Matrix_d*   m,
    Matblock_matrix   orient,
    uint32_t          index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_d* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_d(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_cf
(
    Matblock_cf**      mb_out,
    const Matblock_cf* mb_in,
    const Matrix_cf*   m,
    Matblock_matrix    orient,
    uint32_t           index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_cf* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_cf(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_into_matblock_cd
(
    Matblock_cd**      mb_out,
    const Matblock_cd* mb_in,
    const Matrix_cd*   m,
    Matblock_matrix    orient,
    uint32_t           index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matblock_cd* mb;

    num_mats = mb_in->num_mats;
    num_rows = mb_in->num_rows;
    num_cols = mb_in->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb_in->num_mats;
            if (num_rows != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb_in->num_rows;
            if (num_mats != m->num_rows || num_cols != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
            break;
        default:
            num_indices = mb_in->num_cols;
            if (num_mats != m->num_rows || num_rows != m->num_cols)
            {
                return JWSC_EARG("Incompatible matrix and matblock size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matblock");
    }

    if (*mb_out != mb_in)
    {
        copy_matblock_cd(mb_out, mb_in);
    }
    mb = *mb_out;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ index ][ row ][ col ] = m->elts[ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    mb->elts[ mat ][ index ][ col ] = m->elts[ mat ][ col ];
                }
            }
            break;
        default:
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    mb->elts[ mat ][ row ][ index ] = m->elts[ mat ][ row ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_u32
(
    Matrix_u32**        m_out,
    const Matblock_u32* mb,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_u32* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_u32(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_u32(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_u32(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_i32
(
    Matrix_i32**        m_out,
    const Matblock_i32* mb,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_i32* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_i32(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_i32(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_i32(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_i64
(
    Matrix_i64**        m_out,
    const Matblock_i64* mb,
    Matblock_matrix     orient,
    uint32_t            index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_i64* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_i64(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_i64(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_i64(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_f
(
    Matrix_f**        m_out,
    const Matblock_f* mb,
    Matblock_matrix   orient,
    uint32_t          index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_f* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_f(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_f(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_f(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_d
(
    Matrix_d**        m_out,
    const Matblock_d* mb,
    Matblock_matrix   orient,
    uint32_t          index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_d* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_d(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_d(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_d(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_cf
(
    Matrix_cf**        m_out,
    const Matblock_cf* mb,
    Matblock_matrix    orient,
    uint32_t           index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_cf* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_cf(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_cf(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_cf(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

Error* copy_matrix_from_matblock_cd
(
    Matrix_cd**        m_out,
    const Matblock_cd* mb,
    Matblock_matrix    orient,
    uint32_t           index
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t num_indices;

    Matrix_cd* m = NULL;

    num_mats = mb->num_mats;
    num_rows = mb->num_rows;
    num_cols = mb->num_cols;

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            num_indices = mb->num_mats;
            break;
        case MATBLOCK_ROW_MATRIX:
            num_indices = mb->num_rows;
            break;
        default:
            num_indices = mb->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matblock");
    }

    switch (orient)
    {
        case MATBLOCK_MAT_MATRIX:
            create_matrix_cd(m_out, num_rows, num_cols);
            m = *m_out;
            for (row = 0; row < num_rows; row++)
            {
                for (col = 0; col < num_cols; col++)
                {
                    m->elts[ row ][ col ] = mb->elts[ index ][ row ][ col ];
                }
            }
            break;
        case MATBLOCK_ROW_MATRIX:
            create_matrix_cd(m_out, num_mats, num_cols);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (col = 0; col < num_cols; col++)
                {
                     m->elts[ mat ][ col ] = mb->elts[ mat ][ index ][ col ];
                }
            }
            break;
        default:
            create_matrix_cd(m_out, num_mats, num_rows);
            m = *m_out;
            for (mat = 0; mat < num_mats; mat++)
            {
                for (row = 0; row < num_rows; row++)
                {
                    m->elts[ mat ][ row ] = mb->elts[ mat ][ row ][ index ];
                }
            }
    }

    return NULL;
}

void free_matblock_u8(Matblock_u8* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_u32(Matblock_u32* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_i32(Matblock_i32* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_i64(Matblock_i64* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_f(Matblock_f* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_d(Matblock_d* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_cf(Matblock_cf* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}

void free_matblock_cd(Matblock_cd* m)
{
    if (m == NULL) return;

    free(**(m->elts));
    free(*(m->elts));
    free(m->elts);
    free(m);
}