matrix_math.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/constants.h"
#include "jwsc/math/blas.h"
#include "jwsc/math/lapack.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector_math.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"

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


void add_scalar_to_matrix_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_in, 
    float           a
)
{
    uint32_t  num_elts;
    uint32_t  elt;
    float*    m_in_elts;
    float*    m_out_elts;

    num_elts = m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matrix_f(m_out, 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 ] = a + m_in_elts[ elt ];
    }
}

void add_scalar_to_matrix_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_in, 
    double          a
)
{
    uint32_t  num_elts;
    uint32_t  elt;
    double*   m_in_elts;
    double*   m_out_elts;

    num_elts = m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matrix_d(m_out, 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 ] = a + m_in_elts[ elt ];
    }
}

void add_scalar_to_matrix_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_in, 
    Complex_f        z
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_f* m_in_elts;
    Complex_f* m_out_elts;

    num_elts = m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matrix_cf(m_out, 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 ].r = z.r + m_in_elts[ elt ].r;
        m_out_elts[ elt ].i = z.i + m_in_elts[ elt ].i;
    }
}

void add_scalar_to_matrix_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_in, 
    Complex_d        z
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Complex_d* m_in_elts;
    Complex_d* m_out_elts;

    num_elts = m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matrix_cd(m_out, 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 ].r = z.r + m_in_elts[ elt ].r;
        m_out_elts[ elt ].i = z.i + m_in_elts[ elt ].i;
    }
}

void multiply_matrix_by_scalar_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_in, 
    float           a
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
    float*   m_out_elts;
    float*   m_in_elts;
#endif

    num_elts = m_in->num_rows*m_in->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out != m_in)
    {
        copy_matrix_f(m_out, m_in);
    }

    blas_sscal(num_elts, a, *((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matrix_f(m_out, 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 ] = a * m_in_elts[ elt ];
    }
#endif
}

void multiply_matrix_by_scalar_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_in, 
    double          a
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
    double*  m_out_elts;
    double*  m_in_elts;
#endif

    num_elts = m_in->num_rows*m_in->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out != m_in)
    {
        copy_matrix_d(m_out, m_in);
    }

    blas_dscal(num_elts, a, *((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matrix_d(m_out, 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 ] = a * m_in_elts[ elt ];
    }
#endif
}

void multiply_matrix_by_scalar_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_in, 
    Complex_f        z
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t   elt;
    float      r, i;
    Complex_f* m_out_elts;
    Complex_f* m_in_elts;
#endif

    num_elts = m_in->num_rows*m_in->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out != m_in)
    {
        copy_matrix_cf(m_out, m_in);
    }

    blas_cscal(num_elts, z, *((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matrix_cf(m_out, 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++)
    {
        r = z.r*m_in_elts[ elt ].r - z.i*m_in_elts[ elt ].i;
        i = z.i*m_in_elts[ elt ].r + z.r*m_in_elts[ elt ].i;

        m_out_elts[ elt ].r = r;
        m_out_elts[ elt ].i = i;
    }
#endif
}

void multiply_matrix_by_scalar_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_in, 
    Complex_d        z
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t   elt;
    float      r, i;
    Complex_d* m_out_elts;
    Complex_d* m_in_elts;
#endif

    num_elts = m_in->num_rows*m_in->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out != m_in)
    {
        copy_matrix_cd(m_out, m_in);
    }

    blas_zscal(num_elts, z, *((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matrix_cd(m_out, 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++)
    {
        r = z.r*m_in_elts[ elt ].r - z.i*m_in_elts[ elt ].i;
        i = z.i*m_in_elts[ elt ].r + z.r*m_in_elts[ elt ].i;

        m_out_elts[ elt ].r = r;
        m_out_elts[ elt ].i = i;
    }
#endif
}

Error* multiply_matrices_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_1, 
    const Matrix_f* m_2
)
{
    uint32_t  num_rows, num_cols;
    Matrix_f* m;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t row, col;
    uint32_t elt;
    uint32_t length;
    float    sum;
#endif

    if (m_1->num_cols != m_2->num_rows)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have compatible dimensions");
    }

    num_rows = m_1->num_rows;
    num_cols = m_2->num_cols;

    m = (*m_out == m_1 || *m_out == m_2) ? NULL : *m_out;
    create_matrix_f(&m, num_rows, num_cols);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_sgemm('N', 'N', m_1->num_rows, m_2->num_cols, m_1->num_cols, 1, 
            *(m_1->elts), m_1->num_cols, *(m_2->elts), m_2->num_cols, 0, 
            *(m->elts), m->num_cols);
#else
    length = m_1->num_cols;
    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            sum = 0;
            for (elt = 0; elt < length; elt++)
            {
                sum += m_1->elts[ row ][ elt ] * m_2->elts[ elt ][ col ];
            }
            m->elts[ row ][ col ] = sum;
        }
    }
#endif

    if (*m_out == m_1 || *m_out == m_2)
    {
        copy_matrix_f(m_out, m);
        free_matrix_f(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

Error* multiply_matrices_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_1, 
    const Matrix_d* m_2
)
{
    uint32_t  num_rows, num_cols;
    Matrix_d* m;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t row, col;
    uint32_t elt;
    uint32_t length;
    float    sum;
#endif

    if (m_1->num_cols != m_2->num_rows)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have compatible dimensions");
    }

    num_rows = m_1->num_rows;
    num_cols = m_2->num_cols;

    m = (*m_out == m_1 || *m_out == m_2) ? NULL : *m_out;
    create_matrix_d(&m, num_rows, num_cols);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_dgemm('N', 'N', m_1->num_rows, m_2->num_cols, m_1->num_cols, 1, 
            *(m_1->elts), m_1->num_cols, *(m_2->elts), m_2->num_cols, 0, 
            *(m->elts), m->num_cols);
#else
    length = m_1->num_cols;
    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            sum = 0;
            for (elt = 0; elt < length; elt++)
            {
                sum += m_1->elts[ row ][ elt ] * m_2->elts[ elt ][ col ];
            }
            m->elts[ row ][ col ] = sum;
        }
    }
#endif

    if (*m_out == m_1 || *m_out == m_2)
    {
        copy_matrix_d(m_out, m);
        free_matrix_d(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

Error* multiply_matrices_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_1, 
    const Matrix_cf* m_2
)
{
    uint32_t   num_rows, num_cols;
    Matrix_cf* m;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_f alpha;
    Complex_f beta;
#else
    uint32_t  row, col;
    uint32_t  elt;
    uint32_t  length;
    Complex_f sum;
#endif

    if (m_1->num_cols != m_2->num_rows)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cf(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have compatible dimensions");
    }

    num_rows = m_1->num_rows;
    num_cols = m_2->num_cols;

    m = (*m_out == m_1 || *m_out == m_2) ? NULL : *m_out;
    create_matrix_cf(&m, num_rows, num_cols);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    alpha.r = 1; alpha.i = 1;
    beta.r  = 0; beta.i  = 0;
    blas_cgemm('N', 'N', m_1->num_rows, m_2->num_cols, m_1->num_cols, alpha, 
            *(m_1->elts), m_1->num_cols, *(m_2->elts), m_2->num_cols, beta, 
            *(m->elts), m->num_cols);
#else
    length = m_1->num_cols;
    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            sum.r = 0; 
            sum.i = 0;
            for (elt = 0; elt < length; elt++)
            {
                sum.r += m_1->elts[ row ][ elt ].r * m_2->elts[ elt ][ col ].r -
                         m_1->elts[ row ][ elt ].i * m_2->elts[ elt ][ col ].i;
                sum.i += m_1->elts[ row ][ elt ].i * m_2->elts[ elt ][ col ].r +
                         m_1->elts[ row ][ elt ].r * m_2->elts[ elt ][ col ].i;
            }
            m->elts[ row ][ col ].r = sum.r;
            m->elts[ row ][ col ].i = sum.i;
        }
    }
