matrix.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/vector/vector.h"
#include "jwsc/matrix/matrix.h"

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


void create_matrix_u32
(
    Matrix_u32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    uint32_t    row;
    uint32_t*   elts;
    Matrix_u32* m;

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

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

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

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

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

void create_matrix_i32
(
    Matrix_i32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    uint32_t    row;
    int32_t*    elts;
    Matrix_i32* m;

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

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

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

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

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

void create_matrix_i64
(
    Matrix_i64** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    uint32_t    row;
    int64_t*    elts;
    Matrix_i64* m;

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

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

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

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

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

void create_matrix_f
(
    Matrix_f** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols
)
{
    uint32_t  row;
    float*    elts;
    Matrix_f* m;

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

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

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

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

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

void create_matrix_d
(
    Matrix_d** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols
)
{
    uint32_t  row;
    double*   elts;
    Matrix_d* m;

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

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

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

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

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

void create_matrix_cf
(
    Matrix_cf** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols
)
{
    uint32_t   row;
    Complex_f* elts;
    Matrix_cf* m;

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

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

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

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

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

void create_matrix_cd
(
    Matrix_cd** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols
)
{
    uint32_t   row;
    Complex_d* elts;
    Matrix_cd* m;

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

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

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

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

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

void create_init_matrix_u32
(
    Matrix_u32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    uint32_t     val
)
{
    uint32_t  num_elts;
    uint32_t  elt;
    uint32_t* m_elts;

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

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

void create_init_matrix_i32
(
    Matrix_i32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    int32_t      val
)
{
    uint32_t num_elts;
    uint32_t elt;
    int32_t* m_elts;

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

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

void create_init_matrix_i64
(
    Matrix_i64** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols,
    int64_t      val
)
{
    uint32_t num_elts;
    uint32_t elt;
    int64_t* m_elts;

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

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

void create_init_matrix_f
(
    Matrix_f** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols,
    float      val
)
{
    uint32_t num_elts;
    uint32_t elt;
    float*   m_elts;

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

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

void create_init_matrix_d
(
    Matrix_d** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols,
    double     val
)
{
    uint32_t num_elts;
    uint32_t elt;
    double*  m_elts;

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

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

void create_init_matrix_cf
(
    Matrix_cf** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols,
    Complex_f   val
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_f* m_elts;

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

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

void create_init_matrix_cd
(
    Matrix_cd** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols,
    Complex_d   val
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_d* m_elts;

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

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

void create_identity_matrix_u32(Matrix_u32** m_out, uint32_t size)
{
    uint32_t elt;

    create_zero_matrix_u32(m_out, size, size);

    for (elt = 0; elt < size; elt++)
    {
        (*m_out)->elts[ elt ][ elt ] = 1;
    }
}

void create_identity_matrix_i32(Matrix_i32** m_out, uint32_t size)
{
    uint32_t elt;

    create_zero_matrix_i32(m_out, size, size);

    for (elt = 0; elt < size; elt++)
    {
        (*m_out)->elts[ elt ][ elt ] = 1;
    }
}

void create_identity_matrix_i64(Matrix_i64** m_out, uint32_t size)
{
    uint32_t elt;

    create_zero_matrix_i64(m_out, size, size);

    for (elt = 0; elt < size; elt++)
    {
        (*m_out)->elts[ elt ][ elt ] = 1;
    }
}

void create_identity_matrix_f(Matrix_f** m_out, uint32_t size)
{
    uint32_t elt;

    create_zero_matrix_f(m_out, size, size);

    for (elt = 0; elt < size; elt++)
    {
        (*m_out)->elts[ elt ][ elt ] = 1;
    }
}

void create_identity_matrix_d(Matrix_d** m_out, uint32_t size)
{
    uint32_t elt;

    create_zero_matrix_d(m_out, size, size);

    for (elt = 0; elt < size; elt++)
    {
        (*m_out)->elts[ elt ][ elt ] = 1;
    }
}

void create_zero_matrix_u32
(
    Matrix_u32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    create_init_matrix_u32(m_out, num_rows, num_cols, 0);
}

void create_zero_matrix_i32
(
    Matrix_i32** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    create_init_matrix_i32(m_out, num_rows, num_cols, 0);
}

void create_zero_matrix_i64
(
    Matrix_i64** m_out, 
    uint32_t     num_rows, 
    uint32_t     num_cols
)
{
    create_init_matrix_i64(m_out, num_rows, num_cols, 0);
}

void create_zero_matrix_f
(
    Matrix_f** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols
)
{
    create_init_matrix_f(m_out, num_rows, num_cols, 0);
}

void create_zero_matrix_d
(
    Matrix_d** m_out, 
    uint32_t   num_rows, 
    uint32_t   num_cols
)
{
    create_init_matrix_d(m_out, num_rows, num_cols, 0);
}

void create_zero_matrix_cf
(
    Matrix_cf** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols
)
{
    Complex_f z = {0};
    create_init_matrix_cf(m_out, num_rows, num_cols, z);
}

void create_zero_matrix_cd
(
    Matrix_cd** m_out, 
    uint32_t    num_rows, 
    uint32_t    num_cols
)
{
    Complex_d z = {0};
    create_init_matrix_cd(m_out, num_rows, num_cols, z);
}

void create_diagnonal_matrix_u32
(
    Matrix_u32**      m_out, 
    const Vector_u32* v
)
{
    uint32_t    n, N;
    Matrix_u32* m;

    N = v->num_elts;

    create_zero_matrix_u32(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_i32
(
    Matrix_i32**      m_out, 
    const Vector_i32* v
)
{
    uint32_t    n, N;
    Matrix_i32* m;

    N = v->num_elts;

    create_zero_matrix_i32(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_i64
(
    Matrix_i64**      m_out, 
    const Vector_i64* v
)
{
    uint32_t    n, N;
    Matrix_i64* m;

    N = v->num_elts;

    create_zero_matrix_i64(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_f
(
    Matrix_f**      m_out, 
    const Vector_f* v
)
{
    uint32_t  n, N;
    Matrix_f* m;

    N = v->num_elts;

    create_zero_matrix_f(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_d
(
    Matrix_d**      m_out, 
    const Vector_d* v
)
{
    uint32_t  n, N;
    Matrix_d* m;

    N = v->num_elts;

    create_zero_matrix_d(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_cf
(
    Matrix_cf**      m_out, 
    const Vector_cf* v
)
{
    uint32_t    n, N;
    Matrix_cf* m;

    N = v->num_elts;

    create_zero_matrix_cf(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_diagonal_matrix_cd
(
    Matrix_cd**      m_out, 
    const Vector_cd* v
)
{
    uint32_t    n, N;
    Matrix_cd* m;

    N = v->num_elts;

    create_zero_matrix_cd(m_out, N, N);
    m = *m_out;

    for (n = 0; n < N; n++)
    {
        m->elts[ n ][ n ] = v->elts[ n ];
    }
}

void create_random_matrix_u32
(
    Matrix_u32** m_out, 
    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_matrix_u32(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_i32
(
    Matrix_i32** m_out, 
    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_matrix_i32(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_i64
(
    Matrix_i64** m_out, 
    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_matrix_i64(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_f
(
    Matrix_f** m_out, 
    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_matrix_f(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_d
(
    Matrix_d** m_out, 
    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_matrix_d(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_cf
(
    Matrix_cf** m_out, 
    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_matrix_cf(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_cd
(
    Matrix_cd** m_out, 
    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_matrix_cd(m_out, num_rows, num_cols);
    m_elts = *((*m_out)->elts);
    num_elts = 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_matrix_u32(Matrix_u32** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_u32(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_i32(Matrix_i32** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_i32(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_i64(Matrix_i64** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_i64(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_f(Matrix_f** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_f(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_d(Matrix_d** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_d(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_cf(Matrix_cf** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_cf(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_cd(Matrix_cd** m_out, const Matrix_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 matrix into itself is not supported");
    }

    create_matrix_cd(m_out, m_in->num_rows, m_in->num_cols);
    num_elts = 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_matrix_block_u32
(
    Matrix_u32**      m_out, 
    const Matrix_u32* m_in,
    uint32_t          row_offset,
    uint32_t          col_offset,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t    row, col;
    Matrix_u32* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_u32(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_i32
(
    Matrix_i32**      m_out, 
    const Matrix_i32* m_in,
    uint32_t          row_offset,
    uint32_t          col_offset,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t    row, col;
    Matrix_i32* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_i32(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_i64
(
    Matrix_i64**      m_out, 
    const Matrix_i64* m_in,
    uint32_t          row_offset,
    uint32_t          col_offset,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t    row, col;
    Matrix_i64* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_i64(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_in,
    uint32_t        row_offset,
    uint32_t        col_offset,
    uint32_t        num_rows,
    uint32_t        num_cols
)
{
    uint32_t  row, col;
    Matrix_f* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_f(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_in,
    uint32_t        row_offset,
    uint32_t        col_offset,
    uint32_t        num_rows,
    uint32_t        num_cols
)
{
    uint32_t  row, col;
    Matrix_d* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_d(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_in,
    uint32_t         row_offset,
    uint32_t         col_offset,
    uint32_t         num_rows,
    uint32_t         num_cols
)
{
    uint32_t   row, col;
    Matrix_cf* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_cf(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_in,
    uint32_t         row_offset,
    uint32_t         col_offset,
    uint32_t         num_rows,
    uint32_t         num_cols
)
{
    uint32_t   row, col;
    Matrix_cd* m;

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

    if (m_in->num_rows < (row_offset+num_rows) || 
        m_in->num_cols < (col_offset+num_cols))
    {
        free_matrix_cd(*m_out); m_out = NULL;
        return JWSC_EARG("Section to copy outside bounds of input 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 ] 
                = m_in->elts[ row+row_offset ][ col+col_offset ];
        }
    }

    return NULL;
}

Error* copy_matrix_block_into_matrix_u32
(
    Matrix_u32*       m_1, 
    uint32_t          row_offset_1,
    uint32_t          col_offset_1,
    const Matrix_u32* m_2,
    uint32_t          row_offset_2,
    uint32_t          col_offset_2,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_i32
(
    Matrix_i32*       m_1, 
    uint32_t          row_offset_1,
    uint32_t          col_offset_1,
    const Matrix_i32* m_2,
    uint32_t          row_offset_2,
    uint32_t          col_offset_2,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_i64
(
    Matrix_i64*       m_1, 
    uint32_t          row_offset_1,
    uint32_t          col_offset_1,
    const Matrix_i64* m_2,
    uint32_t          row_offset_2,
    uint32_t          col_offset_2,
    uint32_t          num_rows,
    uint32_t          num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_f
(
    Matrix_f*       m_1, 
    uint32_t        row_offset_1,
    uint32_t        col_offset_1,
    const Matrix_f* m_2,
    uint32_t        row_offset_2,
    uint32_t        col_offset_2,
    uint32_t        num_rows,
    uint32_t        num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_d
(
    Matrix_d*       m_1, 
    uint32_t        row_offset_1,
    uint32_t        col_offset_1,
    const Matrix_d* m_2,
    uint32_t        row_offset_2,
    uint32_t        col_offset_2,
    uint32_t        num_rows,
    uint32_t        num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_cf
(
    Matrix_cf*       m_1, 
    uint32_t         row_offset_1,
    uint32_t         col_offset_1,
    const Matrix_cf* m_2,
    uint32_t         row_offset_2,
    uint32_t         col_offset_2,
    uint32_t         num_rows,
    uint32_t         num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_matrix_block_into_matrix_cd
(
    Matrix_cd*       m_1, 
    uint32_t         row_offset_1,
    uint32_t         col_offset_1,
    const Matrix_cd* m_2,
    uint32_t         row_offset_2,
    uint32_t         col_offset_2,
    uint32_t         num_rows,
    uint32_t         num_cols
)
{
    uint32_t row, col;

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

    /* Test if the size of the block to copy is in the bounds of the input
     * matrix. */
    if (m_2->num_rows < (row_offset_2+num_rows) ||
        m_2->num_cols < (col_offset_2+num_cols) ||
        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 matrix");
    }

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

    return NULL;
}

Error* copy_vector_into_matrix_u32
(
    Matrix_u32**      m_out,
    const Matrix_u32* m_in,
    const Vector_u32* v,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_u32* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_u32(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_i32
(
    Matrix_i32**      m_out,
    const Matrix_i32* m_in,
    const Vector_i32* v,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_i32* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_i32(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_i64
(
    Matrix_i64**      m_out,
    const Matrix_i64* m_in,
    const Vector_i64* v,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_i64* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_i64(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_f
(
    Matrix_f**      m_out,
    const Matrix_f* m_in,
    const Vector_f* v,
    Matrix_vector   orient,
    uint32_t        index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_f* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_f(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_d
(
    Matrix_d**      m_out,
    const Matrix_d* m_in,
    const Vector_d* v,
    Matrix_vector   orient,
    uint32_t        index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_d* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_d(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_cf
(
    Matrix_cf**      m_out,
    const Matrix_cf* m_in,
    const Vector_cf* v,
    Matrix_vector    orient,
    uint32_t         index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_cf* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_cf(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_into_matrix_cd
(
    Matrix_cd**      m_out,
    const Matrix_cd* m_in,
    const Vector_cd* v,
    Matrix_vector    orient,
    uint32_t         index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Matrix_cd* m;

    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m_in->num_rows;
            if (num_cols != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
            break;
        default:
            num_indices = m_in->num_cols;
            if (num_rows != v->num_elts)
            {
                return JWSC_EARG("Incompatible vector and matrix size");
            }
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy into matrix");
    }

    if (*m_out != m_in)
    {
        copy_matrix_cd(m_out, m_in);
    }
    m = *m_out;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            for (col = 0; col < num_cols; col++)
            {
                m->elts[ index ][ col ] = v->elts[ col ];
            }
            break;
        default:
            for (row = 0; row < num_rows; row++)
            {
                m->elts[ row ][ index ] = v->elts[ row ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_u32
(
    Vector_u32**      v_out,
    const Matrix_u32* m,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_u32* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_u32(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_u32(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_i32
(
    Vector_i32**      v_out,
    const Matrix_i32* m,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_i32* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_i32(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_i32(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_i64
(
    Vector_i64**      v_out,
    const Matrix_i64* m,
    Matrix_vector     orient,
    uint32_t          index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_i64* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_i64(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_i64(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_f
(
    Vector_f**      v_out,
    const Matrix_f* m,
    Matrix_vector   orient,
    uint32_t        index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_f* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_f(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_f(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_d
(
    Vector_d**      v_out,
    const Matrix_d* m,
    Matrix_vector   orient,
    uint32_t        index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_d* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_d(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_d(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_cf
(
    Vector_cf**      v_out,
    const Matrix_cf* m,
    Matrix_vector    orient,
    uint32_t         index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_cf* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_cf(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_cf(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

Error* copy_vector_from_matrix_cd
(
    Vector_cd**      v_out,
    const Matrix_cd* m,
    Matrix_vector    orient,
    uint32_t         index
)
{
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    uint32_t num_indices;

    Vector_cd* v;

    num_rows = m->num_rows;
    num_cols = m->num_cols;

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            num_indices = m->num_rows;
            break;
        default:
            num_indices = m->num_cols;
    }
    if (index >= num_indices)
    {
        return JWSC_EARG("Invalid index to copy from matrix");
    }

    switch (orient)
    {
        case MATRIX_ROW_VECTOR:
            create_vector_cd(v_out, num_cols);
            v = *v_out;
            for (col = 0; col < num_cols; col++)
            {
                v->elts[ col ] = m->elts[ index ][ col ];
            }
            break;
        default:
            create_vector_cd(v_out, num_rows);
            v = *v_out;
            for (row = 0; row < num_rows; row++)
            {
                v->elts[ row ] = m->elts[ row ][ index ];
            }
    }

    return NULL;
}

void free_matrix_u32(Matrix_u32* m)
{
    if (m == NULL) return;

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

void free_matrix_i32(Matrix_i32* m)
{
    if (m == NULL) return;

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

void free_matrix_i64(Matrix_i64* m)
{
    if (m == NULL) return;

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

void free_matrix_f(Matrix_f* m)
{
    if (m == NULL) return;

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

void free_matrix_d(Matrix_d* m)
{
    if (m == NULL) return;

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

void free_matrix_cf(Matrix_cf* m)
{
    if (m == NULL) return;

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

void free_matrix_cd(Matrix_cd* m)
{
    if (m == NULL) return;

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