vector_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 <math.h>
#include <assert.h>
#include <inttypes.h>

#include "jwsc/base/error.h"
#include "jwsc/math/complex.h"
#include "jwsc/math/blas.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector_math.h"


void add_scalar_to_vector_f(Vector_f** v_out, const Vector_f* v_in, float a)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Vector_f*  v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_f(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ] = a + v_in->elts[ elt ];
    }
}

void add_scalar_to_vector_d(Vector_d** v_out, const Vector_d* v_in, double a)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Vector_d*  v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_d(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ] = a + v_in->elts[ elt ];
    }
}

void add_scalar_to_vector_cf
(
    Vector_cf**      v_out, 
    const Vector_cf* v_in, 
    Complex_f        z
)
{
    uint32_t    num_elts;
    uint32_t    elt;
    Vector_cf*  v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_cf(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ].r = z.r + v_in->elts[ elt ].r;
        v->elts[ elt ].i = z.i + v_in->elts[ elt ].i;
    }
}

void add_scalar_to_vector_cd
(
    Vector_cd**      v_out, 
    const Vector_cd* v_in, 
    Complex_d        z
)
{
    uint32_t   num_elts;
    uint32_t   elt;
    Vector_cd*  v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_cd(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ].r = z.r + v_in->elts[ elt ].r;
        v->elts[ elt ].i = z.i + v_in->elts[ elt ].i;
    }
}

void multiply_vector_by_scalar_f
(
    Vector_f**      v_out, 
    const Vector_f* v_in, 
    float           a
)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    uint32_t  num_elts;
    Vector_f* v;

    num_elts = v_in->num_elts;

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

    blas_sscal(num_elts, a, v->elts, 1);
#else
    uint32_t  num_elts;
    uint32_t  elt;
    Vector_f* v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_f(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ] = a * v_in->elts[ elt ];
    }
#endif
}

void multiply_vector_by_scalar_d
(
    Vector_d**      v_out, 
    const Vector_d* v_in, 
    double          a
)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    uint32_t  num_elts;
    Vector_d* v;

    num_elts = v_in->num_elts;

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

    blas_dscal(num_elts, a, v->elts, 1);
#else
    uint32_t  num_elts;
    uint32_t  elt;
    Vector_d* v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_d(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        v->elts[ elt ] = a * v_in->elts[ elt ];
    }
#endif
}

void multiply_vector_by_scalar_cf
(
    Vector_cf**      v_out, 
    const Vector_cf* v_in, 
    Complex_f        z
)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    uint32_t   num_elts;
    Vector_cf* v;

    num_elts = v_in->num_elts;

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

    blas_cscal(num_elts, z, v->elts, 1);
#else
    uint32_t   num_elts;
    uint32_t   elt;
    float      r, i;
    Vector_cf* v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_cf(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = z.r*v_in->elts[ elt ].r - z.i*v_in->elts[ elt ].i;
        i = z.i*v_in->elts[ elt ].r + z.r*v_in->elts[ elt ].i;

        v->elts[ elt ].r = r;
        v->elts[ elt ].i = i;
    }
#endif
}

void multiply_vector_by_scalar_cd
(
    Vector_cd**      v_out, 
    const Vector_cd* v_in, 
    Complex_d        z
)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    uint32_t   num_elts;
    Vector_cd* v;

    num_elts = v_in->num_elts;

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

    blas_zscal(num_elts, z, v->elts, 1);
#else
    uint32_t   num_elts;
    uint32_t   elt;
    double     r, i;
    Vector_cd* v;

    num_elts = v_in->num_elts;

    if (*v_out != v_in)
    {
        create_vector_cd(v_out, num_elts);
    }
    v = *v_out;

    for (elt = 0; elt < num_elts; elt++)
    {
        r = z.r*v_in->elts[ elt ].r - z.i*v_in->elts[ elt ].i;
        i = z.i*v_in->elts[ elt ].r + z.r*v_in->elts[ elt ].i;

        v->elts[ elt ].r = r;
        v->elts[ elt ].i = i;
    }
#endif
}