#endif

    if (*m_out == m_1 || *m_out == m_2)
    {
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

Error* multiply_matrices_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_1, 
    const Matrix_cd* m_2
)
{
    uint32_t   num_rows, num_cols;
    Matrix_cd* m;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_d alpha;
    Complex_d beta;
#else
    uint32_t  row, col;
    uint32_t  elt;
    uint32_t  length;
    Complex_d sum;
#endif

    if (m_1->num_cols != m_2->num_rows)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cd(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have compatible dimensions");
    }

    num_rows = m_1->num_rows;
    num_cols = m_2->num_cols;

    m = (*m_out == m_1 || *m_out == m_2) ? NULL : *m_out;
    create_matrix_cd(&m, num_rows, num_cols);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    alpha.r = 1; alpha.i = 1;
    beta.r  = 0; beta.i  = 0;
    blas_zgemm('N', 'N', m_1->num_rows, m_2->num_cols, m_1->num_cols, alpha, 
            *(m_1->elts), m_1->num_cols, *(m_2->elts), m_2->num_cols, beta, 
            *(m->elts), m->num_cols);
#else
    length = m_1->num_cols;
    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            sum.r = 0; 
            sum.i = 0;
            for (elt = 0; elt < length; elt++)
            {
                sum.r += m_1->elts[ row ][ elt ].r * m_2->elts[ elt ][ col ].r -
                         m_1->elts[ row ][ elt ].i * m_2->elts[ elt ][ col ].i;
                sum.i += m_1->elts[ row ][ elt ].i * m_2->elts[ elt ][ col ].r +
                         m_1->elts[ row ][ elt ].r * m_2->elts[ elt ][ col ].i;
            }
            m->elts[ row ][ col ].r = sum.r;
            m->elts[ row ][ col ].i = sum.i;
        }
    }
#endif

    if (*m_out == m_1 || *m_out == m_2)
    {
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        *m_out = m;
    }

    return NULL;
}

Error* multiply_matrices_elt_wise_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_1, 
    const Matrix_f* m_2
)
{
    uint32_t  num_rows, num_cols;
    uint32_t  num_elts;
    uint32_t  elt;
    float*    m_out_elts;
    float*    m_1_elts;
    float*    m_2_elts;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows ||  num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have the same dimensions");
    }

    create_matrix_f(m_out, num_rows, num_cols);

    m_out_elts = *((*m_out)->elts);
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);

    num_elts = num_rows*num_cols;
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_1_elts[ elt ]*m_2_elts[ elt ];
    }

    return NULL;
}

Error* multiply_matrices_elt_wise_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_1, 
    const Matrix_d* m_2
)
{
    uint32_t  num_rows, num_cols;
    uint32_t  num_elts;
    uint32_t  elt;
    double*   m_out_elts;
    double*   m_1_elts;
    double*   m_2_elts;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows ||  num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have the same dimensions");
    }

    create_matrix_d(m_out, num_rows, num_cols);

    m_out_elts = *((*m_out)->elts);
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);

    num_elts = num_rows*num_cols;
    for (elt = 0; elt < num_elts; elt++)
    {
        m_out_elts[ elt ] = m_1_elts[ elt ]*m_2_elts[ elt ];
    }

    return NULL;
}

Error* multiply_matrices_elt_wise_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_1, 
    const Matrix_cf* m_2
)
{
    uint32_t   num_rows, num_cols;
    uint32_t   num_elts;
    uint32_t   elt;
    float      r, i;
    Complex_f* m_out_elts;
    Complex_f* m_1_elts;
    Complex_f* m_2_elts;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows ||  num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cf(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have the same dimensions");
    }

    create_matrix_cf(m_out, num_rows, num_cols);

    m_out_elts = *((*m_out)->elts);
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);

    num_elts = num_rows*num_cols;
    for (elt = 0; elt < num_elts; elt++)
    {
        r = m_1_elts[ elt ].r * m_2_elts[ elt ].r -
            m_1_elts[ elt ].i * m_2_elts[ elt ].i;
        i = m_1_elts[ elt ].i * m_2_elts[ elt ].r +
            m_1_elts[ elt ].r * m_2_elts[ elt ].i;

        m_out_elts[ elt ].r = r;
        m_out_elts[ elt ].i = i;
    }

    return NULL;
}

Error* multiply_matrices_elt_wise_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_1, 
    const Matrix_cd* m_2
)
{
    uint32_t   num_rows, num_cols;
    uint32_t   num_elts;
    uint32_t   elt;
    float      r, i;
    Complex_d* m_out_elts;
    Complex_d* m_1_elts;
    Complex_d* m_2_elts;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows ||  num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cd(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrices do not have the same dimensions");
    }

    create_matrix_cd(m_out, num_rows, num_cols);

    m_out_elts = *((*m_out)->elts);
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);

    num_elts = num_rows*num_cols;
    for (elt = 0; elt < num_elts; elt++)
    {
        r = m_1_elts[ elt ].r * m_2_elts[ elt ].r -
            m_1_elts[ elt ].i * m_2_elts[ elt ].i;
        i = m_1_elts[ elt ].i * m_2_elts[ elt ].r +
            m_1_elts[ elt ].r * m_2_elts[ elt ].i;

        m_out_elts[ elt ].r = r;
        m_out_elts[ elt ].i = i;
    }

    return NULL;
}

