Векторное ускоренное евклидово расстояние в 3D

Мне нужно выполнить очень обычную и простую матричную операцию.
Однако мне это нужно быстро, очень быстро ...

Я уже рассматриваю многопоточную реализацию, но пока я просто хочу посмотреть, насколько быстро я могу получить ее на одном процессоре.

Матричная операция выглядит следующим образом:

Я вычисляю евклидово расстояние между вектором точек (A) и опорной точкой (B).
Точки находятся в трехмерном пространстве, и каждая точка имеет набор координат X, Y и Z.
Следовательно, Вектор точек описывается тремя массивами с плавающей запятой, содержащими координаты X, Y, Z для каждой точки.
Результатом является другой вектор длины N, содержащий расстояния между каждой точкой в ​​массиве и опорной точкой.

Три массива XYZ расположены как столбцы матрицы Nx3.

x[0]      y[0]      z[0]
x[1]      y[1]      z[1]
x[2]      y[2]      z[2]
x[3]      y[3]      z[3]
.         .         .
.         .         .
.         .         .
x[N-1]    y[N-1]    z[N-1]

В памяти матрица расположена в строчном порядке как одномерный массив, содержащий последовательно значения столбцов X, Y и Z.

x[0], x[1], x[2], x[3] . . . x[N-1], y[0], y[1], y[2], y[3] . . . y[N-1], z[0], z[1], z[2], z[3] . . . z[N-1]

Все это немного усложняется тем фактом, что нам нужно добавить скаляр к каждому члену матрицы, прежде чем мы извлечем квадратный корень.

Ниже приводится процедура в простом коде C:

void calculateDistances3D(float *matrix, float Bx, float By, float Bz, float scalar, float *distances, int N)
{
    float *Ax = matrix;
    float *Ay = Ax + N;
    float *Az = Ay + N;
    int i;

    for (i = 0; i < N; i++) {

        float dx = Ax[i] - Bx;
        float dy = Ay[i] - By;
        float dz = Az[i] - Bz;

        float dx2 = dx * dx;
        float dy2 = dy * dy;
        float dz2 = dz * dz;

        float squaredDistance = dx2 + dy2 + dz2;
        float squaredDistancePlusScalar = squaredDistance + scalar;

        distances[i] = sqrt(squaredDistancePlusScalar);
    }
}

… И вот наивная реализация Accelerate (с использованием vDSP и VecLib):
(обратите внимание, что вся обработка выполняется на месте)

void calculateDistances3D_vDSP(float *matrix, float Bx, float By, float Bz, float scalar, float *distances, int N)
{
    float *Ax = matrix;
    float *Ay = Ax + N;
    float *Az = Ay + N;

    // for each point in the array take the difference with the reference point
    Bx = -Bx;
    By = -By;
    Bz = -Bz;
    vDSP_vsadd(Ax, 1, &Bx, Ax, 1, N);
    vDSP_vsadd(Ay, 1, &By, Ay, 1, N);
    vDSP_vsadd(Az, 1, &Bz, Az, 1, N);

    // square each coordinate
    vDSP_vsq(Ax, 1, Ax, 1, N);
    vDSP_vsq(Ay, 1, Ay, 1, N);
    vDSP_vsq(Az, 1, Az, 1, N);

    // reduce XYZ columns to a single column in Ax (reduction by summation)
    vDSP_vadd(Ax, 1, Ay, 1, Ax, 1, N);
    vDSP_vadd(Ax, 1, Az, 1, Ax, 1, N);

    // add scalar
    vDSP_vsadd(Ax, 1, &scalar, Ax, 1, N);

    // take sqrt
    vvsqrtf(distances, Ax, &N);
}

В библиотеке vDSP для вычисления расстояний между векторами можно использовать следующие функции:

vDSP_vdist()
vDSP_distancesq()
vDSP_vpythg()

Возможно, мне что-то не хватает, но насколько я могу судить, ни один из них не поддерживает три входных вектора, необходимых для расчета расстояния в 3D.

Несколько замечаний:
(1) Я не сравниваю расстояния, поэтому я не могу жить с квадратным расстоянием. Мне нужно реальное расстояние, поэтому вычисление квадратного корня абсолютно необходимо.
(2) Получение обратного квадратного корня может быть возможным, если вы действительно думаете, что это сделает код значительно быстрее.

У меня сложилось впечатление, что я не использую фреймворк Accelerate в полной мере.
Я ищу что-то более умное и, возможно, краткое, выполняющее больше работы с меньшим количеством вызовов функций. Также можно было бы переупорядочить память другими способами, однако я думаю, что структура памяти и так довольно хороша.

Я также открыт для предложений о других высокооптимизированных / векторизованных библиотеках линейной алгебры, которые работают на процессорах Intel. Меня не волнует, являются ли они коммерческими решениями или решениями с открытым исходным кодом, если они работают быстро и надежно.

Возникает вопрос: какая функция или комбинация функций в фреймворке Accelerate лучше всего для достижения более быстрого кода, чем указано выше?

Я разрабатываю в Xcode 7 на MacBook Pro (Retina, 15 дюймов, середина 2014 г.) под управлением Mac OS X El Capitan.

Спасибо.