Error* calc_vector_dot_prod_f
(
    float*          dp_out,
    const Vector_f* v_1, 
    const Vector_f* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    num_elts = v_1->num_elts;
    if (num_elts != v_2->num_elts)
    {
        return JWSC_EARG("Dot product requires vectors with equal dimensions");
    }

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *dp_out = blas_sdot(num_elts, v_1->elts, 1, v_2->elts, 1);
#else
    *dp_out = 0;
    for (elt = 0; elt < num_elts; elt++)
    {
        *dp_out += v_1->elts[ elt ] * v_2->elts[ elt ];
    }
#endif

    return NULL;
}

Error* calc_vector_dot_prod_d
(
    double*         dp_out,
    const Vector_d* v_1, 
    const Vector_d* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    num_elts = v_1->num_elts;
    if (num_elts != v_2->num_elts)
    {
        return JWSC_EARG("Dot product requires vectors with equal dimensions");
    }

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *dp_out = blas_ddot(num_elts, v_1->elts, 1, v_2->elts, 1);
#else
    *dp_out = 0;
    for (elt = 0; elt < num_elts; elt++)
    {
        *dp_out += v_1->elts[ elt ] * v_2->elts[ elt ];
    }
#endif

    return NULL;
}

Error* calc_3d_vector_cross_prod_f
(
    Vector_f**      v_out, 
    const Vector_f* v_1, 
    const Vector_f* v_2
)
{
    float cp_1, cp_2, cp_3;

    if (v_1->num_elts != 3 || v_2->num_elts != 3)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_f(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Cross product requires vectors with three-dimensions");
    }

    cp_1  = v_1->elts[1] * v_2->elts[2];
    cp_1 -= v_1->elts[2] * v_2->elts[1];
    cp_2  = v_1->elts[2] * v_2->elts[0];
    cp_2 -= v_1->elts[0] * v_2->elts[2];
    cp_3  = v_1->elts[0] * v_2->elts[1];
    cp_3 -= v_1->elts[1] * v_2->elts[0];

    create_vector_f(v_out, 3);
    (*v_out)->elts[0] = cp_1;
    (*v_out)->elts[1] = cp_2;
    (*v_out)->elts[2] = cp_3;

    return NULL;
}

Error* calc_3d_vector_cross_prod_d
(
    Vector_d**      v_out, 
    const Vector_d* v_1, 
    const Vector_d* v_2
)
{
    double cp_1, cp_2, cp_3;

    if (v_1->num_elts != 3 || v_2->num_elts != 3)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_d(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Cross product requires vectors with three-dimensions");
    }

    cp_1  = v_1->elts[1] * v_2->elts[2];
    cp_1 -= v_1->elts[2] * v_2->elts[1];
    cp_2  = v_1->elts[2] * v_2->elts[0];
    cp_2 -= v_1->elts[0] * v_2->elts[2];
    cp_3  = v_1->elts[0] * v_2->elts[1];
    cp_3 -= v_1->elts[1] * v_2->elts[0];

    create_vector_d(v_out, 3);
    (*v_out)->elts[0] = cp_1;
    (*v_out)->elts[1] = cp_2;
    (*v_out)->elts[2] = cp_3;

    return NULL;
}

void calc_vector_asum_f(float* sum_out, const Vector_f* v)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *sum_out = blas_sasum(v->num_elts, v->elts, 1);
#else
    uint32_t num_elts;
    uint32_t elt;

    num_elts = v->num_elts;
    *sum_out = 0;

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

void calc_vector_asum_d(double* sum_out, const Vector_d* v)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *sum_out = blas_dasum(v->num_elts, v->elts, 1);
#else
    uint32_t num_elts;
    uint32_t elt;

    num_elts = v->num_elts;
    *sum_out = 0;

    for (elt = 0; elt < num_elts; elt++)
    {
        (*sum_out) += fabs(v->elts[ elt ]);
    }
#endif
}

void calc_vector_mag_f(float* mag_out, const Vector_f* v)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *mag_out = blas_snrm2(v->num_elts, v->elts, 1);
#else
    uint32_t num_elts;
    uint32_t elt;

    num_elts = v->num_elts;
    *mag_out = 0;

    for (elt = 0; elt < num_elts; elt++)
    {
        *mag_out += v->elts[ elt ] * v->elts[ elt ];
    }

    *mag_out = (float)sqrt(*mag_out);
#endif
}