void transpose_matrix_f(Matrix_f** m_out, const Matrix_f* m_in)
{
    uint32_t  num_rows, num_cols;
    uint32_t  row, col;
    Matrix_f* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_f(&m, num_rows, num_cols);

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

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

void transpose_matrix_d(Matrix_d** m_out, const Matrix_d* m_in)
{
    uint32_t  num_rows, num_cols;
    uint32_t  row, col;
    Matrix_d* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_d(&m, num_rows, num_cols);

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

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

void transpose_matrix_cf(Matrix_cf** m_out, const Matrix_cf* m_in)
{
    uint32_t   num_rows, num_cols;
    uint32_t   row, col;
    Matrix_cf* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cf(&m, num_rows, num_cols);

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

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

void transpose_matrix_cd(Matrix_cd** m_out, const Matrix_cd* m_in)
{
    uint32_t   num_rows, num_cols;
    uint32_t   row, col;
    Matrix_cd* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matrix_cd(&m, num_rows, num_cols);

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

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

Error* multiply_matrix_and_vector_f
(
    Vector_f**      v_out, 
    const Matrix_f* m_in, 
    const Vector_f* v_in
)
{
    Vector_f* v;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    float    sum;
#endif

    if (m_in->num_cols != v_in->num_elts)
    {
        if (*v_out != v_in)
        {
            free_vector_f(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Matrix and vector do not have compatible dimensions");
    }

    v = (*v_out == v_in) ? NULL : *v_out;
    create_vector_f(&v, m_in->num_rows);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_sgemv('N', m_in->num_rows, m_in->num_cols, 1, *(m_in->elts), 
            m_in->num_cols, v_in->elts, 1, 0, v->elts, 1);
#else
    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    for (row = 0; row < num_rows; row++)
    {
        sum = 0;
        for (col = 0; col < num_cols; col++)
        {
            sum += m_in->elts[ row ][ col ] * v_in->elts[ col ];
        }
        v->elts[ row ] = sum;
    }
#endif

    if (*v_out == v_in)
    {
        copy_vector_f(v_out, v);
        free_vector_f(v);
    }
    else
    {
        *v_out = v;
    }

    return NULL;
}

Error* multiply_matrix_and_vector_d
(
    Vector_d**      v_out, 
    const Matrix_d* m_in, 
    const Vector_d* v_in
)
{
    Vector_d* v;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t num_rows, num_cols;
    uint32_t row, col;
    double   sum;
#endif

    if (m_in->num_cols != v_in->num_elts)
    {
        if (*v_out != v_in)
        {
            free_vector_d(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Matrix and vector do not have compatible dimensions");
    }

    v = (*v_out == v_in) ? NULL : *v_out;
    create_vector_d(&v, m_in->num_rows);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    blas_dgemv('N', m_in->num_rows, m_in->num_cols, 1, *(m_in->elts), 
            m_in->num_cols, v_in->elts, 1, 0, v->elts, 1);
#else
    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    for (row = 0; row < num_rows; row++)
    {
        sum = 0;
        for (col = 0; col < num_cols; col++)
        {
            sum += m_in->elts[ row ][ col ] * v_in->elts[ col ];
        }
        v->elts[ row ] = sum;
    }
#endif

    if (*v_out == v_in)
    {
        copy_vector_d(v_out, v);
        free_vector_d(v);
    }
    else
    {
        *v_out = v;
    }

    return NULL;
}

Error* multiply_matrix_and_vector_cf
(
    Vector_cf**      v_out, 
    const Matrix_cf* m_in, 
    const Vector_cf* v_in
)
{
    Vector_cf* v;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_f alpha;
    Complex_f beta;
#else
    uint32_t  num_rows, num_cols;
    uint32_t  row, col;
    Complex_f sum;
#endif

    if (m_in->num_cols != v_in->num_elts)
    {
        if (*v_out != v_in)
        {
            free_vector_cf(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Matrix and vector do not have compatible dimensions");
    }

    v = (*v_out == v_in) ? NULL : *v_out;
    create_vector_cf(&v, m_in->num_rows);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    alpha.r = 1; alpha.i = 1;
    beta.r  = 0; beta.i  = 0;
    blas_cgemv('N', m_in->num_rows, m_in->num_cols, alpha, *(m_in->elts), 
            m_in->num_cols, v_in->elts, 1, beta, v->elts, 1);
#else
    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    for (row = 0; row < num_rows; row++)
    {
        sum.r = 0;
        sum.i = 0;
        for (col = 0; col < num_cols; col++)
        {
            sum.r += m_in->elts[ row ][ col ].r * v_in->elts[ col ].r -
                     m_in->elts[ row ][ col ].i * v_in->elts[ col ].i;
            sum.i += m_in->elts[ row ][ col ].i * v_in->elts[ col ].r +
                     m_in->elts[ row ][ col ].r * v_in->elts[ col ].i;
        }
        v->elts[ row ].r = sum.r;
        v->elts[ row ].i = sum.i;
    }
#endif

    if (*v_out == v_in)
    {
        copy_vector_cf(v_out, v);
        free_vector_cf(v);
    }
    else
    {
        *v_out = v;
    }

    return NULL;
}

Error* multiply_matrix_and_vector_cd
(
    Vector_cd**      v_out, 
    const Matrix_cd* m_in, 
    const Vector_cd* v_in
)
{
    Vector_cd* v;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_d alpha;
    Complex_d beta;
#else
    uint32_t  num_rows, num_cols;
    uint32_t  row, col;
    Complex_d sum;
#endif

    if (m_in->num_cols != v_in->num_elts)
    {
        if (*v_out != v_in)
        {
            free_vector_cd(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Matrix and vector do not have compatible dimensions");
    }

    v = (*v_out == v_in) ? NULL : *v_out;
    create_vector_cd(&v, m_in->num_rows);

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    alpha.r = 1; alpha.i = 1;
    beta.r  = 0; beta.i  = 0;
    blas_zgemv('N', m_in->num_rows, m_in->num_cols, alpha, *(m_in->elts), 
            m_in->num_cols, v_in->elts, 1, beta, v->elts, 1);
#else
    num_rows = m_in->num_rows;
    num_cols = m_in->num_cols;

    for (row = 0; row < num_rows; row++)
    {
        sum.r = 0;
        sum.i = 0;
        for (col = 0; col < num_cols; col++)
        {
            sum.r += m_in->elts[ row ][ col ].r * v_in->elts[ col ].r -
                     m_in->elts[ row ][ col ].i * v_in->elts[ col ].i;
            sum.i += m_in->elts[ row ][ col ].i * v_in->elts[ col ].r +
                     m_in->elts[ row ][ col ].r * v_in->elts[ col ].i;
        }
        v->elts[ row ].r = sum.r;
        v->elts[ row ].i = sum.i;
    }
#endif

    if (*v_out == v_in)
    {
        copy_vector_cd(v_out, v);
        free_vector_cd(v);
    }
    else
    {
        *v_out = v;
    }

    return NULL;
}




Error* add_matrices_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_1, 
    const Matrix_f* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
    float*   m_1_elts;
    float*   m_2_elts;
    float*   m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_saxpy(num_elts, 1, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        blas_saxpy(num_elts, 1, *(m_1->elts), 1, *(m_2->elts), 1);
    }
    else
    {
        copy_matrix_f(m_out, m_2);
        blas_saxpy(num_elts, 1, *(m_1->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ] += m_2_elts[ elt ];
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ] += m_1_elts[ elt ];
        }
    }
    else
    {
        create_matrix_f(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ] = m_1_elts[ elt ] + m_2_elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* add_matrices_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_1, 
    const Matrix_d* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
    double*  m_1_elts;
    double*  m_2_elts;
    double*  m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_daxpy(num_elts, 1, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        blas_daxpy(num_elts, 1, *(m_1->elts), 1, *(m_2->elts), 1);
    }
    else
    {
        copy_matrix_d(m_out, m_2);
        blas_daxpy(num_elts, 1, *(m_1->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ] += m_2_elts[ elt ];
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ] += m_1_elts[ elt ];
        }
    }
    else
    {
        create_matrix_d(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ] = m_1_elts[ elt ] + m_2_elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* add_matrices_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_1, 
    const Matrix_cf* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_f z = {1, 0};
#else
    uint32_t elt;
    Complex_f* m_1_elts;
    Complex_f* m_2_elts;
    Complex_f* m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cf(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_caxpy(num_elts, z, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        blas_caxpy(num_elts, z, *(m_1->elts), 1, *(m_2->elts), 1);
    }
    else
    {
        copy_matrix_cf(m_out, m_2);
        blas_caxpy(num_elts, z, *(m_1->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_1_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_2_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
    else
    {
        create_matrix_cf(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_out_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* add_matrices_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_1, 
    const Matrix_cd* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_d z = {1, 0};
#else
    uint32_t elt;
    Complex_d* m_1_elts;
    Complex_d* m_2_elts;
    Complex_d* m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cd(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_zaxpy(num_elts, z, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        blas_zaxpy(num_elts, z, *(m_1->elts), 1, *(m_2->elts), 1);
    }
    else
    {
        copy_matrix_cd(m_out, m_2);
        blas_zaxpy(num_elts, z, *(m_1->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_1_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_2_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
    else
    {
        create_matrix_cd(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ].r = m_1_elts[ elt ].r + m_2_elts[ elt ].r;
            m_out_elts[ elt ].i = m_1_elts[ elt ].i + m_2_elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* subtract_matrices_f
(
    Matrix_f**      m_out, 
    const Matrix_f* m_1, 
    const Matrix_f* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matrix_f* m = NULL;
#else
    uint32_t elt;
    float*   m_1_elts;
    float*   m_2_elts;
    float*   m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_saxpy(num_elts, -1, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        copy_matrix_f(&m, m_1);
        blas_saxpy(num_elts, -1, *(m_2->elts), 1, *(m->elts), 1);
        copy_matrix_f(m_out, m);
        free_matrix_f(m);
    }
    else
    {
        copy_matrix_f(m_out, m_1);
        blas_saxpy(num_elts, -1, *(m_2->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
    else
    {
        create_matrix_f(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* subtract_matrices_d
(
    Matrix_d**      m_out, 
    const Matrix_d* m_1, 
    const Matrix_d* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matrix_d* m = NULL;
#else
    uint32_t elt;
    double*  m_1_elts;
    double*  m_2_elts;
    double*  m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_daxpy(num_elts, -1, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        copy_matrix_d(&m, m_1);
        blas_daxpy(num_elts, -1, *(m_2->elts), 1, *(m->elts), 1);
        copy_matrix_d(m_out, m);
        free_matrix_d(m);
    }
    else
    {
        copy_matrix_d(m_out, m_1);
        blas_daxpy(num_elts, -1, *(m_2->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
    else
    {
        create_matrix_d(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ] = m_1_elts[ elt ] - m_2_elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* subtract_matrices_cf
(
    Matrix_cf**      m_out, 
    const Matrix_cf* m_1, 
    const Matrix_cf* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matrix_cf* m = NULL;
    Complex_f  z = {-1, 0};
#else
    uint32_t   elt;
    Complex_f* m_1_elts;
    Complex_f* m_2_elts;
    Complex_f* m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cf(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_caxpy(num_elts, z, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        copy_matrix_cf(&m, m_1);
        blas_caxpy(num_elts, z, *(m_2->elts), 1, *(m->elts), 1);
        copy_matrix_cf(m_out, m);
        free_matrix_cf(m);
    }
    else
    {
        copy_matrix_cf(m_out, m_1);
        blas_caxpy(num_elts, z, *(m_2->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_1_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_2_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
    else
    {
        create_matrix_cf(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_out_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* subtract_matrices_cd
(
    Matrix_cd**      m_out, 
    const Matrix_cd* m_1, 
    const Matrix_cd* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matrix_cd* m = NULL;
    Complex_d  z = {-1, 0};
#else
    uint32_t   elt;
    Complex_d* m_1_elts;
    Complex_d* m_2_elts;
    Complex_d* m_out_elts;
#endif

    if (m_1->num_rows != m_2->num_rows || m_1->num_cols != m_2->num_cols)
    {
        if (*m_out != m_1 && *m_out != m_2)
        {
            free_matrix_cd(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matrices requires equal dimensions");
    }

    num_elts = m_1->num_rows*m_1->num_cols;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*m_out == m_1)
    {
        blas_zaxpy(num_elts, z, *(m_2->elts), 1, *(m_1->elts), 1);
    }
    else if (*m_out == m_2)
    {
        copy_matrix_cd(&m, m_1);
        blas_zaxpy(num_elts, z, *(m_2->elts), 1, *(m->elts), 1);
        copy_matrix_cd(m_out, m);
        free_matrix_cd(m);
    }
    else
    {
        copy_matrix_cd(m_out, m_1);
        blas_zaxpy(num_elts, z, *(m_2->elts), 1, *((*m_out)->elts), 1);
    }
#else
    m_1_elts   = *(m_1->elts);
    m_2_elts   = *(m_2->elts);
    m_out_elts = *((*m_out)->elts);

    if (*m_out == m_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_1_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_1_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
    else if (*m_out == m_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            m_2_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_2_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
    else
    {
        create_matrix_cd(m_out, m_1->num_rows, m_1->num_cols);
        for (elt = 0; elt < num_elts; elt++)
        {
            m_out_elts[ elt ].r = m_1_elts[ elt ].r - m_2_elts[ elt ].r;
            m_out_elts[ elt ].i = m_1_elts[ elt ].i - m_2_elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

void calc_matrix_asum_f(float* sum_out, const Matrix_f* m)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *sum_out = blas_sasum(m->num_rows*m->num_cols, *(m->elts), 1);
#else
    uint32_t num_elts;
    uint32_t elt;
    float*   m_elts;

    num_elts = m->num_rows*m->num_cols;
    m_elts   = *(m->elts);
    *sum_out = 0;

    for (elt = 0; elt < num_elts; elt++)
    {
        (*sum_out) += (float)fabs(m_elts[ elt ]);
    }
#endif
}

void calc_matrix_asum_d(double* sum_out, const Matrix_d* m)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *sum_out = blas_dasum(m->num_rows*m->num_cols, *(m->elts), 1);
#else
    uint32_t num_elts;
    uint32_t elt;
    double*  m_elts;

    num_elts = m->num_rows*m->num_cols;
    m_elts   = *(m->elts);
    *sum_out = 0;

    for (elt = 0; elt < num_elts; elt++)
    {
        (*sum_out) += fabs(m_elts[ elt ]);
    }
#endif
}

void normalize_matrix_sum_f(Matrix_f** m_out, const Matrix_f* m_in)
{
    float sum;

    calc_matrix_asum_f(&sum, m_in);

    if (sum > 0)
    {
        multiply_matrix_by_scalar_f(m_out, m_in, 1.0f / sum);
    } 
    else 
    {
        copy_matrix_f(m_out, m_in);
    }
}

void normalize_matrix_sum_d(Matrix_d** m_out, const Matrix_d* m_in)
{
    double sum;

    calc_matrix_asum_d(&sum, m_in);

    if (sum > 0)
    {
        multiply_matrix_by_scalar_d(m_out, m_in, 1.0 / sum);
    } 
    else 
    {
        copy_matrix_d(m_out, m_in);
    }
}

void calc_frobenius_norm_f(float* norm_out, const Matrix_f* m)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    float    norm;

    num_rows = m->num_rows;
    num_cols = m->num_cols;
    norm     = 0;

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

    *norm_out = sqrt(norm);
}

void calc_frobenius_norm_d(double* norm_out, const Matrix_d* m)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    double   norm;

    num_rows = m->num_rows;
    num_cols = m->num_cols;
    norm     = 0;

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

    *norm_out = sqrt(norm);
}

Error* calc_frobenius_norm_diff_f
(
    float*          norm_out, 
    const Matrix_f* m_1,
    const Matrix_f* m_2
)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    float    norm;
    float    d;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows || num_cols != m_2->num_cols)
    {
        return JWSC_EARG("Matrix dimensions do not match");
    }

    norm = 0;

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            d = m_1->elts[ row ][ col ] - m_2->elts[ row ][ col ];

            norm += d * d;
        }
    }

    *norm_out = sqrt(norm);

    return NULL;
}

Error* calc_frobenius_norm_diff_d
(
    double*         norm_out, 
    const Matrix_d* m_1,
    const Matrix_d* m_2
)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    double   norm;
    double   d;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows || num_cols != m_2->num_cols)
    {
        return JWSC_EARG("Matrix dimensions do not match");
    }

    norm = 0;

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            d = m_1->elts[ row ][ col ] - m_2->elts[ row ][ col ];

            norm += d * d;
        }
    }

    *norm_out = sqrt(norm);

    return NULL;
}

Error* calc_frobenius_norm_abs_diff_f
(
    float*          norm_out, 
    const Matrix_f* m_1,
    const Matrix_f* m_2
)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    float    norm;
    float    d;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows || num_cols != m_2->num_cols)
    {
        return JWSC_EARG("Matrix dimensions do not match");
    }

    norm = 0;

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            d = fabs(m_1->elts[ row ][ col ]) - fabs(m_2->elts[ row ][ col ]);

            norm += d * d;
        }
    }

    *norm_out = sqrt(norm);

    return NULL;
}

Error* calc_frobenius_norm_abs_diff_d
(
    double*         norm_out, 
    const Matrix_d* m_1,
    const Matrix_d* m_2
)
{
    uint32_t row, num_rows;
    uint32_t col, num_cols;
    double   norm;
    double   d;

    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

    if (num_rows != m_2->num_rows || num_cols != m_2->num_cols)
    {
        return JWSC_EARG("Matrix dimensions do not match");
    }

    norm = 0;

    for (row = 0; row < num_rows; row++)
    {
        for (col = 0; col < num_cols; col++)
        {
            d = fabs(m_1->elts[ row ][ col ]) - fabs(m_2->elts[ row ][ col ]);

            norm += d * d;
        }
    }

    *norm_out = sqrt(norm);

    return NULL;
}

#if defined JWSC_HAVE_LAPACK

Error* lu_decompose_matrix_f
(
    Matrix_f**      m_out, 
    Vector_i32**    pivot_out,
    const Matrix_f* m_in
)
{
    int info;
    uint32_t pivot_len;
    uint32_t num_rows;
    uint32_t num_cols;

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

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

    pivot_len = (num_rows < num_cols) ? num_rows : num_cols;
    create_vector_i32(pivot_out, pivot_len);

    lapack_sgetrf(num_rows, num_cols, *((*m_out)->elts), num_rows, 
            (*pivot_out)->elts, &info);

    assert(info >= 0);
    if (info > 0)
    {
        if (*m_out != m_in)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        free_vector_i32(*pivot_out); *pivot_out = NULL;
        return JWSC_ECALC("Singular matrix");
    }

    return NULL;
}

Error* lu_decompose_matrix_d
(
    Matrix_d**      m_out, 
    Vector_i32**    pivot_out,
    const Matrix_d* m_in
)
{
    int info;
    uint32_t pivot_len;
    uint32_t num_rows;
    uint32_t num_cols;

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

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

    pivot_len = (num_rows < num_cols) ? num_rows : num_cols;
    create_vector_i32(pivot_out, pivot_len);

    lapack_dgetrf(num_rows, num_cols, *((*m_out)->elts), num_rows, 
            (*pivot_out)->elts, &info);

    assert(info >= 0);
    if (info > 0)
    {
        if (*m_out != m_in)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        free_vector_i32(*pivot_out); *pivot_out = NULL;
        return JWSC_ECALC("Singular matrix");
    }

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* invert_matrix_f(Matrix_f** m_out, const Matrix_f* m_in)
{
    Vector_i32* pivot = NULL;
    Error* e;

    if ((e = lu_decompose_matrix_f(m_out, &pivot, m_in)) != NULL)
    {
        return e;
    }

    if ((e = invert_lu_matrix_f(m_out, *m_out, pivot)) != NULL)
    {
        free_vector_i32(pivot);
        return e;
    }

    free_vector_i32(pivot);

    return NULL;
}

Error* invert_matrix_d(Matrix_d** m_out, const Matrix_d* m_in)
{
    Vector_i32* pivot = NULL;
    Error* e;

    if ((e = lu_decompose_matrix_d(m_out, &pivot, m_in)) != NULL)
    {
        return e;
    }

    if ((e = invert_lu_matrix_d(m_out, *m_out, pivot)) != NULL)
    {
        free_vector_i32(pivot);
        return e;
    }

    free_vector_i32(pivot);

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* invert_lu_matrix_f
(
    Matrix_f**        m_out, 
    const Matrix_f*   m_in, 
    const Vector_i32* pivot
)
{
    int      info;
    uint32_t N;
    float*   work;

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

    N = m_in->num_rows;

    if (N != m_in->num_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrix not square");
    }

    assert(work = malloc(N*sizeof(float)));

    lapack_sgetri(N, *((*m_out)->elts), N, pivot->elts, 
            work, N, &info);

    free(work);

    assert(info >= 0);
    if (info > 0)
    {
        if (*m_out != m_in)
        {
            free_matrix_f(*m_out); *m_out = NULL;
        }
        return JWSC_ECALC("Singular matrix");
    }

    return NULL;
}

Error* invert_lu_matrix_d
(
    Matrix_d**        m_out, 
    const Matrix_d*   m_in, 
    const Vector_i32* pivot
)
{
    int      info;
    uint32_t N;
    double*  work;

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

    N = m_in->num_rows;

    if (N != m_in->num_cols)
    {
        if (*m_out != m_in)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Matrix not square");
    }

    assert(work = malloc(N*sizeof(double)));

    lapack_dgetri(N, *((*m_out)->elts), N, pivot->elts, 
            work, N, &info);

    free(work);

    assert(info >= 0);
    if (info > 0)
    {
        if (*m_out != m_in)
        {
            free_matrix_d(*m_out); *m_out = NULL;
        }
        return JWSC_ECALC("Singular matrix");
    }

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* calc_matrix_determinant_f(float* det_out, const Matrix_f* m)
{
    Matrix_f*   lu    = NULL;
    Vector_i32* pivot = NULL;
    Error*      e;

    if ((e = lu_decompose_matrix_f(&lu, &pivot, m)) != NULL)
    {
        return e;
    }

    if ((e = calc_lu_matrix_determinant_f(det_out, lu, pivot)) != NULL) 
    {
        free_matrix_f(lu);
        free_vector_i32(pivot);
        return e;
    }

    free_matrix_f(lu);
    free_vector_i32(pivot);

    return NULL;
}

Error* calc_matrix_determinant_d(double* det_out, const Matrix_d* m)
{
    Matrix_d*   lu    = NULL;
    Vector_i32* pivot = NULL;
    Error*      e;

    if ((e = lu_decompose_matrix_d(&lu, &pivot, m)) != NULL)
    {
        return e;
    }

    if ((e = calc_lu_matrix_determinant_d(det_out, lu, pivot)) != NULL) 
    {
        free_matrix_d(lu);
        free_vector_i32(pivot);
        return e;
    }

    free_matrix_d(lu);
    free_vector_i32(pivot);

    return NULL;
}

#endif




Error* calc_lu_matrix_determinant_f
(
    float*            det_out,
    const Matrix_f*   m, 
    const Vector_i32* pivot
)
{
    uint32_t N;
    uint32_t n;

    N = m->num_rows;

    if (N != m->num_cols)
    {
        return JWSC_EARG("Matrix not square");
    }

    *det_out = 1.0f;

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

        if (pivot->elts[ n ] != (n+1))
        {
            *det_out *= -1.0f;
        }
    }

    return NULL;
}

Error* calc_lu_matrix_determinant_d
(
    double*           det_out,
    const Matrix_d*   m, 
    const Vector_i32* pivot
)
{
    uint32_t N;
    uint32_t n;

    N = m->num_rows;

    if (N != m->num_cols)
    {
        return JWSC_EARG("Matrix not square");
    }

    *det_out = 1.0;

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

        if (pivot->elts[ n ] != (n+1))
        {
            *det_out *= -1.0;
        }
    }

    return NULL;
}

#if defined JWSC_HAVE_LAPACK

Error* singular_value_decompose_matrix_f
(
    Matrix_f**      U_out,
    Vector_f**      S_out, 
    Matrix_f**      Vt_out,
    const Matrix_f* M
)
{
    uint32_t num_rows;
    uint32_t num_cols;

    float* work;
    int    lwork;
    int*   iwork;
    int    info;
    int    max_MN;
    int    min_MN;

    Matrix_f* M_copy = NULL;

    copy_matrix_f(&M_copy, M);

    num_rows = M->num_rows;
    num_cols = M->num_cols;

    min_MN = num_rows;
    max_MN = num_cols;
    if (min_MN > max_MN)
    {
        min_MN = num_cols;
        max_MN = num_rows;
    }

    create_matrix_f(U_out, num_rows, num_rows);
    create_vector_f(S_out, min_MN);
    create_matrix_f(Vt_out, num_cols, num_cols);

    if (max_MN >= 4*min_MN*(min_MN+1))
        lwork = 3*min_MN*min_MN + max_MN;
    else
        lwork = 3*min_MN*min_MN + 4*min_MN*(min_MN+1);

    assert(work = malloc(lwork*sizeof(float)));
    assert(iwork = malloc(8*min_MN*sizeof(int)));

    lapack_sgesdd('A', num_rows, num_cols, *(M_copy->elts), num_rows, 
            (*S_out)->elts, *((*U_out)->elts), num_rows, *((*Vt_out)->elts), 
            num_cols, work, lwork, iwork, &info);

    free(work);
    free(iwork);

    free_matrix_f(M_copy);

    assert(info >= 0);
    if (info > 0)
    {
        free_matrix_f(*U_out);  *U_out  = NULL;
        free_vector_f(*S_out);  *S_out  = NULL;
        free_matrix_f(*Vt_out); *Vt_out = NULL;
        return JWSC_ECALC("SVD calculation failed");
    }

    return NULL;
}

Error* singular_value_decompose_matrix_d
(
    Matrix_d**      U_out,
    Vector_d**      S_out, 
    Matrix_d**      Vt_out,
    const Matrix_d* M
)
{
    uint32_t num_rows;
    uint32_t num_cols;

    double* work;
    int     lwork;
    int*    iwork;
    int     info;
    int     max_MN;
    int     min_MN;

    Matrix_d* M_copy = NULL;

    copy_matrix_d(&M_copy, M);

    num_rows = M->num_rows;
    num_cols = M->num_cols;

    min_MN = num_rows;
    max_MN = num_cols;
    if (min_MN > max_MN)
    {
        min_MN = num_cols;
        max_MN = num_rows;
    }

    create_matrix_d(U_out, num_rows, num_rows);
    create_vector_d(S_out, min_MN);
    create_matrix_d(Vt_out, num_cols, num_cols);

    if (max_MN >= 4*min_MN*(min_MN+1))
        lwork = 3*min_MN*min_MN + max_MN;
    else
        lwork = 3*min_MN*min_MN + 4*min_MN*(min_MN+1);

    assert(work = malloc(lwork*sizeof(double)));
    assert(iwork = malloc(8*min_MN*sizeof(int)));

    lapack_dgesdd('A', num_rows, num_cols, *(M_copy->elts), num_rows, 
            (*S_out)->elts, *((*U_out)->elts), num_rows, *((*Vt_out)->elts), 
            num_cols, work, lwork, iwork, &info);

    free(work);
    free(iwork);

    free_matrix_d(M_copy);

    assert(info >= 0);
    if (info > 0)
    {
        free_matrix_d(*U_out);  *U_out  = NULL;
        free_vector_d(*S_out);  *S_out  = NULL;
        free_matrix_d(*Vt_out); *Vt_out = NULL;
        return JWSC_ECALC("SVD calculation failed");
    }

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* solve_real_symmetric_matrix_eigenproblem_f
(
    Matrix_f**      V_out,
    Vector_f**      D_out, 
    const Matrix_f* M
)
{
    uint32_t num_rows;
    uint32_t num_cols;

    float*  work;
    int     lwork;
    int     info;

    num_rows = M->num_rows;
    num_cols = M->num_cols;

    if (num_rows != num_cols)
    {
        return JWSC_EARG("Matrix not symmetric");
    }

    copy_matrix_f(V_out, M);
    create_vector_f(D_out, num_rows);

    lwork = 3*num_rows;

    assert(work = malloc(lwork*sizeof(float)));

    lapack_ssyev('V', 'U', num_rows, *((*V_out)->elts), num_rows, 
            (*D_out)->elts, work, lwork, &info);

    free(work);

    assert(info >= 0);
    if (info > 0)
    {
        free_matrix_f(*V_out);  *V_out  = NULL;
        free_vector_f(*D_out);  *D_out  = NULL;
        return JWSC_ECALC("Eigenproblem solver failed");
    }

    return NULL;
}

Error* solve_real_symmetric_matrix_eigenproblem_d
(
    Matrix_d**      V_out,
    Vector_d**      D_out, 
    const Matrix_d* M
)
{
    uint32_t num_rows;
    uint32_t num_cols;

    double* work;
    int     lwork;
    int     info;

    num_rows = M->num_rows;
    num_cols = M->num_cols;

    if (num_rows != num_cols)
    {
        return JWSC_EARG("Matrix not symmetric");
    }

    copy_matrix_d(V_out, M);
    create_vector_d(D_out, num_rows);

    lwork = 3*num_rows;

    assert(work = malloc(lwork*sizeof(double)));

    lapack_dsyev('V', 'U', num_rows, *((*V_out)->elts), num_rows, 
            (*D_out)->elts, work, lwork, &info);

    free(work);

    assert(info >= 0);
    if (info > 0)
    {
        free_matrix_d(*V_out);  *V_out  = NULL;
        free_vector_d(*D_out);  *D_out  = NULL;
        return JWSC_ECALC("Eigenproblem solver failed");
    }

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* pseudoinvert_matrix_no_svd_f
(
    Matrix_f**      M_out,
    const Matrix_f* M_in
)
{
    Matrix_f* Mt = NULL;
    Error*    err;

    transpose_matrix_f(&Mt, M_in);

    assert(multiply_matrices_f(M_out, Mt, M_in) == NULL);

    if ((err = invert_matrix_f(M_out, *M_out)) != NULL)
    {
        if (*M_out != M_in)
        {
            free_matrix_f(*M_out); *M_out = NULL;
        }
        free_error(err);
        return JWSC_ECALC("Pseudoinverse calculation failed");
    }

    assert(multiply_matrices_f(M_out, *M_out, Mt) == NULL);

    free_matrix_f(Mt);

    return NULL;
}

Error* pseudoinvert_matrix_no_svd_d
(
    Matrix_d**      M_out,
    const Matrix_d* M_in
)
{
    Matrix_d* Mt = NULL;
    Error*    err;

    transpose_matrix_d(&Mt, M_in);

    assert(multiply_matrices_d(M_out, Mt, M_in) == NULL);

    if ((err = invert_matrix_d(M_out, *M_out)) != NULL)
    {
        if (*M_out != M_in)
        {
            free_matrix_d(*M_out); *M_out = NULL;
        }
        free_error(err);
        return JWSC_ECALC("Pseudoinverse calculation failed");
    }

    assert(multiply_matrices_d(M_out, *M_out, Mt) == NULL);

    free_matrix_d(Mt);

    return NULL;
}

#endif




#if defined JWSC_HAVE_LAPACK

Error* pseudoinvert_matrix_f
(
    Matrix_f**      M_out,
    const Matrix_f* M_in
)
{
    uint32_t i ;
    Error*   err;

    Matrix_f* U = NULL;
    Vector_f* S = NULL;
    Matrix_f* T = NULL;
    Matrix_f* V = NULL;

    if ((err = singular_value_decompose_matrix_f(&U, &S, &V, M_in)) != NULL)
    {
        if (*M_out != M_in)
        {
            free_matrix_f(*M_out); *M_out = NULL;
        }
        free_error(err);
        return JWSC_ECALC("Pseudoinverse calculation with SVD failed");
    }

    transpose_matrix_f(&U, U);
    transpose_matrix_f(&V, V);
    create_zero_matrix_f(&T, M_in->num_cols, M_in->num_rows);

    for (i = 0; i < S->num_elts; i++)
    {
        if (S->elts[ i ] > 1.0e-16f)
        {
            T->elts[ i ][ i ] = 1.0f / S->elts[ i ];
        }
    }

    assert(multiply_matrices_f(M_out, T, U) == NULL);
    assert(multiply_matrices_f(M_out, V, *M_out) == NULL);

    free_matrix_f(U);
    free_vector_f(S);
    free_matrix_f(T);
    free_matrix_f(V);

    return NULL;
}

Error* pseudoinvert_matrix_d
(
    Matrix_d**      M_out,
    const Matrix_d* M_in
)
{
    uint32_t i ;
    Error*   err;

    Matrix_d* U = NULL;
    Vector_d* S = NULL;
    Matrix_d* T = NULL;
    Matrix_d* V = NULL;

    if ((err = singular_value_decompose_matrix_d(&U, &S, &V, M_in)) != NULL)
    {
        if (*M_out != M_in)
        {
            free_matrix_d(*M_out); *M_out = NULL;
        }
        free_error(err);
        return JWSC_ECALC("Pseudoinverse calculation with SVD failed");
    }

    transpose_matrix_d(&U, U);
    transpose_matrix_d(&V, V);
    create_zero_matrix_d(&T, M_in->num_cols, M_in->num_rows);

    for (i = 0; i < S->num_elts; i++)
    {
        if (S->elts[ i ] > 1.0e-24)
        {
            T->elts[ i ][ i ] = 1.0 / S->elts[ i ];
        }
    }

    assert(multiply_matrices_d(M_out, T, U) == NULL);
    assert(multiply_matrices_d(M_out, V, *M_out) == NULL);

    free_matrix_d(U);
    free_vector_d(S);
    free_matrix_d(T);
    free_matrix_d(V);

    return NULL;
}

#endif




void create_euler_rotation_matrix_f
(
    Matrix_f** m_out,
    float      phi,
    float      theta,
    float      psi
)
{
    float cos_phi, sin_phi;
    float cos_theta, sin_theta;
    float cos_psi, sin_psi;

    cos_phi   = (float)cos(phi);
    sin_phi   = (float)sin(phi);
    cos_theta = (float)cos(theta);
    sin_theta = (float)sin(theta);
    cos_psi   = (float)cos(psi);
    sin_psi   = (float)sin(psi);

    create_matrix_f(m_out, 3, 3);
    (*m_out)->elts[0][0] = cos_psi*cos_phi - cos_theta*sin_phi*sin_psi;
    (*m_out)->elts[0][1] = cos_psi*sin_phi + cos_theta*cos_phi*sin_psi;
    (*m_out)->elts[0][2] = sin_psi*sin_theta;
    (*m_out)->elts[1][0] = -sin_psi*cos_phi - cos_theta*sin_phi*cos_psi;
    (*m_out)->elts[1][1] = -sin_psi*sin_phi + cos_theta*cos_phi*cos_psi;
    (*m_out)->elts[1][2] = cos_psi*sin_theta;
    (*m_out)->elts[2][0] = sin_theta*sin_phi;
    (*m_out)->elts[2][1] = -sin_theta*cos_phi;
    (*m_out)->elts[2][2] = cos_theta;
}


void create_euler_rotation_matrix_d
(
    Matrix_d** m_out,
    float      phi,
    float      theta,
    float      psi
)
{
    double cos_phi, sin_phi;
    double cos_theta, sin_theta;
    double cos_psi, sin_psi;

    cos_phi   = cos(phi);
    sin_phi   = sin(phi);
    cos_theta = cos(theta);
    sin_theta = sin(theta);
    cos_psi   = cos(psi);
    sin_psi   = sin(psi);

    create_matrix_d(m_out, 3, 3);
    (*m_out)->elts[0][0] = cos_psi*cos_phi - cos_theta*sin_phi*sin_psi;
    (*m_out)->elts[0][1] = cos_psi*sin_phi + cos_theta*cos_phi*sin_psi;
    (*m_out)->elts[0][2] = sin_psi*sin_theta;
    (*m_out)->elts[1][0] = -sin_psi*cos_phi - cos_theta*sin_phi*cos_psi;
    (*m_out)->elts[1][1] = -sin_psi*sin_phi + cos_theta*cos_phi*cos_psi;
    (*m_out)->elts[1][2] = cos_psi*sin_theta;
    (*m_out)->elts[2][0] = sin_theta*sin_phi;
    (*m_out)->elts[2][1] = -sin_theta*cos_phi;
    (*m_out)->elts[2][2] = cos_theta;
}

Error* create_3d_rotation_matrix_1_f
(
    Matrix_f**      m_out, 
    float           phi,
    const Vector_f* v
)
{
    if (v->num_elts != 3)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D rotation must have three dimensions");
    }

    return create_3d_rotation_matrix_2_f(m_out, phi, 
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_rotation_matrix_1_d
(
    Matrix_d**      m_out, 
    double          phi,
    const Vector_d* v
)
{
    if (v->num_elts != 3)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D rotation must have three dimensions");
    }

    return create_3d_rotation_matrix_2_d(m_out, phi, 
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_rotation_matrix_2_f
(
    Matrix_f** m_out, 
    float      phi,
    float      x,
    float      y,
    float      z
)
{
    float cos_phi, cos_theta, cos_psi;
    float sin_phi, sin_theta, sin_psi;
    float mag_xyz, mag_xy;

    Matrix_f* R_phi   = NULL;
    Matrix_f* R_theta = NULL;
    Matrix_f* R_psi   = NULL;

    mag_xyz = sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8f)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Magnitude of rotation vector too small");
    }

    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8f)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0f;
        sin_psi = 1.0f;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = (float)cos((double)phi);
    sin_phi = (float)sin((double)phi);

    create_identity_matrix_f(&R_phi, 3);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;

    create_identity_matrix_f(&R_theta, 3);
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_identity_matrix_f(&R_psi, 3);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;

    /* *m_out = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_f(m_out, R_psi, R_theta);
    multiply_matrices_f(m_out, *m_out, R_phi);
    transpose_matrix_f(&R_theta, R_theta);
    transpose_matrix_f(&R_psi, R_psi);
    multiply_matrices_f(m_out, *m_out, R_theta);
    multiply_matrices_f(m_out, *m_out, R_psi);

    free_matrix_f(R_phi);
    free_matrix_f(R_theta);
    free_matrix_f(R_psi);

    return NULL;
}


Error* create_3d_rotation_matrix_2_d
(
    Matrix_d** m_out, 
    double     phi,
    double     x,
    double     y,
    double     z
)
{
    double cos_phi, cos_theta, cos_psi;
    double sin_phi, sin_theta, sin_psi;
    double mag_xyz, mag_xy;

    Matrix_d* R_phi   = NULL;
    Matrix_d* R_theta = NULL;
    Matrix_d* R_psi   = NULL;

    mag_xyz = sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Magnitude of rotation vector too small");
    }

    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0;
        sin_psi = 1.0;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    create_identity_matrix_d(&R_phi, 3);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;

    create_identity_matrix_d(&R_theta, 3);
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_identity_matrix_d(&R_psi, 3);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;

    /* *m_out = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_d(m_out, R_psi, R_theta);
    multiply_matrices_d(m_out, *m_out, R_phi);
    transpose_matrix_d(&R_theta, R_theta);
    transpose_matrix_d(&R_psi, R_psi);
    multiply_matrices_d(m_out, *m_out, R_theta);
    multiply_matrices_d(m_out, *m_out, R_psi);

    free_matrix_d(R_phi);
    free_matrix_d(R_theta);
    free_matrix_d(R_psi);

    return NULL;
}

Error* create_3d_homo_rotation_matrix_1_f
(
    Matrix_f**      m_out, 
    float           phi,
    const Vector_f* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous rotation is not 3D");
    }

    return create_3d_homo_rotation_matrix_2_f(m_out, phi, 
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_rotation_matrix_1_d
(
    Matrix_d**      m_out, 
    double          phi,
    const Vector_d* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous rotation is not 3D");
    }

    return create_3d_homo_rotation_matrix_2_d(m_out, phi, 
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_rotation_matrix_2_f
(
    Matrix_f** m_out, 
    float      phi,
    float      x,
    float      y,
    float      z
)
{
    float cos_phi, cos_theta, cos_psi;
    float sin_phi, sin_theta, sin_psi;
    float mag_xyz, mag_xy;

    Matrix_f* R_phi   = NULL;
    Matrix_f* R_theta = NULL;
    Matrix_f* R_psi   = NULL;

    mag_xyz = sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8f)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Magnitude of rotation vector too small");
    }

    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8f)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0f;
        sin_psi = 1.0f;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = (float)cos((double)phi);
    sin_phi = (float)sin((double)phi);

    create_identity_matrix_f(&R_phi, 4);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;

    create_identity_matrix_f(&R_theta, 4);
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_identity_matrix_f(&R_psi, 4);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;

    /* *m_out = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_f(m_out, R_psi, R_theta);
    multiply_matrices_f(m_out, *m_out, R_phi);
    transpose_matrix_f(&R_theta, R_theta);
    transpose_matrix_f(&R_psi, R_psi);
    multiply_matrices_f(m_out, *m_out, R_theta);
    multiply_matrices_f(m_out, *m_out, R_psi);

    free_matrix_f(R_phi);
    free_matrix_f(R_theta);
    free_matrix_f(R_psi);

    return NULL;
}


Error* create_3d_homo_rotation_matrix_2_d
(
    Matrix_d** m_out, 
    double     phi,
    double     x,
    double     y,
    double     z
)
{
    double cos_phi, cos_theta, cos_psi;
    double sin_phi, sin_theta, sin_psi;
    double mag_xyz, mag_xy;

    Matrix_d* R_phi   = NULL;
    Matrix_d* R_theta = NULL;
    Matrix_d* R_psi   = NULL;

    mag_xyz = sqrt(x*x + y*y + z*z);

    if (mag_xyz < 1.0e-8)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Magnitude of rotation vector too small");
    }

    /* psi used to rotate around z-axis to put vector on the yz-plane. */
    mag_xy = sqrt(x*x + y*y);
    if (mag_xy > 1.0e-8)
    {
        cos_psi = y / mag_xy;
        sin_psi = x / mag_xy;
    }
    else
    {
        cos_psi = 0.0;
        sin_psi = 1.0;
    }

    /* theta used to rotate around x-axis to put vector onto the z-axis. */
    cos_theta = z / mag_xyz;
    sin_theta = mag_xy / mag_xyz;

    /* phi used to rotate around z-axis. This is actually the angle to rotate
     * around the vector. */
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    create_identity_matrix_d(&R_phi, 4);
    R_phi->elts[ 0 ][ 0 ] = cos_phi;
    R_phi->elts[ 0 ][ 1 ] = sin_phi;
    R_phi->elts[ 1 ][ 0 ] = -sin_phi;
    R_phi->elts[ 1 ][ 1 ] = cos_phi;

    create_identity_matrix_d(&R_theta, 4);
    R_theta->elts[ 1 ][ 1 ] = cos_theta;
    R_theta->elts[ 1 ][ 2 ] = sin_theta;
    R_theta->elts[ 2 ][ 1 ] = -sin_theta;
    R_theta->elts[ 2 ][ 2 ] = cos_theta;

    create_identity_matrix_d(&R_psi, 4);
    R_psi->elts[ 0 ][ 0 ] = cos_psi;
    R_psi->elts[ 0 ][ 1 ] = sin_psi;
    R_psi->elts[ 1 ][ 0 ] = -sin_psi;
    R_psi->elts[ 1 ][ 1 ] = cos_psi;

    /* *m_out = R_psi * R_theta * R_phi * R_theta^T * R_psi^T */
    multiply_matrices_d(m_out, R_psi, R_theta);
    multiply_matrices_d(m_out, *m_out, R_phi);
    transpose_matrix_d(&R_theta, R_theta);
    transpose_matrix_d(&R_psi, R_psi);
    multiply_matrices_d(m_out, *m_out, R_theta);
    multiply_matrices_d(m_out, *m_out, R_psi);

    free_matrix_d(R_phi);
    free_matrix_d(R_theta);
    free_matrix_d(R_psi);

    return NULL;
}

Error* create_3d_scaling_matrix_1_f
(
    Matrix_f**      m_out, 
    const Vector_f* v
)
{
    if (v->num_elts != 3)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D scaling is not 3D");
    }

    return create_3d_scaling_matrix_2_f(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_scaling_matrix_1_d
(
    Matrix_d**      m_out, 
    const Vector_d* v
)
{
    if (v->num_elts != 3)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D scaling is not 3D");
    }

    return create_3d_scaling_matrix_2_d(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_scaling_matrix_2_f
(
    Matrix_f** m_out, 
    float      x,
    float      y,
    float      z
)
{
    create_zero_matrix_f(m_out, 3, 3);
    (*m_out)->elts[ 0 ][ 0 ] = x;
    (*m_out)->elts[ 1 ][ 1 ] = y;
    (*m_out)->elts[ 2 ][ 2 ] = z;

    return NULL;
}


Error* create_3d_scaling_matrix_2_d
(
    Matrix_d** m_out, 
    double     x,
    double     y,
    double     z
)
{
    create_zero_matrix_d(m_out, 3, 3);
    (*m_out)->elts[ 0 ][ 0 ] = x;
    (*m_out)->elts[ 1 ][ 1 ] = y;
    (*m_out)->elts[ 2 ][ 2 ] = z;

    return NULL;
}

Error* create_3d_homo_scaling_matrix_1_f
(
    Matrix_f**      m_out, 
    const Vector_f* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous scaling is not 3D");
    }

    return create_3d_homo_scaling_matrix_2_f(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_scaling_matrix_1_d
(
    Matrix_d**      m_out, 
    const Vector_d* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous scaling is not 3D");
    }

    return create_3d_homo_scaling_matrix_2_d(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_scaling_matrix_2_f
(
    Matrix_f** m_out, 
    float      x,
    float      y,
    float      z
)
{
    create_identity_matrix_f(m_out, 4);
    (*m_out)->elts[ 0 ][ 0 ] = x;
    (*m_out)->elts[ 1 ][ 1 ] = y;
    (*m_out)->elts[ 2 ][ 2 ] = z;

    return NULL;
}


Error* create_3d_homo_scaling_matrix_2_d
(
    Matrix_d** m_out, 
    double     x,
    double     y,
    double     z
)
{
    create_identity_matrix_d(m_out, 4);
    (*m_out)->elts[ 0 ][ 0 ] = x;
    (*m_out)->elts[ 1 ][ 1 ] = y;
    (*m_out)->elts[ 2 ][ 2 ] = z;

    return NULL;
}

Error* create_3d_homo_translation_matrix_1_f
(
    Matrix_f**      m_out, 
    const Vector_f* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_f(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous translation is not 3D");
    }

    return create_3d_homo_translation_matrix_2_f(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_translation_matrix_1_d
(
    Matrix_d**      m_out, 
    const Vector_d* v
)
{
    if (v->num_elts < 3)
    {
        free_matrix_d(*m_out); *m_out = NULL;
        return JWSC_EARG("Vector for 3D homogeneous translation is not 3D");
    }

    return create_3d_homo_translation_matrix_2_d(m_out,
            v->elts[0], v->elts[1], v->elts[2]);
}


Error* create_3d_homo_translation_matrix_2_f
(
    Matrix_f** m_out, 
    float      x,
    float      y,
    float      z
)
{
    create_identity_matrix_f(m_out, 4);
    (*m_out)->elts[ 0 ][ 3 ] = x;
    (*m_out)->elts[ 1 ][ 3 ] = y;
    (*m_out)->elts[ 2 ][ 3 ] = z;

    return NULL;
}


Error* create_3d_homo_translation_matrix_2_d
(
    Matrix_d** m_out, 
    double     x,
    double     y,
    double     z
)
{
    create_identity_matrix_d(m_out, 4);
    (*m_out)->elts[ 0 ][ 3 ] = x;
    (*m_out)->elts[ 1 ][ 3 ] = y;
    (*m_out)->elts[ 2 ][ 3 ] = z;

    return NULL;
}