matblock_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/blas.h"
#include "jwsc/matblock/matblock.h"
#include "jwsc/matblock/matblock_math.h"


void add_scalar_to_matblock_f
(
    Matblock_f**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matblock_f(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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

void add_scalar_to_matblock_d
(
    Matblock_d**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matblock_d(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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

void add_scalar_to_matblock_cf
(
    Matblock_cf**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matblock_cf(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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

void add_scalar_to_matblock_cd
(
    Matblock_cd**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

    if (*m_out != m_in)
    {
        create_matblock_cd(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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

void multiply_matblock_by_scalar_f
(
    Matblock_f**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

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

    blas_sscal(num_elts, a, **((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matblock_f(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }

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

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

void multiply_matblock_by_scalar_d
(
    Matblock_d**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

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

    blas_dscal(num_elts, a, **((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matblock_d(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }

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

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

void multiply_matblock_by_scalar_cf
(
    Matblock_cf**      m_out, 
    const Matblock_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_mats*m_in->num_rows*m_in->num_cols;

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

    blas_cscal(num_elts, z, **((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matblock_cf(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }

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

    for (elt = 0; elt < num_elts; elt++)
    {
        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_matblock_by_scalar_cd
(
    Matblock_cd**      m_out, 
    const Matblock_cd* m_in, 
    Complex_d          z
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t   elt;
    double     r, i;
    Complex_d* m_out_elts;
    Complex_d* m_in_elts;
#endif

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

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

    blas_zscal(num_elts, z, **((*m_out)->elts), 1);
#else
    if (*m_out != m_in)
    {
        create_matblock_cd(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }

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

    for (elt = 0; elt < num_elts; elt++)
    {
        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_matblocks_elt_wise_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_1, 
    const Matblock_f* m_2
)
{
    uint32_t  num_mats, num_rows, num_cols;
    uint32_t  num_elts;
    uint32_t  elt;
    float*    m_out_elts;
    float*    m_1_elts;
    float*    m_2_elts;

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

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

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

    num_elts = num_mats*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_matblocks_elt_wise_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_1, 
    const Matblock_d* m_2
)
{
    uint32_t  num_mats, num_rows, num_cols;
    uint32_t  num_elts;
    uint32_t  elt;
    double*   m_out_elts;
    double*   m_1_elts;
    double*   m_2_elts;

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

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

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

    num_elts = num_mats*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_matblocks_elt_wise_cf
(
    Matblock_cf**      m_out, 
    const Matblock_cf* m_1, 
    const Matblock_cf* m_2
)
{
    uint32_t   num_mats, 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_mats = m_1->num_mats;
    num_rows = m_1->num_rows;
    num_cols = m_1->num_cols;

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

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

    num_elts = num_mats*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_matblocks_elt_wise_cd
(
    Matblock_cd**      m_out, 
    const Matblock_cd* m_1, 
    const Matblock_cd* m_2
)
{
    uint32_t   num_mats, num_rows, num_cols;
    uint32_t   num_elts;
    uint32_t   elt;
    double     r, i;
    Complex_d* m_out_elts;
    Complex_d* m_1_elts;
    Complex_d* m_2_elts;

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

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

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

    num_elts = num_mats*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* add_matblocks_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_1, 
    const Matblock_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_mats != m_2->num_mats ||
        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_matblock_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matblocks requires equal dimensions");
    }

    num_elts = m_1->num_mats*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_matblock_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_matblock_f(m_out, m_1->num_cols, 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_matblocks_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_1, 
    const Matblock_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_mats != m_2->num_mats ||
        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_matblock_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Adding matblocks requires equal dimensions");
    }

    num_elts = m_1->num_mats*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_matblock_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_matblock_d(m_out, m_1->num_cols, 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_matblocks_f
(
    Matblock_f**      m_out, 
    const Matblock_f* m_1, 
    const Matblock_f* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matblock_f* m = NULL;
#else
    uint32_t elt;
    float*   m_1_elts;
    float*   m_2_elts;
    float*   m_out_elts;
#endif

    if (m_1->num_mats != m_2->num_mats ||
        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_matblock_f(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matblocks requires equal dimensions");
    }

    num_elts = m_1->num_mats*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_matblock_f(&m, m_1);
        blas_saxpy(num_elts, -1, **(m_2->elts), 1, **(m->elts), 1);
        copy_matblock_f(m_out, m);
        free_matblock_f(m);
    }
    else
    {
        copy_matblock_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_matblock_f(m_out, m_1->num_mats, 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_matblocks_d
(
    Matblock_d**      m_out, 
    const Matblock_d* m_1, 
    const Matblock_d* m_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Matblock_d* m = NULL;
#else
    uint32_t elt;
    double*  m_1_elts;
    double*  m_2_elts;
    double*  m_out_elts;
#endif

    if (m_1->num_mats != m_2->num_mats ||
        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_matblock_d(*m_out); *m_out = NULL;
        }
        return JWSC_EARG("Subtracting matblocks requires equal dimensions");
    }

    num_elts = m_1->num_mats*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_matblock_d(&m, m_1);
        blas_daxpy(num_elts, -1, **(m_2->elts), 1, **(m->elts), 1);
        copy_matblock_d(m_out, m);
        free_matblock_d(m);
    }
    else
    {
        copy_matblock_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_matblock_d(m_out, m_1->num_mats, 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;
}

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

    num_elts = m->num_mats*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_matblock_asum_d(double* sum_out, const Matblock_d* m)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *sum_out = blas_dasum(m->num_mats*m->num_rows*m->num_cols, **(m->elts), 1);
#else
    uint32_t num_elts;
    uint32_t elt;
    double*  m_elts;

    num_elts = m->num_mats*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_matblock_sum_f(Matblock_f** m_out, const Matblock_f* m_in)
{
    float sum;

    calc_matblock_asum_f(&sum, m_in);

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

void normalize_matblock_sum_d(Matblock_d** m_out, const Matblock_d* m_in)
{
    double sum;

    calc_matblock_asum_d(&sum, m_in);

    if (sum > 0)
    {
        multiply_matblock_by_scalar_d(m_out, m_in, 1.0 / sum);
    } 
    else 
    {
        copy_matblock_d(m_out, m_in);
    }
}

void abs_matblock_f(Matblock_f** m_out, const Matblock_f* m_in)
{
    uint32_t  num_elts;
    uint32_t  elt;
    float*    m_in_elts;
    float*    m_out_elts;

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

    if (*m_out != m_in)
    {
        create_matblock_f(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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


void abs_matblock_d(Matblock_d** m_out, const Matblock_d* m_in)
{
    uint32_t  num_elts;
    uint32_t  elt;
    double*    m_in_elts;
    double*    m_out_elts;

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

    if (*m_out != m_in)
    {
        create_matblock_d(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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


void threshold_matblock_f(Matblock_f** m_out, const Matblock_f* m_in, float threshold)
{
    uint32_t  num_elts;
    uint32_t  elt;
    float*    m_in_elts;
    float*    m_out_elts;

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

    if (*m_out != m_in)
    {
        create_matblock_f(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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

void threshold_matblock_d(Matblock_d** m_out, const Matblock_d* m_in, double threshold)
{
    uint32_t  num_elts;
    uint32_t  elt;
    double *  m_in_elts;
    double *  m_out_elts;

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

    if (*m_out != m_in)
    {
        create_matblock_d(m_out, m_in->num_mats, m_in->num_rows, 
                m_in->num_cols);
    }
    m_out_elts = **((*m_out)->elts);
    m_in_elts  = **(m_in->elts);

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