void calc_vector_mag_d(double* mag_out, const Vector_d* v)
{
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    *mag_out = blas_dnrm2(v->num_elts, v->elts, 1);
#else
    uint32_t num_elts;
    uint32_t elt;

    num_elts = v->num_elts;
    *mag_out = 0;

    for (elt = 0; elt < num_elts; elt++)
    {
        *mag_out += v->elts[ elt ] * v->elts[ elt ];
    }

    *mag_out = sqrt(*mag_out);
#endif
}

void normalize_vector_sum_f(Vector_f** v_out, const Vector_f* v_in)
{
    float sum;

    calc_vector_asum_f(&sum, v_in);

    if (sum > 0)
    {
        multiply_vector_by_scalar_f(v_out, v_in, 1.0f / sum);
    } 
    else if (v_in != *v_out)
    {
        copy_vector_f(v_out, v_in);
    }
}

void normalize_vector_sum_d(Vector_d** v_out, const Vector_d* v_in)
{
    double sum;

    calc_vector_asum_d(&sum, v_in);

    if (sum > 0)
    {
        multiply_vector_by_scalar_d(v_out, v_in, 1.0 / sum);
    } 
    else if (v_in != *v_out)
    {
        copy_vector_d(v_out, v_in);
    }
}

void normalize_vector_mag_f(Vector_f** v_out, const Vector_f* v_in)
{
    float mag;

    calc_vector_mag_f(&mag, v_in);

    if (mag > 0)
    {
        multiply_vector_by_scalar_f(v_out, v_in, 1.0f / mag);
    }
    else if (v_in != *v_out)
    {
        copy_vector_f(v_out, v_in);
    }
}

void normalize_vector_mag_d(Vector_d** v_out, const Vector_d* v_in)
{
    double mag;

    calc_vector_mag_d(&mag, v_in);

    if (mag > 0)
    {
        multiply_vector_by_scalar_d(v_out, v_in, 1.0 / mag);
    }
    else if (v_in != *v_out)
    {
        copy_vector_d(v_out, v_in);
    }
}

