matblock_fft.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 <inttypes.h>
#include <assert.h>
#include <math.h>
#ifdef JWSC_HAVE_FFTW3_H
#include <fftw3.h>
#endif

#include "jwsc/matblock/matblock.h"
#include "jwsc/matblock/matblock_math.h"
#include "jwsc/matblock/matblock_fft.h"


#ifdef JWSC_HAVE_FFTW3F
void real_to_complex_matblock_dft_f
(
    Matblock_cf**     m_out,
    const Matblock_f* m_in
)
{
    Matblock_f* s = NULL;

    real_to_complex_matblock_dft_scratch_f(m_out, m_in, &s);

    free_matblock_f(s);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void real_to_complex_matblock_dft_d
(
    Matblock_cd**     m_out,
    const Matblock_d* m_in
)
{
    Matblock_d* s = NULL;

    real_to_complex_matblock_dft_scratch_d(m_out, m_in, &s);

    free_matblock_d(s);
}
#endif

#ifdef JWSC_HAVE_FFTW3F
void real_to_complex_matblock_dft_scratch_f
(
    Matblock_cf**     m_out,
    const Matblock_f* m_in,
    Matblock_f**      s
)
{
    uint32_t    num_mats, num_rows, num_cols;
    fftwf_plan  plan;

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

    create_matblock_cf(m_out, num_mats, num_rows, (num_cols/2)+1);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    create_matblock_f(s, num_mats, num_rows, num_cols);

    assert((plan = fftwf_plan_dft_r2c_3d(num_mats, num_rows, num_cols,
                    **((*s)->elts), (fftwf_complex*)**((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matblock_f(s, m_in);
    fftwf_execute(plan);

    fftwf_destroy_plan(plan);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void real_to_complex_matblock_dft_scratch_d
(
    Matblock_cd**     m_out,
    const Matblock_d* m_in,
    Matblock_d**      s
)
{
    uint32_t    num_mats, num_rows, num_cols;
    fftw_plan   plan;

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

    create_matblock_cd(m_out, num_mats, num_rows, (num_cols/2)+1);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    create_matblock_d(s, num_mats, num_rows, num_cols);

    assert((plan = fftw_plan_dft_r2c_3d(num_mats, num_rows, num_cols,
                    **((*s)->elts), (fftw_complex*)**((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matblock_d(s, m_in);
    fftw_execute(plan);

    fftw_destroy_plan(plan);
}
#endif

#ifdef JWSC_HAVE_FFTW3F
void complex_to_real_matblock_dft_f
(
    Matblock_f**       m_out,
    const Matblock_cf* m_in,
    uint8_t            ncols_even
)
{
    Matblock_cf* s = NULL;

    complex_to_real_matblock_dft_scratch_f(m_out, m_in, ncols_even, &s);

    free_matblock_cf(s);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void complex_to_real_matblock_dft_d
(
    Matblock_d**       m_out,
    const Matblock_cd* m_in,
    uint8_t            ncols_even
)
{
    Matblock_cd* s = NULL;;

    complex_to_real_matblock_dft_scratch_d(m_out, m_in, ncols_even, &s);

    free_matblock_cd(s);
}
#endif

#ifdef JWSC_HAVE_FFTW3F
void complex_to_real_matblock_dft_scratch_f
(
    Matblock_f**       m_out,
    const Matblock_cf* m_in,
    uint8_t            ncols_even,
    Matblock_cf**      s
)
{
    uint32_t     num_mats, num_rows, num_cols;
    fftwf_plan   plan;

    num_mats = m_in->num_mats;
    num_rows = m_in->num_rows;
    num_cols = (ncols_even) ? (m_in->num_cols-1)*2 : (m_in->num_cols-1)*2 + 1;

    create_matblock_f(m_out, num_mats, num_rows, num_cols);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    create_matblock_cf(s, m_in->num_mats, m_in->num_rows, m_in->num_cols);

    assert((plan = fftwf_plan_dft_c2r_3d(num_mats, num_rows, num_cols,
                    (fftwf_complex*)**((*s)->elts), **((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matblock_cf(s, m_in);
    fftwf_execute(plan);

    fftwf_destroy_plan(plan);
}
#endif

#ifdef JWSC_HAVE_FFTW3
void complex_to_real_matblock_dft_scratch_d
(
    Matblock_d**       m_out,
    const Matblock_cd* m_in,
    uint8_t            ncols_even,
    Matblock_cd**      s
)
{
    uint32_t     num_mats, num_rows, num_cols;
    fftw_plan    plan;

    num_mats = m_in->num_mats;
    num_rows = m_in->num_rows;
    num_cols = (ncols_even) ? (m_in->num_cols-1)*2 : (m_in->num_cols-1)*2 + 1;

    create_matblock_d(m_out, num_mats, num_rows, num_cols);

    /* For fftw planning, input array does not need to be initialized, hence we
     * allow fftw to destroy the input to (possibly) get better efficiency. */
    create_matblock_cd(s, m_in->num_mats, m_in->num_rows, m_in->num_cols);

    assert((plan = fftw_plan_dft_c2r_3d(num_mats, num_rows, num_cols,
                    (fftw_complex*)**((*s)->elts), **((*m_out)->elts), 
                    FFTW_MEASURE | FFTW_DESTROY_INPUT)) != NULL);

    /* Copy the input, now that planning is over, and execute. */
    copy_matblock_cd(s, m_in);
    fftw_execute(plan);

    fftw_destroy_plan(plan);
}
#endif

void expand_symmetric_matblock_dft_f
(
    Matblock_cf**      m_out, 
    const Matblock_cf* m_in, 
    uint8_t            ncols_even
)
{
    uint32_t     num_mats, num_rows, num_cols;
    uint32_t     mat, row, col;
    Matblock_cf* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    if (ncols_even)
    {
        create_matblock_cf(&m, num_mats, num_rows, (m_in->num_cols-1)*2);
    }
    else
    {
        create_matblock_cf(&m, num_mats, num_rows, (m_in->num_cols-1)*2 + 1);
    }

    assert(copy_matblock_block_into_matblock_cf(m, 0, 0, 0, m_in, 0, 0, 0, 
                m_in->num_mats, m_in->num_rows, m_in->num_cols) == NULL);

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

    if (*m_out == m_in)
    {
        copy_matblock_cf(m_out, m);
        free_matblock_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void expand_symmetric_matblock_dft_d
(
    Matblock_cd**      m_out, 
    const Matblock_cd* m_in, 
    uint8_t            ncols_even
)
{
    uint32_t     num_mats, num_rows, num_cols;
    uint32_t     mat, row, col;
    Matblock_cd* m;

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

    m = (*m_out == m_in) ? NULL : *m_out;
    if (ncols_even)
    {
        create_matblock_cd(&m, num_mats, num_rows, (m_in->num_cols-1)*2);
    }
    else
    {
        create_matblock_cd(&m, num_mats, num_rows, (m_in->num_cols-1)*2 + 1);
    }

    assert(copy_matblock_block_into_matblock_cd(m, 0, 0, 0, m_in, 0, 0, 0, 
                m_in->num_mats, m_in->num_rows, m_in->num_cols) == NULL);

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

    if (*m_out == m_in)
    {
        copy_matblock_cd(m_out, m);
        free_matblock_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void shrink_symmetric_matblock_dft_f
(
    Matblock_cf**      m_out, 
    const Matblock_cf* m_in
)
{
    Matblock_cf* m;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matblock_cf(&m, m_in->num_mats, m_in->num_rows, 
            (m_in->num_cols/2)+1);

    assert(copy_matblock_block_into_matblock_cf(m, 0, 0, 0, m_in, 0, 0, 0, 
                m->num_mats, m->num_rows, m->num_cols) == NULL);

    if (*m_out == m_in)
    {
        copy_matblock_cf(m_out, m);
        free_matblock_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void shrink_symmetric_matblock_dft_d
(
    Matblock_cd**      m_out, 
    const Matblock_cd* m_in
)
{
    Matblock_cd* m;

    m = (*m_out == m_in) ? NULL : *m_out;
    create_matblock_cd(&m, m_in->num_mats, m_in->num_rows, 
            (m_in->num_cols/2)+1);

    assert(copy_matblock_block_into_matblock_cd(m, 0, 0, 0, m_in, 0, 0, 0, 
                m->num_mats, m->num_rows, m->num_cols) == NULL);

    if (*m_out == m_in)
    {
        copy_matblock_cd(m_out, m);
        free_matblock_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void scale_matblock_dft_f(Matblock_f** m_out, const Matblock_f* m_in)
{
    float scale;

    scale = 1.0f / (float)(m_in->num_mats*m_in->num_rows*m_in->num_cols);
    multiply_matblock_by_scalar_f(m_out, m_in, scale);
}

void scale_matblock_dft_d(Matblock_d** m_out, const Matblock_d* m_in)
{
    float scale;

    scale = 1.0 / (double)(m_in->num_mats*m_in->num_rows*m_in->num_cols);
    multiply_matblock_by_scalar_d(m_out, m_in, scale);
}

void shift_matblock_f
(
    Matblock_f**       m_out,
    const Matblock_f*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_f* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats + (num_mats % 2))
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_f(m_out, m);
        free_matblock_f(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matblock_d
(
    Matblock_d**       m_out,
    const Matblock_d*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_d* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats + (num_mats % 2))
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_d(m_out, m);
        free_matblock_d(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matblock_cf
(
    Matblock_cf**       m_out,
    const Matblock_cf*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_cf* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats + (num_mats % 2))
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_cf(m_out, m);
        free_matblock_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void shift_matblock_cd
(
    Matblock_cd**       m_out,
    const Matblock_cd*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_cd* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats + (num_mats % 2))
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows + (num_rows % 2))
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row + half_rows ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols + (num_cols % 2))
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col + half_cols ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats - (num_mats % 2) ]
                                [ row - half_rows - (num_rows % 2) ]
                                [ col - half_cols - (num_cols % 2) ];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_cd(m_out, m);
        free_matblock_cd(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matblock_f
(
    Matblock_f**       m_out,
    const Matblock_f*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_f* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats)
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_f(m_out, m);
        free_matblock_f(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matblock_d
(
    Matblock_d**       m_out,
    const Matblock_d*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_d* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats)
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_mats ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_d(m_out, m);
        free_matblock_d(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matblock_cf
(
    Matblock_cf**       m_out,
    const Matblock_cf*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_cf* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats)
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_cf(m_out, m);
        free_matblock_cf(m);
    }
    else
    {
        *m_out = m;
    }
}

void unshift_matblock_cd
(
    Matblock_cd**       m_out,
    const Matblock_cd*  m_in
)
{
    uint32_t num_mats, num_rows, num_cols;
    uint32_t mat, row, col;
    uint32_t half_mats, half_rows, half_cols;

    Matblock_cd* m = NULL;

    num_mats  = m_in->num_mats;
    num_rows  = m_in->num_rows;
    num_cols  = m_in->num_cols;
    half_mats = num_mats / 2;
    half_rows = num_rows / 2;
    half_cols = num_cols / 2;

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

    for (mat = 0; mat < num_mats; mat++)
    {
        if (mat < half_mats)
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat + half_mats + (num_mats % 2) ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
        else
        {
            for (row = 0; row < num_rows; row++)
            {
                if (row < half_rows)
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row + half_rows + (num_rows % 2) ]
                                [ col - half_cols ];
                        }
                    }
                }
                else
                {
                    for (col = 0; col < num_cols; col++)
                    {
                        if (col < half_cols)
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row - half_rows ]
                                [ col + half_cols + (num_cols % 2) ];
                        }
                        else
                        {
                            m->elts[ mat ][ row ][ col ] = m_in->elts
                                [ mat - half_cols ]
                                [ row - half_rows]
                                [ col - half_cols];
                        }
                    }
                }
            }
        }
    }

    if (*m_out == m_in)
    {
        copy_matblock_cd(m_out, m);
        free_matblock_cd(m);
    }
    else
    {
        *m_out = m;
    }
}