Loading...
Searching...
No Matches
volk_32f_stddev_and_mean_32f_x2.h
Go to the documentation of this file.
/* -*- c++ -*- */
/*
* Copyright 2012, 2014, 2021 Free Software Foundation, Inc.
*
* This file is part of VOLK
*
* SPDX-License-Identifier: LGPL-3.0-or-later
*/
#ifndef INCLUDED_volk_32f_stddev_and_mean_32f_x2_a_H
#define INCLUDED_volk_32f_stddev_and_mean_32f_x2_a_H
#include <inttypes.h>
#include <math.h>
#include <volk/volk_common.h>
// Youngs and Cramer's Algorithm for calculating std and mean
// Using the methods discussed here:
// https://doi.org/10.1145/3221269.3223036
#ifdef LV_HAVE_GENERIC
static inline void volk_32f_stddev_and_mean_32f_x2_generic(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
const float* in_ptr = inputBuffer;
if (num_points == 0) {
return;
} else if (num_points == 1) {
*stddev = 0.f;
*mean = (*in_ptr);
return;
}
float Sum[2];
float SquareSum[2] = { 0.f, 0.f };
Sum[0] = (*in_ptr++);
Sum[1] = (*in_ptr++);
uint32_t half_points = num_points / 2;
for (uint32_t number = 1; number < half_points; number++) {
float Val0 = (*in_ptr++);
float Val1 = (*in_ptr++);
float n = (float)number;
float n_plus_one = n + 1.f;
float r = 1.f / (n * n_plus_one);
Sum[0] += Val0;
Sum[1] += Val1;
SquareSum[0] += r * powf(n_plus_one * Val0 - Sum[0], 2);
SquareSum[1] += r * powf(n_plus_one * Val1 - Sum[1], 2);
}
SquareSum[0] += SquareSum[1] + .5f / half_points * pow(Sum[0] - Sum[1], 2);
Sum[0] += Sum[1];
uint32_t points_done = half_points * 2;
for (; points_done < num_points; points_done++) {
float Val = (*in_ptr++);
float n = (float)points_done;
float n_plus_one = n + 1.f;
float r = 1.f / (n * n_plus_one);
Sum[0] += Val;
SquareSum[0] += r * powf(n_plus_one * Val - Sum[0], 2);
}
*stddev = sqrtf(SquareSum[0] / num_points);
*mean = Sum[0] / num_points;
}
#endif /* LV_HAVE_GENERIC */
static inline float update_square_sum_1_val(const float SquareSum,
const float Sum,
const uint32_t len,
const float val)
{
// Updates a sum of squares calculated over len values with the value val
float n = (float)len;
float n_plus_one = n + 1.f;
return SquareSum +
1.f / (n * n_plus_one) * (n_plus_one * val - Sum) * (n_plus_one * val - Sum);
}
static inline float add_square_sums(const float SquareSum0,
const float Sum0,
const float SquareSum1,
const float Sum1,
const uint32_t len)
{
// Add two sums of squares calculated over the same number of values, len
float n = (float)len;
return SquareSum0 + SquareSum1 + .5f / n * (Sum0 - Sum1) * (Sum0 - Sum1);
}
static inline void accrue_result(float* PartialSquareSums,
float* PartialSums,
const uint32_t NumberOfPartitions,
const uint32_t PartitionLen)
{
// Add all partial sums and square sums into the first element of the arrays
uint32_t accumulators = NumberOfPartitions;
uint32_t stages = 0;
uint32_t offset = 1;
uint32_t partition_len = PartitionLen;
while (accumulators >>= 1) {
stages++;
} // Integer log2
accumulators = NumberOfPartitions;
for (uint32_t s = 0; s < stages; s++) {
accumulators /= 2;
uint32_t idx = 0;
for (uint32_t a = 0; a < accumulators; a++) {
PartialSquareSums[idx] = add_square_sums(PartialSquareSums[idx],
PartialSums[idx],
PartialSquareSums[idx + offset],
PartialSums[idx + offset],
partition_len);
PartialSums[idx] += PartialSums[idx + offset];
idx += 2 * offset;
}
offset *= 2;
partition_len *= 2;
}
}
#ifdef LV_HAVE_NEON
#include <arm_neon.h>
#include <volk/volk_neon_intrinsics.h>
static inline void volk_32f_stddev_and_mean_32f_x2_neon(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
if (num_points < 8) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
const float* in_ptr = inputBuffer;
__VOLK_ATTR_ALIGNED(32) float SumLocal[8] = { 0.f };
__VOLK_ATTR_ALIGNED(32) float SquareSumLocal[8] = { 0.f };
const uint32_t eigth_points = num_points / 8;
float32x4_t Sum0, Sum1;
Sum0 = vld1q_f32((const float32_t*)in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
Sum1 = vld1q_f32((const float32_t*)in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
float32x4_t SquareSum0 = { 0.f };
float32x4_t SquareSum1 = { 0.f };
float32x4_t Values0, Values1;
float32x4_t Aux0, Aux1;
float32x4_t Reciprocal;
for (uint32_t number = 1; number < eigth_points; number++) {
Values0 = vld1q_f32(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
Values1 = vld1q_f32(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
float n = (float)number;
float n_plus_one = n + 1.f;
Reciprocal = vdupq_n_f32(1.f / (n * n_plus_one));
Sum0 = vaddq_f32(Sum0, Values0);
Aux0 = vdupq_n_f32(n_plus_one);
SquareSum0 =
_neon_accumulate_square_sum_f32(SquareSum0, Sum0, Values0, Reciprocal, Aux0);
Sum1 = vaddq_f32(Sum1, Values1);
Aux1 = vdupq_n_f32(n_plus_one);
SquareSum1 =
_neon_accumulate_square_sum_f32(SquareSum1, Sum1, Values1, Reciprocal, Aux1);
}
vst1q_f32(&SumLocal[0], Sum0);
vst1q_f32(&SumLocal[4], Sum1);
vst1q_f32(&SquareSumLocal[0], SquareSum0);
vst1q_f32(&SquareSumLocal[4], SquareSum1);
accrue_result(SquareSumLocal, SumLocal, 8, eigth_points);
uint32_t points_done = eigth_points * 8;
for (; points_done < num_points; points_done++) {
float val = (*in_ptr++);
SumLocal[0] += val;
SquareSumLocal[0] =
update_square_sum_1_val(SquareSumLocal[0], SumLocal[0], points_done, val);
}
*stddev = sqrtf(SquareSumLocal[0] / num_points);
*mean = SumLocal[0] / num_points;
}
#endif /* LV_HAVE_NEON */
#ifdef LV_HAVE_SSE
#include <volk/volk_sse_intrinsics.h>
#include <xmmintrin.h>
static inline void volk_32f_stddev_and_mean_32f_x2_u_sse(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
if (num_points < 8) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
const float* in_ptr = inputBuffer;
__VOLK_ATTR_ALIGNED(16) float SumLocal[8] = { 0.f };
__VOLK_ATTR_ALIGNED(16) float SquareSumLocal[8] = { 0.f };
const uint32_t eigth_points = num_points / 8;
__m128 Sum0 = _mm_loadu_ps(in_ptr);
in_ptr += 4;
__m128 Sum1 = _mm_loadu_ps(in_ptr);
in_ptr += 4;
__m128 SquareSum0 = _mm_setzero_ps();
__m128 SquareSum1 = _mm_setzero_ps();
__m128 Values0, Values1;
__m128 Aux0, Aux1;
__m128 Reciprocal;
for (uint32_t number = 1; number < eigth_points; number++) {
Values0 = _mm_loadu_ps(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
Values1 = _mm_loadu_ps(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
float n = (float)number;
float n_plus_one = n + 1.f;
Reciprocal = _mm_set_ps1(1.f / (n * n_plus_one));
Sum0 = _mm_add_ps(Sum0, Values0);
Aux0 = _mm_set_ps1(n_plus_one);
SquareSum0 =
_mm_accumulate_square_sum_ps(SquareSum0, Sum0, Values0, Reciprocal, Aux0);
Sum1 = _mm_add_ps(Sum1, Values1);
Aux1 = _mm_set_ps1(n_plus_one);
SquareSum1 =
_mm_accumulate_square_sum_ps(SquareSum1, Sum1, Values1, Reciprocal, Aux1);
}
_mm_store_ps(&SumLocal[0], Sum0);
_mm_store_ps(&SumLocal[4], Sum1);
_mm_store_ps(&SquareSumLocal[0], SquareSum0);
_mm_store_ps(&SquareSumLocal[4], SquareSum1);
accrue_result(SquareSumLocal, SumLocal, 8, eigth_points);
uint32_t points_done = eigth_points * 8;
for (; points_done < num_points; points_done++) {
float val = (*in_ptr++);
SumLocal[0] += val;
SquareSumLocal[0] =
update_square_sum_1_val(SquareSumLocal[0], SumLocal[0], points_done, val);
}
*stddev = sqrtf(SquareSumLocal[0] / num_points);
*mean = SumLocal[0] / num_points;
}
#endif /* LV_HAVE_SSE */
#ifdef LV_HAVE_AVX
#include <immintrin.h>
#include <volk/volk_avx_intrinsics.h>
static inline void volk_32f_stddev_and_mean_32f_x2_u_avx(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
if (num_points < 16) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
const float* in_ptr = inputBuffer;
__VOLK_ATTR_ALIGNED(32) float SumLocal[16] = { 0.f };
__VOLK_ATTR_ALIGNED(32) float SquareSumLocal[16] = { 0.f };
const unsigned int sixteenth_points = num_points / 16;
__m256 Sum0 = _mm256_loadu_ps(in_ptr);
in_ptr += 8;
__m256 Sum1 = _mm256_loadu_ps(in_ptr);
in_ptr += 8;
__m256 SquareSum0 = _mm256_setzero_ps();
__m256 SquareSum1 = _mm256_setzero_ps();
__m256 Values0, Values1;
__m256 Aux0, Aux1;
__m256 Reciprocal;
for (uint32_t number = 1; number < sixteenth_points; number++) {
Values0 = _mm256_loadu_ps(in_ptr);
in_ptr += 8;
__VOLK_PREFETCH(in_ptr + 8);
Values1 = _mm256_loadu_ps(in_ptr);
in_ptr += 8;
__VOLK_PREFETCH(in_ptr + 8);
float n = (float)number;
float n_plus_one = n + 1.f;
Reciprocal = _mm256_set1_ps(1.f / (n * n_plus_one));
Sum0 = _mm256_add_ps(Sum0, Values0);
Aux0 = _mm256_set1_ps(n_plus_one);
SquareSum0 =
_mm256_accumulate_square_sum_ps(SquareSum0, Sum0, Values0, Reciprocal, Aux0);
Sum1 = _mm256_add_ps(Sum1, Values1);
Aux1 = _mm256_set1_ps(n_plus_one);
SquareSum1 =
_mm256_accumulate_square_sum_ps(SquareSum1, Sum1, Values1, Reciprocal, Aux1);
}
_mm256_store_ps(&SumLocal[0], Sum0);
_mm256_store_ps(&SumLocal[8], Sum1);
_mm256_store_ps(&SquareSumLocal[0], SquareSum0);
_mm256_store_ps(&SquareSumLocal[8], SquareSum1);
accrue_result(SquareSumLocal, SumLocal, 16, sixteenth_points);
uint32_t points_done = sixteenth_points * 16;
for (; points_done < num_points; points_done++) {
float val = (*in_ptr++);
SumLocal[0] += val;
SquareSumLocal[0] =
update_square_sum_1_val(SquareSumLocal[0], SumLocal[0], points_done, val);
}
*stddev = sqrtf(SquareSumLocal[0] / num_points);
*mean = SumLocal[0] / num_points;
}
#endif /* LV_HAVE_AVX */
#ifdef LV_HAVE_SSE
#include <xmmintrin.h>
static inline void volk_32f_stddev_and_mean_32f_x2_a_sse(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
if (num_points < 8) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
const float* in_ptr = inputBuffer;
__VOLK_ATTR_ALIGNED(16) float SumLocal[8] = { 0.f };
__VOLK_ATTR_ALIGNED(16) float SquareSumLocal[8] = { 0.f };
const uint32_t eigth_points = num_points / 8;
__m128 Sum0 = _mm_load_ps(in_ptr);
in_ptr += 4;
__m128 Sum1 = _mm_load_ps(in_ptr);
in_ptr += 4;
__m128 SquareSum0 = _mm_setzero_ps();
__m128 SquareSum1 = _mm_setzero_ps();
__m128 Values0, Values1;
__m128 Aux0, Aux1;
__m128 Reciprocal;
for (uint32_t number = 1; number < eigth_points; number++) {
Values0 = _mm_load_ps(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
Values1 = _mm_load_ps(in_ptr);
in_ptr += 4;
__VOLK_PREFETCH(in_ptr + 4);
float n = (float)number;
float n_plus_one = n + 1.f;
Reciprocal = _mm_set_ps1(1.f / (n * n_plus_one));
Sum0 = _mm_add_ps(Sum0, Values0);
Aux0 = _mm_set_ps1(n_plus_one);
SquareSum0 =
_mm_accumulate_square_sum_ps(SquareSum0, Sum0, Values0, Reciprocal, Aux0);
Sum1 = _mm_add_ps(Sum1, Values1);
Aux1 = _mm_set_ps1(n_plus_one);
SquareSum1 =
_mm_accumulate_square_sum_ps(SquareSum1, Sum1, Values1, Reciprocal, Aux1);
}
_mm_store_ps(&SumLocal[0], Sum0);
_mm_store_ps(&SumLocal[4], Sum1);
_mm_store_ps(&SquareSumLocal[0], SquareSum0);
_mm_store_ps(&SquareSumLocal[4], SquareSum1);
accrue_result(SquareSumLocal, SumLocal, 8, eigth_points);
uint32_t points_done = eigth_points * 8;
for (; points_done < num_points; points_done++) {
float val = (*in_ptr++);
SumLocal[0] += val;
SquareSumLocal[0] =
update_square_sum_1_val(SquareSumLocal[0], SumLocal[0], points_done, val);
}
*stddev = sqrtf(SquareSumLocal[0] / num_points);
*mean = SumLocal[0] / num_points;
}
#endif /* LV_HAVE_SSE */
#ifdef LV_HAVE_AVX
#include <immintrin.h>
static inline void volk_32f_stddev_and_mean_32f_x2_a_avx(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
if (num_points < 16) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
const float* in_ptr = inputBuffer;
__VOLK_ATTR_ALIGNED(32) float SumLocal[16] = { 0.f };
__VOLK_ATTR_ALIGNED(32) float SquareSumLocal[16] = { 0.f };
const unsigned int sixteenth_points = num_points / 16;
__m256 Sum0 = _mm256_load_ps(in_ptr);
in_ptr += 8;
__m256 Sum1 = _mm256_load_ps(in_ptr);
in_ptr += 8;
__m256 SquareSum0 = _mm256_setzero_ps();
__m256 SquareSum1 = _mm256_setzero_ps();
__m256 Values0, Values1;
__m256 Aux0, Aux1;
__m256 Reciprocal;
for (uint32_t number = 1; number < sixteenth_points; number++) {
Values0 = _mm256_load_ps(in_ptr);
in_ptr += 8;
__VOLK_PREFETCH(in_ptr + 8);
Values1 = _mm256_load_ps(in_ptr);
in_ptr += 8;
__VOLK_PREFETCH(in_ptr + 8);
float n = (float)number;
float n_plus_one = n + 1.f;
Reciprocal = _mm256_set1_ps(1.f / (n * n_plus_one));
Sum0 = _mm256_add_ps(Sum0, Values0);
Aux0 = _mm256_set1_ps(n_plus_one);
SquareSum0 =
_mm256_accumulate_square_sum_ps(SquareSum0, Sum0, Values0, Reciprocal, Aux0);
Sum1 = _mm256_add_ps(Sum1, Values1);
Aux1 = _mm256_set1_ps(n_plus_one);
SquareSum1 =
_mm256_accumulate_square_sum_ps(SquareSum1, Sum1, Values1, Reciprocal, Aux1);
}
_mm256_store_ps(&SumLocal[0], Sum0);
_mm256_store_ps(&SumLocal[8], Sum1);
_mm256_store_ps(&SquareSumLocal[0], SquareSum0);
_mm256_store_ps(&SquareSumLocal[8], SquareSum1);
accrue_result(SquareSumLocal, SumLocal, 16, sixteenth_points);
uint32_t points_done = sixteenth_points * 16;
for (; points_done < num_points; points_done++) {
float val = (*in_ptr++);
SumLocal[0] += val;
SquareSumLocal[0] =
update_square_sum_1_val(SquareSumLocal[0], SumLocal[0], points_done, val);
}
*stddev = sqrtf(SquareSumLocal[0] / num_points);
*mean = SumLocal[0] / num_points;
}
#endif /* LV_HAVE_AVX */
#ifdef LV_HAVE_RVV
#include <riscv_vector.h>
static inline void volk_32f_stddev_and_mean_32f_x2_rvv(float* stddev,
float* mean,
const float* inputBuffer,
unsigned int num_points)
{
size_t vlmax = __riscv_vsetvlmax_e32m4();
if (num_points < vlmax) {
volk_32f_stddev_and_mean_32f_x2_generic(stddev, mean, inputBuffer, num_points);
return;
}
vfloat32m4_t vsum = __riscv_vle32_v_f32m4(inputBuffer, vlmax);
inputBuffer += vlmax;
vfloat32m4_t vsumsq = __riscv_vfmv_v_f_f32m4(0, vlmax);
size_t partLen = num_points / vlmax;
for (size_t i = 1; i < partLen; ++i, inputBuffer += vlmax) {
vfloat32m4_t v = __riscv_vle32_v_f32m4(inputBuffer, vlmax);
vsum = __riscv_vfadd(vsum, v, vlmax);
vfloat32m4_t vaux = __riscv_vfmsub(v, i + 1.0f, vsum, vlmax);
vaux = __riscv_vfmul(vaux, vaux, vlmax);
vaux = __riscv_vfmul(vaux, 1.0f / (i * (i + 1.0f)), vlmax);
vsumsq = __riscv_vfadd(vsumsq, vaux, vlmax);
}
size_t vl = __riscv_vsetvlmax_e32m2();
vfloat32m2_t vsum2 =
__riscv_vfadd(__riscv_vget_f32m2(vsum, 0), __riscv_vget_f32m2(vsum, 1), vl);
vfloat32m2_t vfix2 =
__riscv_vfsub(__riscv_vget_f32m2(vsum, 0), __riscv_vget_f32m2(vsum, 1), vl);
vfix2 = __riscv_vfmul(vfix2, vfix2, vl);
vfloat32m2_t vsumsq2 =
__riscv_vfadd(__riscv_vget_f32m2(vsumsq, 0), __riscv_vget_f32m2(vsumsq, 1), vl);
vsumsq2 = __riscv_vfmacc(vsumsq2, 0.5f / (num_points / vlmax), vfix2, vl);
vl = __riscv_vsetvlmax_e32m1();
vfloat32m1_t vsum1 =
__riscv_vfadd(__riscv_vget_f32m1(vsum2, 0), __riscv_vget_f32m1(vsum2, 1), vl);
vfloat32m1_t vfix1 =
__riscv_vfsub(__riscv_vget_f32m1(vsum2, 0), __riscv_vget_f32m1(vsum2, 1), vl);
vfix1 = __riscv_vfmul(vfix1, vfix1, vl);
vfloat32m1_t vsumsq1 =
__riscv_vfadd(__riscv_vget_f32m1(vsumsq2, 0), __riscv_vget_f32m1(vsumsq2, 1), vl);
vsumsq1 = __riscv_vfmacc(vsumsq1, 0.5f / (num_points / vlmax * 2), vfix1, vl);
for (size_t n = num_points / vlmax * 4, vl = vlmax >> 2; vl >>= 1; n *= 2) {
vfloat32m1_t vsumdown = __riscv_vslidedown(vsum1, vl, vl);
vfix1 = __riscv_vfsub(vsum1, vsumdown, vl);
vfix1 = __riscv_vfmul(vfix1, vfix1, vl);
vsum1 = __riscv_vfadd(vsum1, vsumdown, vl);
vsumsq1 = __riscv_vfadd(vsumsq1, __riscv_vslidedown(vsumsq1, vl, vl), vl);
vsumsq1 = __riscv_vfmacc(vsumsq1, 0.5f / n, vfix1, vl);
}
float sum = __riscv_vfmv_f(vsum1);
float sumsq = __riscv_vfmv_f(vsumsq1);
for (size_t i = partLen * vlmax; i < num_points; ++i) {
float in = *inputBuffer++;
sum += in;
sumsq = update_square_sum_1_val(sumsq, sum, i, in);
}
*stddev = sqrtf(sumsq / num_points);
*mean = sum / num_points;
}
#endif /*LV_HAVE_RVV*/
#endif /* INCLUDED_volk_32f_stddev_and_mean_32f_x2_a_H */