Error* add_vectors_f
(
    Vector_f**      v_out, 
    const Vector_f* v_1, 
    const Vector_f* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_f(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Adding vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_saxpy(num_elts, 1, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        blas_saxpy(num_elts, 1, v_1->elts, 1, v_2->elts, 1);
    }
    else
    {
        copy_vector_f(v_out, v_2);
        blas_saxpy(num_elts, 1, v_1->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ] += v_2->elts[ elt ];
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ] += v_1->elts[ elt ];
        }
    }
    else
    {
        create_vector_f(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ] = v_1->elts[ elt ] + v_2->elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* add_vectors_d
(
    Vector_d**      v_out, 
    const Vector_d* v_1, 
    const Vector_d* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_d(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Adding vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_daxpy(num_elts, 1, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        blas_daxpy(num_elts, 1, v_1->elts, 1, v_2->elts, 1);
    }
    else
    {
        copy_vector_d(v_out, v_2);
        blas_daxpy(num_elts, 1, v_1->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ] += v_2->elts[ elt ];
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ] += v_1->elts[ elt ];
        }
    }
    else
    {
        create_vector_d(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ] = v_1->elts[ elt ] + v_2->elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* add_vectors_cf
(
    Vector_cf**      v_out, 
    const Vector_cf* v_1, 
    const Vector_cf* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_f z = {1, 0};
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_cf(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Adding vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_caxpy(num_elts, z, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        blas_caxpy(num_elts, z, v_1->elts, 1, v_2->elts, 1);
    }
    else
    {
        copy_vector_cf(v_out, v_2);
        blas_caxpy(num_elts, z, v_1->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            v_1->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            v_2->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
    else
    {
        create_vector_cf(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            (*v_out)->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* add_vectors_cd
(
    Vector_cd**      v_out, 
    const Vector_cd* v_1, 
    const Vector_cd* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Complex_d z = {1, 0};
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_cd(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Adding vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_zaxpy(num_elts, z, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        blas_zaxpy(num_elts, z, v_1->elts, 1, v_2->elts, 1);
    }
    else
    {
        copy_vector_cd(v_out, v_2);
        blas_zaxpy(num_elts, z, v_1->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            v_1->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            v_2->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
    else
    {
        create_vector_cd(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ].r = v_1->elts[ elt ].r + v_2->elts[ elt ].r;
            (*v_out)->elts[ elt ].i = v_1->elts[ elt ].i + v_2->elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* subtract_vectors_f
(
    Vector_f**      v_out, 
    const Vector_f* v_1, 
    const Vector_f* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Vector_f* v = NULL;
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_f(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Subtracting vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_saxpy(num_elts, -1, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        copy_vector_f(&v, v_1);
        blas_saxpy(num_elts, -1, v_2->elts, 1, v->elts, 1);
        copy_vector_f(v_out, v);
        free_vector_f(v);
    }
    else
    {
        copy_vector_f(v_out, v_1);
        blas_saxpy(num_elts, -1, v_2->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
    else
    {
        create_vector_f(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* subtract_vectors_d
(
    Vector_d**      v_out, 
    const Vector_d* v_1, 
    const Vector_d* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Vector_d* v = NULL;
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_d(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Subtracting vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_daxpy(num_elts, -1, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        copy_vector_d(&v, v_1);
        blas_daxpy(num_elts, -1, v_2->elts, 1, v->elts, 1);
        copy_vector_d(v_out, v);
        free_vector_d(v);
    }
    else
    {
        copy_vector_d(v_out, v_1);
        blas_daxpy(num_elts, -1, v_2->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
    else
    {
        create_vector_d(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ] = v_1->elts[ elt ] - v_2->elts[ elt ];
        }
    }
#endif

    return NULL;
}

Error* subtract_vectors_cf
(
    Vector_cf**      v_out, 
    const Vector_cf* v_1, 
    const Vector_cf* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Vector_cf* v = NULL;
    Complex_f  z = {-1, 0};
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_cf(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Subtracting vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_caxpy(num_elts, z, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        copy_vector_cf(&v, v_1);
        blas_caxpy(num_elts, z, v_2->elts, 1, v->elts, 1);
        copy_vector_cf(v_out, v);
        free_vector_cf(v);
    }
    else
    {
        copy_vector_cf(v_out, v_1);
        blas_caxpy(num_elts, z, v_2->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            v_1->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            v_2->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
    else
    {
        create_vector_cf(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            (*v_out)->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
#endif

    return NULL;
}

Error* subtract_vectors_cd
(
    Vector_cd**      v_out, 
    const Vector_cd* v_1, 
    const Vector_cd* v_2
)
{
    uint32_t num_elts;
#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    Vector_cd* v = NULL;
    Complex_d  z = {-1, 0};
#else
    uint32_t elt;
#endif

    if (v_1->num_elts != v_2->num_elts)
    {
        if (*v_out != v_1 && *v_out != v_2)
        {
            free_vector_cd(*v_out); *v_out = NULL;
        }
        return JWSC_EARG("Subtracting vectors requires equal dimensions");
    }

    num_elts = v_1->num_elts;

#if defined JWSC_HAVE_BLAS || defined JWSC_HAVE_ACCELERATE
    if (*v_out == v_1)
    {
        blas_zaxpy(num_elts, z, v_2->elts, 1, v_1->elts, 1);
    }
    else if (*v_out == v_2)
    {
        copy_vector_cd(&v, v_1);
        blas_zaxpy(num_elts, z, v_2->elts, 1, v->elts, 1);
        copy_vector_cd(v_out, v);
        free_vector_cd(v);
    }
    else
    {
        copy_vector_cd(v_out, v_1);
        blas_zaxpy(num_elts, z, v_2->elts, 1, (*v_out)->elts, 1);
    }
#else
    if (*v_out == v_1)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_1->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            v_1->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
    else if (*v_out == v_2)
    {
        for (elt = 0; elt < num_elts; elt++)
        {
            v_2->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            v_2->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
    else
    {
        create_vector_cd(v_out, num_elts);
        for (elt = 0; elt < num_elts; elt++)
        {
            (*v_out)->elts[ elt ].r = v_1->elts[ elt ].r - v_2->elts[ elt ].r;
            (*v_out)->elts[ elt ].i = v_1->elts[ elt ].i - v_2->elts[ elt ].i;
        }
    }
#endif

    return NULL;
}