From 058db8a3917555260fac2d8505eb071af84973a8 Mon Sep 17 00:00:00 2001 From: willyborn Date: Sat, 2 Aug 2025 23:24:45 +0200 Subject: [PATCH 1/2] sum on CPU now has the same precision as other platforms --- src/backend/cpu/kernel/reduce.hpp | 136 ++++++++++++++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/src/backend/cpu/kernel/reduce.hpp b/src/backend/cpu/kernel/reduce.hpp index de685b426a..e06d57d18e 100644 --- a/src/backend/cpu/kernel/reduce.hpp +++ b/src/backend/cpu/kernel/reduce.hpp @@ -17,6 +17,27 @@ namespace arrayfire { namespace cpu { namespace kernel { +// Since only 1 pass is used, the rounding errors will become +// important because the longer the serie the larger the difference +// between the 2 numbers to sum See: +// https://en.wikipedia.org/wiki/Kahan_summation_algorithm to reduce +// this error +template +To stable_add(To &correction, To out_val, Ti in_val) { + // correction is zero for the first time around + To y = in_val - correction; + // Alas out_val is big, y small, so low-order digits of y are lost + To t = out_val + y; + // (t - out_val) cancels the high order part of y + // subtracting y recovers negative (low part of y) + correction = (t - out_val) - y; + // Algebraically, correction should always be zero. Beware of + // overly-agressive optimizing compilers! + return t; + // Next time around, the lost low part will be added to y in a fresh + // attempt +} + template struct reduce_dim { void operator()(Param out, const dim_t outOffset, CParam in, @@ -62,6 +83,32 @@ struct reduce_dim { } }; +template +struct reduce_dim { + common::Transform, compute_t, af_add_t> transform; + common::Binary, af_add_t> reduce; + void operator()(Param out, const dim_t outOffset, CParam in, + const dim_t inOffset, const int dim, bool change_nan, + double nanval) { + const af::dim4 istrides = in.strides(); + const af::dim4 idims = in.dims(); + + data_t *const outPtr = out.get() + outOffset; + data_t const *const inPtr = in.get() + inOffset; + dim_t stride = istrides[dim]; + + compute_t out_val = common::Binary, af_add_t>::init(); + compute_t correction = scalar>(0); + for (dim_t i = 0; i < idims[dim]; i++) { + compute_t in_val = transform(inPtr[i * stride]); + if (change_nan) in_val = IS_NAN(in_val) ? nanval : in_val; + out_val = stable_add(correction, out_val, in_val); + } + + *outPtr = data_t(out_val); + } +}; + template void n_reduced_keys(Param okeys, int *n_reduced, CParam keys) { const af::dim4 kdims = keys.dims(); @@ -159,6 +206,55 @@ struct reduce_dim_by_key { } }; +template +struct reduce_dim_by_key { + common::Transform, compute_t, af_add_t> transform; + void operator()(Param ovals, const dim_t ovOffset, CParam keys, + CParam vals, const dim_t vOffset, int *n_reduced, + const int dim, bool change_nan, double nanval) { + const af::dim4 vstrides = vals.strides(); + const af::dim4 vdims = vals.dims(); + + const af::dim4 ovstrides = ovals.strides(); + + data_t const *const inKeysPtr = keys.get(); + data_t const *const inValsPtr = vals.get(); + data_t *const outValsPtr = ovals.get(); + + int keyidx = 0; + compute_t current_key = compute_t(inKeysPtr[0]); + compute_t out_val = common::Binary, af_add_t>::init(); + compute_t correction = scalar>(0); + + dim_t istride = vstrides[dim]; + dim_t ostride = ovstrides[dim]; + + for (dim_t i = 0; i < vdims[dim]; i++) { + compute_t keyval = inKeysPtr[i]; + + if (keyval == current_key) { + compute_t in_val = + transform(inValsPtr[vOffset + (i * istride)]); + if (change_nan) in_val = IS_NAN(in_val) ? nanval : in_val; + out_val = stable_add(correction, out_val, in_val); + + } else { + outValsPtr[ovOffset + (keyidx * ostride)] = out_val; + + current_key = keyval; + correction = scalar>(0); + out_val = transform(inValsPtr[vOffset + (i * istride)]); + if (change_nan) out_val = IS_NAN(out_val) ? nanval : out_val; + ++keyidx; + } + + if (i == (vdims[dim] - 1)) { + outValsPtr[ovOffset + (keyidx * ostride)] = out_val; + } + } + } +}; + template struct reduce_all { common::Transform, compute_t, op> transform; @@ -199,6 +295,46 @@ struct reduce_all { } }; +template +struct reduce_all { + common::Transform, compute_t, af_add_t> transform; + void operator()(Param out, CParam in, bool change_nan, + double nanval) { + // Decrement dimension of select dimension + af::dim4 dims = in.dims(); + af::dim4 strides = in.strides(); + const data_t *inPtr = in.get(); + data_t *const outPtr = out.get(); + + compute_t out_val = common::Binary, af_add_t>::init(); + compute_t correction = scalar>(0); + + for (dim_t l = 0; l < dims[3]; l++) { + dim_t off3 = l * strides[3]; + + for (dim_t k = 0; k < dims[2]; k++) { + dim_t off2 = k * strides[2]; + + for (dim_t j = 0; j < dims[1]; j++) { + dim_t off1 = j * strides[1]; + + for (dim_t i = 0; i < dims[0]; i++) { + dim_t idx = i + off1 + off2 + off3; + + compute_t in_val = transform(inPtr[idx]); + if (change_nan) { + in_val = IS_NAN(in_val) ? nanval : in_val; + } + out_val = stable_add(correction, out_val, in_val); + } + } + } + } + + *outPtr = data_t(out_val); + } +}; + } // namespace kernel } // namespace cpu } // namespace arrayfire From 4a8dca4d27a27b6a872a3b3221a7be13ec31fc89 Mon Sep 17 00:00:00 2001 From: willyborn Date: Fri, 8 Aug 2025 20:01:55 +0200 Subject: [PATCH 2/2] Fixes to mean, on all platforms --- src/backend/cpu/kernel/mean.hpp | 45 +++- src/backend/cpu/mean.cpp | 43 +-- src/backend/cuda/kernel/mean.hpp | 122 ++++++--- src/backend/oneapi/kernel/mean.hpp | 187 +++++++++---- src/backend/opencl/kernel/mean.hpp | 185 ++++++++++--- src/backend/opencl/kernel/mean_dim.cl | 18 +- src/backend/opencl/kernel/mean_ops.cl | 35 ++- test/mean.cpp | 371 ++++++++++++++++++++++---- 8 files changed, 775 insertions(+), 231 deletions(-) diff --git a/src/backend/cpu/kernel/mean.hpp b/src/backend/cpu/kernel/mean.hpp index c15773687e..0ed4f424a2 100644 --- a/src/backend/cpu/kernel/mean.hpp +++ b/src/backend/cpu/kernel/mean.hpp @@ -38,6 +38,45 @@ struct MeanOp { } }; +template +struct MeanOpWithCorrection { + common::Transform transform; + To runningMean; + Tw runningCount; + To correction; + MeanOpWithCorrection(Ti mean, Tw count) + : transform() + , runningMean(transform(mean)) + , runningCount(count) + , correction(scalar(0)) {} + + void operator()(Ti newMean, Tw newCount) { + runningCount += newCount; + if ((newCount != 0) || (runningCount != 0)) { // + // Since only 1 pass is used, the rounding errors will become + // important because the longer the serie the larger the + // difference between the 2 numbers to sum See: + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm to + // reduce this error + // correction is zero for the first time around + To y = + (transform(newMean) - runningMean) * (newCount / runningCount) - + correction; + // Alas, runningMean is big, y small, so low-order digits of y + // are lost + To t = runningMean + y; + // (t - runningMean) cancels the high order part of y + // subtracting y recovers negative (low part of y) + correction = (t - runningMean) - y; + // Algebraically, correction should always be zero. Beware + // overly-agressive optimizing compilers! + runningMean = t; + // Next time around, the lost low part will be added to y in a + // fresh attempt + } + } +}; + template struct mean_weighted_dim { void operator()(Param output, const dim_t outOffset, @@ -74,7 +113,8 @@ struct mean_weighted_dim { dim_t istride = istrides[dim]; dim_t wstride = wstrides[dim]; - MeanOp, compute_t, compute_t> Op(0, 0); + MeanOpWithCorrection, compute_t, compute_t> Op(0, + 0); for (dim_t i = 0; i < idims[dim]; i++) { Op(compute_t(in[inOffset + i * istride]), compute_t(wt[wtOffset + i * wstride])); @@ -113,7 +153,8 @@ struct mean_dim { dim_t istride = istrides[dim]; dim_t end = inOffset + idims[dim] * istride; - MeanOp, compute_t, compute_t> Op(0, 0); + MeanOpWithCorrection, compute_t, compute_t> Op(0, + 0); for (dim_t i = inOffset; i < end; i += istride) { Op(compute_t(in[i]), 1); } diff --git a/src/backend/cpu/mean.cpp b/src/backend/cpu/mean.cpp index 2323442110..8bb59a6caa 100644 --- a/src/backend/cpu/mean.cpp +++ b/src/backend/cpu/mean.cpp @@ -63,32 +63,35 @@ Array mean(const Array &in, const Array &wt, const int dim) { template T mean(const Array &in, const Array &wt) { - using MeanOpT = kernel::MeanOp, compute_t, compute_t>; + using MeanOpT = + kernel::MeanOpWithCorrection, compute_t, compute_t>; in.eval(); wt.eval(); getQueue().sync(); - af::dim4 dims = in.dims(); - af::dim4 strides = in.strides(); - const T *inPtr = in.get(); - const Tw *wtPtr = wt.get(); - - auto input = compute_t(inPtr[0]); - auto weight = compute_t(wtPtr[0]); - MeanOpT Op(input, weight); + const af::dim4 &dims = in.dims(); + const af::dim4 &istrides = in.strides(); + const af::dim4 &wstrides = wt.strides(); + const T *inPtr = in.get(); + const Tw *wtPtr = wt.get(); + // Split workload in equal parts, to improve accuracy or larger arrays + MeanOpT Op(0, 0); for (dim_t l = 0; l < dims[3]; l++) { - dim_t off3 = l * strides[3]; + dim_t ioff3 = l * istrides[3]; + dim_t woff3 = l * wstrides[3]; for (dim_t k = 0; k < dims[2]; k++) { - dim_t off2 = k * strides[2]; + dim_t ioff2 = k * istrides[2]; + dim_t woff2 = k * wstrides[2]; for (dim_t j = 0; j < dims[1]; j++) { - dim_t off1 = j * strides[1]; + dim_t ioff1 = j * istrides[1]; + dim_t woff1 = j * wstrides[1]; for (dim_t i = 0; i < dims[0]; i++) { - dim_t idx = i + off1 + off2 + off3; - Op(compute_t(inPtr[idx]), compute_t(wtPtr[idx])); + Op(compute_t(inPtr[i + ioff1 + ioff2 + ioff3]), + compute_t(wtPtr[i + woff1 + woff2 + woff3])); } } } @@ -99,16 +102,16 @@ T mean(const Array &in, const Array &wt) { template To mean(const Array &in) { - using MeanOpT = kernel::MeanOp, compute_t, compute_t>; + using MeanOpT = kernel::MeanOpWithCorrection, compute_t, + compute_t>; in.eval(); getQueue().sync(); - af::dim4 dims = in.dims(); - af::dim4 strides = in.strides(); - const Ti *inPtr = in.get(); + const af::dim4 &dims = in.dims(); + const af::dim4 &strides = in.strides(); + const Ti *inPtr = in.get(); MeanOpT Op(0, 0); - for (dim_t l = 0; l < dims[3]; l++) { dim_t off3 = l * strides[3]; @@ -120,7 +123,7 @@ To mean(const Array &in) { for (dim_t i = 0; i < dims[0]; i++) { dim_t idx = i + off1 + off2 + off3; - Op(compute_t(inPtr[idx]), 1); + Op(compute_t(inPtr[idx]), 1.); } } } diff --git a/src/backend/cuda/kernel/mean.hpp b/src/backend/cuda/kernel/mean.hpp index a26eeac7fd..7c7818adf5 100644 --- a/src/backend/cuda/kernel/mean.hpp +++ b/src/backend/cuda/kernel/mean.hpp @@ -39,13 +39,40 @@ namespace kernel { template __device__ __host__ void stable_mean(To *lhs, Tw *l_wt, To rhs, Tw r_wt) { - if (((*l_wt) != (Tw)0) || (r_wt != (Tw)0)) { - Tw l_scale = (*l_wt); - (*l_wt) += r_wt; + Tw l_scale = (*l_wt); + (*l_wt) += (r_wt); + if (((*l_wt) != (Tw)0) || ((r_wt) != (Tw)0)) { l_scale = l_scale / (*l_wt); - Tw r_scale = r_wt / (*l_wt); - (*lhs) = (l_scale * *lhs) + (r_scale * rhs); + Tw r_scale = (r_wt) / (*l_wt); + (*lhs) = (l_scale * (*lhs)) + (r_scale * (rhs)); + } +} + +template +__device__ __host__ void stable_mean(To *c, To *lhs, Tw *l_wt, To rhs, + Tw r_wt) { + (*l_wt) += (r_wt); + if (((*l_wt) != (Tw)0) || ((r_wt) != (Tw)0)) { + // Since only 1 pass is used, the rounding errors will become important + // because the longer the serie the larger the difference between the 2 + // numbers to sum See: + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm to reduce + // this error + // (*lhs) = (*lhs) + (rhs - (*lhs)) * r_wt / (*l_wt); + + // *c is zero for the first time around + To y = ((r_wt) / (*l_wt)) * ((rhs) - (*lhs)) - (*c); + // Alas, (*lhs) is big, y small, so low-order digits of y are lost + To t = (*lhs) + y; + // (t - (*lhs)) cancels the high-order part of y + // subtracting y recovers negative (low part of y) + (*c) = (t - (*lhs)) - y; + // Algebraically, *c should always be zero. Beware overly-agressive + // optimizing compilers! + (*lhs) = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt } } @@ -74,24 +101,32 @@ __global__ static void mean_dim_kernel(Param out, Param owt, int ooffset = ids[3] * out.strides[3] + ids[2] * out.strides[2] + ids[1] * out.strides[1] + ids[0]; + optr += ooffset; + if (owptr != NULL) { + int owoffset = ids[3] * owt.strides[3] + ids[2] * owt.strides[2] + + ids[1] * owt.strides[1] + ids[0]; + owptr += owoffset; + } + // There is only one element per block for out // There are blockDim.y elements per block for in // Hence increment ids[dim] just after offseting out and before offsetting // in - optr += ooffset; - if (owptr != NULL) owptr += ooffset; - const uint blockIdx_dim = ids[dim]; - - ids[dim] = ids[dim] * blockDim.y + tidy; + ids[dim] = ids[dim] * blockDim.y + tidy; int ioffset = ids[3] * in.strides[3] + ids[2] * in.strides[2] + ids[1] * in.strides[1] + ids[0]; iptr += ioffset; - if (iwptr != NULL) iwptr += ioffset; + if (iwptr != NULL) { + int iwoffset = ids[3] * iwt.strides[3] + ids[2] * iwt.strides[2] + + ids[1] * iwt.strides[1] + ids[0]; + iwptr += iwoffset; + } - const uint id_dim_in = ids[dim]; - const uint istride_dim = in.strides[dim]; + const uint id_dim_in = ids[dim]; + const uint istride_dim = in.strides[dim]; + const uint iwstride_dim = iwt.strides[dim]; bool is_valid = (ids[0] < in.dims[0]) && (ids[1] < in.dims[1]) && (ids[2] < in.dims[2]) && (ids[3] < in.dims[3]); @@ -106,7 +141,7 @@ __global__ static void mean_dim_kernel(Param out, Param owt, if (iwptr != NULL) { weight = *iwptr; } else { - weight = (Tw)1; + weight = compute_t(1); } } @@ -115,16 +150,17 @@ __global__ static void mean_dim_kernel(Param out, Param owt, __shared__ compute_t s_val[THREADS_X * DIMY]; __shared__ compute_t s_idx[THREADS_X * DIMY]; + compute_t correction = scalar>(0); for (int id = id_dim_in_start; is_valid && (id < in.dims[dim]); id += offset_dim * blockDim.y) { iptr = iptr + offset_dim * blockDim.y * istride_dim; if (iwptr != NULL) { - iwptr = iwptr + offset_dim * blockDim.y * istride_dim; - stable_mean(&val, &weight, transform(*iptr), compute_t(*iwptr)); + iwptr = iwptr + offset_dim * blockDim.y * iwstride_dim; + stable_mean(&correction, &val, &weight, transform(*iptr), + compute_t(*iwptr)); } else { - // Faster version of stable_mean when iwptr is NULL - val = val + (transform(*iptr) - val) / (weight + (Tw)1); - weight = weight + (Tw)1; + stable_mean(&correction, &val, &weight, transform(*iptr), + compute_t(1)); } } @@ -299,16 +335,16 @@ __global__ static void mean_first_kernel(Param out, Param owt, __shared__ compute_t s_val[THREADS_PER_BLOCK]; __shared__ compute_t s_idx[THREADS_PER_BLOCK]; + compute_t correction = scalar>(0); if (iwptr != NULL) { for (int id = xid + DIMX; id < lim; id += DIMX) { - stable_mean(&val, &weight, transform(iptr[id]), + stable_mean(&correction, &val, &weight, transform(iptr[id]), compute_t(iwptr[id])); } } else { for (int id = xid + DIMX; id < lim; id += DIMX) { - // Faster version of stable_mean when iwptr is NULL - val = val + (transform(iptr[id]) - val) / (weight + (Tw)1); - weight = weight + (Tw)1; + stable_mean(&correction, &val, &weight, transform(iptr[id]), + compute_t(1)); } } @@ -406,8 +442,7 @@ void mean_first(Param out, CParam in, CParam iwt) { threads_x); if (blocks_x > 1) { - Param owt; - owt.ptr = NULL; + Array owt = createEmptyArray(dim4()); mean_first_launcher(out, owt, tmpOut, tmpWt, 1, blocks_y, threads_x); } @@ -425,7 +460,7 @@ void mean_weighted(Param out, CParam in, CParam iwt, int dim) { template void mean(Param out, CParam in, int dim) { - Param dummy_weight; + Array dummy_weight = createEmptyArray(dim4()); mean_weighted(out, in, dummy_weight, dim); } @@ -433,17 +468,16 @@ template T mean_all_weighted(CParam in, CParam iwt) { int in_elements = in.dims[0] * in.dims[1] * in.dims[2] * in.dims[3]; - // FIXME: Use better heuristics to get to the optimum number - if (in_elements > 4096) { - bool in_is_linear = (in.strides[0] == 1); - bool wt_is_linear = (iwt.strides[0] == 1); - for (int k = 1; k < 4; k++) { - in_is_linear &= - (in.strides[k] == (in.strides[k - 1] * in.dims[k - 1])); - wt_is_linear &= - (iwt.strides[k] == (iwt.strides[k - 1] * iwt.dims[k - 1])); - } + bool in_is_linear = (in.strides[0] == 1); + bool wt_is_linear = (iwt.strides[0] == 1); + for (int k = 1; k < 4; k++) { + in_is_linear &= (in.strides[k] == (in.strides[k - 1] * in.dims[k - 1])); + wt_is_linear &= + (iwt.strides[k] == (iwt.strides[k - 1] * iwt.dims[k - 1])); + } + // FIXME: Use better heuristics to get to the optimum number + if (in_elements > 4096 || !in_is_linear || !wt_is_linear) { if (in_is_linear && wt_is_linear) { in.dims[0] = in_elements; for (int k = 1; k < 4; k++) { @@ -487,9 +521,11 @@ T mean_all_weighted(CParam in, CParam iwt) { compute_t val = static_cast>(h_ptr[0]); compute_t weight = static_cast>(h_wptr[0]); + compute_t correction = + common::Binary, af_add_t>::init(); for (int i = 1; i < tmp_elements; i++) { - stable_mean(&val, &weight, compute_t(h_ptr[i]), + stable_mean(&correction, &val, &weight, compute_t(h_ptr[i]), compute_t(h_wptr[i])); } @@ -508,8 +544,10 @@ T mean_all_weighted(CParam in, CParam iwt) { compute_t val = static_cast>(h_ptr[0]); compute_t weight = static_cast>(h_wptr[0]); + compute_t correction = + common::Binary, af_add_t>::init(); for (int i = 1; i < in_elements; i++) { - stable_mean(&val, &weight, compute_t(h_ptr[i]), + stable_mean(&correction, &val, &weight, compute_t(h_ptr[i]), compute_t(h_wptr[i])); } @@ -566,9 +604,11 @@ To mean_all(CParam in) { compute_t val = static_cast>(h_ptr[0]); compute_t weight = static_cast>(h_cptr[0]); + compute_t correction = + common::Binary, af_add_t>::init(); for (int i = 1; i < tmp_elements; i++) { - stable_mean(&val, &weight, compute_t(h_ptr[i]), + stable_mean(&correction, &val, &weight, compute_t(h_ptr[i]), compute_t(h_cptr[i])); } @@ -583,11 +623,13 @@ To mean_all(CParam in) { common::Transform, af_add_t> transform; compute_t count = static_cast>(1); + compute_t correction = + common::Binary, af_add_t>::init(); compute_t val = transform(h_ptr[0]); compute_t weight = count; for (int i = 1; i < in_elements; i++) { - stable_mean(&val, &weight, transform(h_ptr[i]), count); + stable_mean(&correction, &val, &weight, transform(h_ptr[i]), count); } return static_cast(val); diff --git a/src/backend/oneapi/kernel/mean.hpp b/src/backend/oneapi/kernel/mean.hpp index 4c8533b1ec..84e613ff1f 100644 --- a/src/backend/oneapi/kernel/mean.hpp +++ b/src/backend/oneapi/kernel/mean.hpp @@ -32,15 +32,70 @@ namespace oneapi { namespace kernel { +// Calculates the weighted mean of 2 items with similar weight +// INPUT +// *lhs: value (init 0. or first value) +// *l_wt: weight (init 0. or first weight) +// rhs : extra value +// r_wt: extra weight +// OUTPUT +// *lhs: cumulative mean +// *l_wt: cumulative weight +// template void stable_mean(To *lhs, Tw *l_wt, To rhs, Tw r_wt) { - if (((*l_wt) != (Tw)0) || (r_wt != (Tw)0)) { - Tw l_scale = (*l_wt); - (*l_wt) += r_wt; + Tw l_scale = (*l_wt); + (*l_wt) += (r_wt); + if (((*l_wt) != (Tw)0) || ((r_wt) != (Tw)0)) { l_scale = l_scale / (*l_wt); - Tw r_scale = r_wt / (*l_wt); - (*lhs) = (l_scale * *lhs) + (r_scale * rhs); + Tw r_scale = (r_wt) / (*l_wt); + (*lhs) = (l_scale * (*lhs)) + (r_scale * (rhs)); + + /* Alternative, although bigger error for values of same size + if (r_wt != (Tw)0) { + (*l_wt) += r_wt; + Tw r_scale = r_wt / (*l_wt); + (*lhs) = (*lhs) + (rhs - (*lhs)) * r_scale; + } + */ + } +} + +// Calculates the weighted mean of a series of items in 1 pass +// INPUT +// *c: estimated correction from previous item (init=0.) +// *lhs: value (init=0.0 or first value) +// *l_wt: weight (init=0.0 or first weight) +// rhs : extra value +// r_wt: extra weight +// OUTPUT +// *c: estimated correction for next item +// *lhs: cumulative mean +// *l_wt: cumulative weight +// +// Since only 1 pass is used, the rounding errors will become important because +// the longer the serie the larger the difference between the 2 numbers to sum +// See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm to reduce this +// error +template +void stable_mean(To *c, To *lhs, Tw *l_wt, To rhs, Tw r_wt) { + (*l_wt) += (r_wt); + if (((*l_wt) != (Tw)0) || ((r_wt) != (Tw)0)) { + // (*lhs) = (*lhs) + (rhs - (*lhs)) * r_wt / (*l_wt); + + // *c is zero for the first time around + To y = ((rhs) - (*lhs)) * ((r_wt) / (*l_wt)) - (*c); + // Alas, (*lhs) is big, y small, so low-order digits of y are lost + To t = (*lhs) + y; + // (t - (*lhs)) cancels the high-order part of y + // subtracting y recovers negative (low part of y) + (*c) = (t - (*lhs)) - y; + // Algebraically, *c should always be zero. Beware overly-agressive + // optimizing compilers! + (*lhs) = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt } } @@ -92,11 +147,19 @@ class meanDimKernelSMEM { uint ooffset = ids[3] * oInfo_.strides[3] + ids[2] * oInfo_.strides[2] + ids[1] * oInfo_.strides[1] + ids[0] + oInfo_.offset; + optr += ooffset; + + Tw *owptr = nullptr; + if (output_weight_) { + uint owoffset = + ids[3] * owInfo_.strides[3] + ids[2] * owInfo_.strides[2] + + ids[1] * owInfo_.strides[1] + ids[0] + owInfo_.offset; + owptr = owt_.get_pointer() + owoffset; + } // There is only one element per block for out // There are blockDim.y elements per block for in // Hence increment ids[dim] just after offseting out and before // offsetting in - optr += ooffset; const uint blockIdx_dim = ids[dim]; ids[dim] = ids[dim] * g.get_local_range(1) + lidy; @@ -106,13 +169,16 @@ class meanDimKernelSMEM { iptr += ioffset; const Tw *iwptr = nullptr; - Tw *owptr = nullptr; - - if (output_weight_) owptr = owt_.get_pointer() + ooffset; - if (input_weight_) iwptr = iwt_.get_pointer() + ioffset; + if (input_weight_) { + uint iwoffset = + ids[3] * iwInfo_.strides[3] + ids[2] * iwInfo_.strides[2] + + ids[1] * iwInfo_.strides[1] + ids[0] + iwInfo_.offset; + iwptr = iwt_.get_pointer() + iwoffset; + } - const uint id_dim_in = ids[dim]; - const uint istride_dim = iInfo_.strides[dim]; + const uint id_dim_in = ids[dim]; + const uint istride_dim = iInfo_.strides[dim]; + const uint iwstride_dim = iwInfo_.strides[dim]; bool is_valid = (ids[0] < iInfo_.dims[0]) && (ids[1] < iInfo_.dims[1]) && @@ -125,7 +191,7 @@ class meanDimKernelSMEM { if (is_valid && id_dim_in < iInfo_.dims[dim]) { val = transform(*iptr); - if (iwptr) { + if (input_weight_) { weight = *iwptr; } else { weight = (Tw)1; @@ -134,19 +200,19 @@ class meanDimKernelSMEM { const uint id_dim_in_start = id_dim_in + offset_dim_ * g.get_local_range(1); + compute_t correction = scalar>(0); for (int id = id_dim_in_start; is_valid && (id < iInfo_.dims[dim]); id += offset_dim_ * g.get_local_range(1)) { iptr = iptr + offset_dim_ * g.get_local_range(1) * istride_dim; if (input_weight_) { iwptr = - iwptr + offset_dim_ * g.get_local_range(1) * istride_dim; - stable_mean(&val, &weight, transform(*iptr), + iwptr + offset_dim_ * g.get_local_range(1) * iwstride_dim; + stable_mean(&correction, &val, &weight, transform(*iptr), compute_t(*iwptr)); } else { - // Faster version of stable_mean when iwptr is NULL - val = val + (transform(*iptr) - val) / (weight + (Tw)1); - weight = weight + (Tw)1; + stable_mean(&correction, &val, &weight, transform(*iptr), + compute_t(1)); } } @@ -540,8 +606,7 @@ void mean_first(Param out, Param in, Param iwt) { threads_x); if (blocks_x > 1) { - Param owt; - owt.data = nullptr; + Array owt = createEmptyArray(dim4()); mean_first_launcher(out, owt, tmpOut, tmpWt, 1, blocks_y, threads_x); } @@ -559,7 +624,7 @@ void mean_weighted(Param out, Param in, Param iwt, int dim) { template void mean(Param out, Param in, int dim) { - Param dummy_weight; + Array dummy_weight = createEmptyArray(dim4()); mean_weighted(out, in, dummy_weight, dim); } @@ -567,17 +632,17 @@ template T mean_all_weighted(Param in, Param iwt) { uintl in_elements = in.info.dims[0] * in.info.dims[1] * in.info.dims[2] * in.info.dims[3]; - // FIXME: Use better heuristics to get to the optimum number - if (in_elements > 4096) { - bool in_is_linear = (in.info.strides[0] == 1); - bool wt_is_linear = (iwt.info.strides[0] == 1); - for (int k = 1; k < 4; k++) { - in_is_linear &= (in.info.strides[k] == - (in.info.strides[k - 1] * in.info.dims[k - 1])); - wt_is_linear &= (iwt.info.strides[k] == - (iwt.info.strides[k - 1] * iwt.info.dims[k - 1])); - } + bool in_is_linear = in.info.strides[0] == 1; + bool wt_is_linear = iwt.info.strides[0] == 1; + for (int k = 1; k < 4; k++) { + in_is_linear &= (in.info.strides[k] == + (in.info.strides[k - 1] * in.info.dims[k - 1])); + wt_is_linear &= (iwt.info.strides[k] == + (iwt.info.strides[k - 1] * iwt.info.dims[k - 1])); + } + // FIXME: Use better heuristics to get to the optimum number + if (in_elements > 4096 || !in_is_linear || !wt_is_linear) { if (in_is_linear && wt_is_linear) { in.info.dims[0] = in_elements; for (int k = 1; k < 4; k++) { @@ -610,21 +675,21 @@ T mean_all_weighted(Param in, Param iwt) { compute_t val; auto tmpOut_get = tmpOut.get(); - auto tmpWt_get = tmpWt.get(); + auto tmpWt_get = tmpWt.get(); getQueue() .submit([&](sycl::handler &h) { - auto acc_in = - tmpOut_get->get_host_access(h, sycl::read_only); - auto acc_wt = - tmpWt_get->get_host_access(h, sycl::read_only); + auto acc_in = tmpOut_get->get_host_access(h, sycl::read_only); + auto acc_wt = tmpWt_get->get_host_access(h, sycl::read_only); h.host_task([acc_in, acc_wt, tmp_elements, &val] { val = static_cast>(acc_in[0]); compute_t weight = static_cast>(acc_wt[0]); + compute_t correction = static_cast>(0); for (int i = 1; i < tmp_elements; i++) { - stable_mean(&val, &weight, compute_t(acc_in[i]), + stable_mean(&correction, &val, &weight, + compute_t(acc_in[i]), compute_t(acc_wt[i])); } }); @@ -639,13 +704,18 @@ T mean_all_weighted(Param in, Param iwt) { h, sycl::range{in_elements}, sycl::read_only); auto acc_wt = iwt.data->get_host_access( h, sycl::range{in_elements}, sycl::read_only); - - h.host_task([acc_in, acc_wt, in_elements, &val]() { - val = acc_in[0]; - compute_t weight = acc_wt[0]; + dim_t offset_in = in.info.offset; + dim_t offset_wt = iwt.info.offset; + + h.host_task([acc_in, offset_in, acc_wt, offset_wt, in_elements, + &val]() { + val = acc_in[offset_in]; + compute_t weight = acc_wt[offset_wt]; + compute_t correction = static_cast>(0); for (int i = 1; i < in_elements; i++) { - stable_mean(&val, &weight, compute_t(acc_in[i]), - compute_t(acc_wt[i])); + stable_mean(&correction, &val, &weight, + compute_t(acc_in[i + offset_in]), + compute_t(acc_wt[i + offset_wt])); } }); }) @@ -688,7 +758,7 @@ To mean_all(Param in) { Array tmpOut = createEmptyArray(outDims); Array tmpCt = createEmptyArray(outDims); - Param iwt; + Array iwt = createEmptyArray(dim4()); mean_first_launcher(tmpOut, tmpCt, in, iwt, blocks_x, blocks_y, threads_x); @@ -696,20 +766,20 @@ To mean_all(Param in) { compute_t val; auto tmpOut_get = tmpOut.get(); - auto tmpCt_get = tmpCt.get(); + auto tmpCt_get = tmpCt.get(); getQueue() .submit([&](sycl::handler &h) { - auto out = - tmpOut_get->get_host_access(h, sycl::read_only); - auto ct = - tmpCt_get->get_host_access(h, sycl::read_only); + auto out = tmpOut_get->get_host_access(h, sycl::read_only); + auto ct = tmpCt_get->get_host_access(h, sycl::read_only); h.host_task([out, ct, tmp_elements, &val] { val = static_cast>(out[0]); compute_t weight = static_cast>(ct[0]); + compute_t correction = scalar>(0); for (int i = 1; i < tmp_elements; i++) { - stable_mean(&val, &weight, compute_t(out[i]), + stable_mean(&correction, &val, &weight, + compute_t(out[i]), compute_t(ct[i])); } }); @@ -720,16 +790,19 @@ To mean_all(Param in) { compute_t val; getQueue() .submit([&](sycl::handler &h) { - auto acc_in = - in.data->get_host_access(h, sycl::read_only); - h.host_task([acc_in, in_elements, &val]() { + auto acc_in = in.data->get_host_access(h, sycl::read_only); + dim_t offset_in = in.info.offset; + h.host_task([acc_in, offset_in, in_elements, &val]() { common::Transform, af_add_t> transform; compute_t count = static_cast>(1); - val = transform(acc_in[0]); - compute_t weight = count; - for (int i = 1; i < in_elements; i++) { - stable_mean(&val, &weight, transform(acc_in[i]), count); + val = transform(acc_in[offset_in]); + compute_t weight = count; + compute_t correction = static_cast>(0); + for (int i = 1 + offset_in; i < in_elements + offset_in; + i++) { + stable_mean(&correction, &val, &weight, + transform(acc_in[i]), count); } }); }) diff --git a/src/backend/opencl/kernel/mean.hpp b/src/backend/opencl/kernel/mean.hpp index bc80a23be9..06cf6cd781 100644 --- a/src/backend/opencl/kernel/mean.hpp +++ b/src/backend/opencl/kernel/mean.hpp @@ -38,10 +38,10 @@ struct MeanOp { MeanOp(T mean, Tw count) : runningMean(mean), runningCount(count) {} void operator()(T newMean, Tw newCount) { + Tw runningScale = runningCount; + Tw newScale = newCount; + runningCount += newCount; if ((newCount != 0) || (runningCount != 0)) { - Tw runningScale = runningCount; - Tw newScale = newCount; - runningCount += newCount; runningScale = runningScale / runningCount; newScale = newScale / (Tw)runningCount; runningMean = (runningScale * runningMean) + (newScale * newMean); @@ -56,10 +56,10 @@ struct MeanOp { MeanOp(cfloat mean, float count) : runningMean(mean), runningCount(count) {} void operator()(cfloat newMean, float newCount) { + float runningScale = runningCount; + float newScale = newCount; + runningCount += newCount; if ((newCount != 0) || (runningCount != 0)) { - float runningScale = runningCount; - float newScale = newCount; - runningCount += newCount; runningScale = runningScale / runningCount; newScale = newScale / (float)runningCount; runningMean.s[0] = @@ -78,10 +78,10 @@ struct MeanOp { : runningMean(mean), runningCount(count) {} void operator()(cdouble newMean, double newCount) { + double runningScale = runningCount; + double newScale = newCount; + runningCount += newCount; if ((newCount != 0) || (runningCount != 0)) { - double runningScale = runningCount; - double newScale = newCount; - runningCount += newCount; runningScale = runningScale / runningCount; newScale = newScale / (double)runningCount; runningMean.s[0] = @@ -92,6 +92,106 @@ struct MeanOp { } }; +template +struct MeanOpWithCorrection { + T runningMean; + Tw runningCount; + T correction; + MeanOpWithCorrection(T mean, Tw count) + : runningMean(mean), runningCount(count), correction(scalar(0)) {} + + void operator()(T newMean, Tw newCount) { + runningCount += newCount; + if ((newCount != 0) || (runningCount != 0)) { + // correction is zero for the first time around + T y = + (newMean - runningMean) * newCount / runningCount - correction; + // Alas, runningMean is big, y small, so low-order digits of y are + // lost + T t = runningMean + y; + // (t - runningMean) cancels the high order part of y + // subtracting y recovers negative (low part of y) + correction = (t - runningMean) - y; + // Algebraically, correction should always be zero. Beware + // overly-agressive optimizing compilers! + runningMean = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt + } + } +}; + +template<> +struct MeanOpWithCorrection { + cfloat runningMean; + float runningCount; + cfloat correction; + MeanOpWithCorrection(cfloat mean, float count) + : runningMean(mean) + , runningCount(count) + , correction(scalar(0)) {} + + void operator()(cfloat newMean, float newCount) { + runningCount += newCount; + if ((newCount != 0) || (runningCount != 0)) { + // correction is zero for the first time around + cfloat y{ + (newMean.s[0] - runningMean.s[0]) * newCount / runningCount - + correction.s[0], + (newMean.s[1] - runningMean.s[1]) * newCount / runningCount - + correction.s[1]}; + // Alas, runningMean is big, y small, so low-order digits of y are + // lost + cfloat t{runningMean.s[0] + y.s[0], runningMean.s[1] + y.s[1]}; + // (t - runningMean) cancels the high order part of y + // subtracting y recovers negative (low part of y) + correction.s[0] = (t.s[0] - runningMean.s[0]) - y.s[0]; + correction.s[1] = (t.s[1] - runningMean.s[1]) - y.s[1]; + // Algebraically, correction should always be zero. Beware + // overly-agressive optimizing compilers! + runningMean = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt + } + } +}; + +template<> +struct MeanOpWithCorrection { + cdouble runningMean; + double runningCount; + cdouble correction; + MeanOpWithCorrection(cdouble mean, double count) + : runningMean(mean) + , runningCount(count) + , correction(scalar(0)) {} + + void operator()(cdouble newMean, double newCount) { + runningCount += newCount; + if ((newCount != 0) || (runningCount != 0)) { + // correction is zero for the first time around + cdouble y{ + (newMean.s[0] - runningMean.s[0]) * newCount / runningCount - + correction.s[0], + (newMean.s[1] - runningMean.s[1]) * newCount / runningCount - + correction.s[1]}; + // Alas, runningMean is big, y small, so low-order digits of y are + // lost + cdouble t{t.s[0] = runningMean.s[0] + y.s[0], + t.s[1] = runningMean.s[1] + y.s[1]}; + // (t - runningMean) cancels the high order part of y + // subtracting y recovers negative (low part of y) + correction.s[0] = (t.s[0] - runningMean.s[0]) - y.s[0]; + correction.s[1] = (t.s[1] - runningMean.s[1]) - y.s[1]; + // Algebraically, correction should always be zero. Beware + // overly-agressive optimizing compilers! + runningMean = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt + } + } +}; + template void meanDimLauncher(Param out, Param owt, Param in, Param inWeight, const int dim, const int threads_y, @@ -175,12 +275,12 @@ void meanDim(Param out, Param in, Param inWeight, int dim) { meanDimLauncher(tmpOut, tmpWeight, in, inWeight, dim, threads_y, groups_all); - Param owt; + Array owt = createEmptyArray(dim4()); groups_all[dim] = 1; meanDimLauncher(out, owt, tmpOut, tmpWeight, dim, threads_y, groups_all); } else { - Param tmpWeight; + Array tmpWeight = createEmptyArray(dim4()); meanDimLauncher(out, tmpWeight, in, inWeight, dim, threads_y, groups_all); } @@ -258,18 +358,10 @@ void meanFirst(Param out, Param in, Param inWeight) { uint groups_x = divup(in.info.dims[0], threads_x * REPEAT); uint groups_y = divup(in.info.dims[1], threads_y); - Param tmpOut = out; - Param noWeight; - noWeight.info.offset = 0; - for (int k = 0; k < 4; ++k) { - noWeight.info.dims[k] = 0; - noWeight.info.strides[k] = 0; - } - // Does not matter what the value is it will not be used. Just needs to be - // valid. - noWeight.data = inWeight.data; + Param tmpOut(out); + Param noWeight(nullptr, {0}); - Param tmpWeight = noWeight; + Param tmpWeight(noWeight); if (groups_x > 1) { tmpOut.data = bufferAlloc(groups_x * in.info.dims[1] * in.info.dims[2] * @@ -307,7 +399,7 @@ void meanWeighted(Param out, Param in, Param inWeight, int dim) { template void mean(Param out, Param in, int dim) { - Param noWeight; + Array noWeight = createEmptyArray(dim4()); meanWeighted(out, in, noWeight, dim); } @@ -316,25 +408,28 @@ T meanAllWeighted(Param in, Param inWeight) { int in_elements = in.info.dims[0] * in.info.dims[1] * in.info.dims[2] * in.info.dims[3]; - // FIXME: Use better heuristics to get to the optimum number - if (in_elements > 4096) { - bool in_is_linear = (in.info.strides[0] == 1); - bool wt_is_linear = (in.info.strides[0] == 1); - for (int k = 1; k < 4; k++) { - in_is_linear &= (in.info.strides[k] == - (in.info.strides[k - 1] * in.info.dims[k - 1])); - wt_is_linear &= - (inWeight.info.strides[k] == - (inWeight.info.strides[k - 1] * inWeight.info.dims[k - 1])); - } + bool in_is_linear = (in.info.strides[0] == 1); + bool wt_is_linear = (inWeight.info.strides[0] == 1); + for (int k = 1; k < 4; k++) { + in_is_linear &= (in.info.strides[k] == + (in.info.strides[k - 1] * in.info.dims[k - 1])); + wt_is_linear &= + (inWeight.info.strides[k] == + (inWeight.info.strides[k - 1] * inWeight.info.dims[k - 1])); + } + // FIXME: Use better heuristics to get to the optimum number + if (in_elements > 4096 || !in_is_linear || !wt_is_linear) { if (in_is_linear && wt_is_linear) { in.info.dims[0] = in_elements; for (int k = 1; k < 4; k++) { in.info.dims[k] = 1; in.info.strides[k] = in_elements; } - inWeight.info = in.info; + for (int k = 0; k < 4; ++k) { + inWeight.info.dims[k] = in.info.dims[k]; + inWeight.info.strides[k] = in.info.strides[k]; + } } uint threads_x = nextpow2(std::max(32u, (uint)in.info.dims[0])); @@ -344,8 +439,10 @@ T meanAllWeighted(Param in, Param inWeight) { uint groups_x = divup(in.info.dims[0], threads_x * REPEAT); uint groups_y = divup(in.info.dims[1], threads_y); - Array tmpOut = createEmptyArray(groups_x); - Array tmpWeight = createEmptyArray(groups_x); + dim4 outDims(groups_x, in.info.dims[1], in.info.dims[2], + in.info.dims[3]); + Array tmpOut = createEmptyArray(outDims); + Array tmpWeight = createEmptyArray(outDims); meanFirstLauncher(tmpOut, tmpWeight, in, inWeight, threads_x, groups_x, groups_y); @@ -362,8 +459,8 @@ T meanAllWeighted(Param in, Param inWeight) { compute_t initial = static_cast>(h_ptr[0]); compute_t w = static_cast>(h_wptr[0]); - MeanOp, compute_t> Op(initial, w); - for (int i = 1; i < (int)tmpOut.elements(); i++) { + MeanOpWithCorrection, compute_t> Op(initial, w); + for (int i = 1; i < (int)h_ptr.size(); i++) { Op(compute_t(h_ptr[i]), compute_t(h_wptr[i])); } @@ -381,7 +478,7 @@ T meanAllWeighted(Param in, Param inWeight) { compute_t initial = static_cast>(h_ptr[0]); compute_t w = static_cast>(h_wptr[0]); - MeanOp, compute_t> Op(initial, w); + MeanOpWithCorrection, compute_t> Op(initial, w); for (int i = 1; i < (int)in_elements; i++) { Op(compute_t(h_ptr[i]), compute_t(h_wptr[i])); } @@ -422,7 +519,7 @@ To meanAll(Param in) { Array tmpOut = createEmptyArray(outDims); Array tmpCt = createEmptyArray(outDims); - Param iWt; + Array iWt = createEmptyArray(dim4()); meanFirstLauncher(tmpOut, tmpCt, in, iWt, threads_x, groups_x, groups_y); @@ -438,7 +535,7 @@ To meanAll(Param in) { compute_t initial = static_cast>(h_ptr[0]); compute_t w = static_cast>(h_cptr[0]); - MeanOp, compute_t> Op(initial, w); + MeanOpWithCorrection, compute_t> Op(initial, w); for (int i = 1; i < (int)h_ptr.size(); i++) { Op(compute_t(h_ptr[i]), compute_t(h_cptr[i])); } @@ -454,8 +551,8 @@ To meanAll(Param in) { // TODO : MeanOp with (Tw)1 common::Transform, af_add_t> transform; common::Transform, af_add_t> transform_weight; - MeanOp, compute_t> Op(transform(h_ptr[0]), - transform_weight(1)); + MeanOpWithCorrection, compute_t> Op( + transform(h_ptr[0]), transform_weight(1)); for (int i = 1; i < (int)in_elements; i++) { Op(transform(h_ptr[i]), transform_weight(1)); } diff --git a/src/backend/opencl/kernel/mean_dim.cl b/src/backend/opencl/kernel/mean_dim.cl index 9448486391..bb899cffdc 100644 --- a/src/backend/opencl/kernel/mean_dim.cl +++ b/src/backend/opencl/kernel/mean_dim.cl @@ -37,8 +37,8 @@ kernel void meanDim(global To *oData, KParam oInfo, ids[1] * oInfo.strides[1] + ids[0] + oInfo.offset; #ifdef OUTPUT_WEIGHT - owData += ids[3] * oInfo.strides[3] + ids[2] * oInfo.strides[2] + - ids[1] * oInfo.strides[1] + ids[0] + oInfo.offset; + owData += ids[3] * owInfo.strides[3] + ids[2] * owInfo.strides[2] + + ids[1] * owInfo.strides[1] + ids[0] + owInfo.offset; #endif const uint id_dim_out = ids[kDim]; @@ -48,8 +48,10 @@ kernel void meanDim(global To *oData, KParam oInfo, ids[1] * iInfo.strides[1] + ids[0] + iInfo.offset; #ifdef INPUT_WEIGHT - iwData += ids[3] * iInfo.strides[3] + ids[2] * iInfo.strides[2] + - ids[1] * iInfo.strides[1] + ids[0] + iInfo.offset; + iwData += ids[3] * iwInfo.strides[3] + ids[2] * iwInfo.strides[2] + + ids[1] * iwInfo.strides[1] + ids[0] + iwInfo.offset; + + const uint iwstride_dim = iwInfo.strides[kDim]; #endif const uint id_dim_in = ids[kDim]; @@ -76,17 +78,19 @@ kernel void meanDim(global To *oData, KParam oInfo, const uint id_dim_in_start = id_dim_in + group_dim * get_local_size(1); #ifdef INPUT_WEIGHT + To correction = (To)(0); for (int id = id_dim_in_start; is_valid && (id < iInfo.dims[kDim]); id += group_dim * get_local_size(1)) { iData = iData + group_dim * get_local_size(1) * istride_dim; - iwData = iwData + group_dim * get_local_size(1) * istride_dim; - binOp(&out_val, &out_wt, transform(*iData), *iwData); + iwData = iwData + group_dim * get_local_size(1) * iwstride_dim; + binOpWithCorr(&correction, &out_val, &out_wt, transform(*iData), *iwData); } #else + To correction = (To)(0); for (int id = id_dim_in_start; is_valid && (id < iInfo.dims[kDim]); id += group_dim * get_local_size(1)) { iData = iData + group_dim * get_local_size(1) * istride_dim; - binOp(&out_val, &out_wt, transform(*iData), one_Tw); + binOpWithCorr(&correction, &out_val, &out_wt, transform(*iData), one_Tw); } #endif diff --git a/src/backend/opencl/kernel/mean_ops.cl b/src/backend/opencl/kernel/mean_ops.cl index b4f104f6ba..795ac133a4 100644 --- a/src/backend/opencl/kernel/mean_ops.cl +++ b/src/backend/opencl/kernel/mean_ops.cl @@ -10,12 +10,37 @@ To transform(Ti in) { return (To)(in); } void binOp(To *lhs, Tw *l_wt, To rhs, Tw r_wt) { - if (((*l_wt) != 0) || (r_wt != 0)) { - Tw l_scale = (*l_wt); - (*l_wt) += r_wt; + Tw l_scale = (*l_wt); + (*l_wt) += (r_wt); + if (((*l_wt) != 0) || ((r_wt) != 0)) { l_scale = l_scale / (*l_wt); - Tw r_scale = r_wt / (*l_wt); - (*lhs) = (l_scale * (*lhs)) + (r_scale * rhs); + Tw r_scale = (r_wt) / (*l_wt); + (*lhs) = (l_scale * (*lhs)) + (r_scale * (rhs)); + } +} + +// +// Since only 1 pass is used, the rounding errors will become important because +// the longer the serie the larger the difference between the 2 numbers to sum +// See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm to reduce this +// error +void binOpWithCorr(To *c, To *lhs, Tw *l_wt, To rhs, Tw r_wt) { + (*l_wt) += (r_wt); + if (((*l_wt) != 0) || ((r_wt) != 0)) { + // (*lhs) = (*lhs) + ((rhs) - (*lhs)) * ((r_wt) / (*l_wt)); + + // *c is zero for the first time around + To y = ((rhs) - (*lhs)) * ((r_wt) / (*l_wt)) - (*c); + // Alas, (*lhs) is big, y small, so low-order digits of y are lost + To t = (*lhs) + y; + // (t - (*lhs)) cancels the high-order part of y + // subtracting y recovers negative (low part of y) + (*c) = (t - (*lhs)) - y; + // Algebraically, *c should always be zero. Beware overly-agressive + // optimizing compilers! + (*lhs) = t; + // Next time around, the lost low part will be added to y in a fresh + // attempt } } diff --git a/test/mean.cpp b/test/mean.cpp index 79dd76db2d..5539b7effb 100644 --- a/test/mean.cpp +++ b/test/mean.cpp @@ -82,8 +82,9 @@ void meanDimTest(string pFileName, dim_t dim, bool isWeighted = false) { SUPPORTED_TYPE_CHECK(T); SUPPORTED_TYPE_CHECK(outType); - double tol = 1.0e-3; - if ((af_dtype)af::dtype_traits::af_type == f16) tol = 4.e-3; + // The precision in the test files goes up to 1e-4. The mean functions + // achieve 5e-5 for all types except half. + double tol = ((af_dtype)af::dtype_traits::af_type == f16) ? 4e-3 : 1e-4; vector numDims; vector> in; vector> tests; @@ -102,7 +103,7 @@ void meanDimTest(string pFileName, dim_t dim, bool isWeighted = false) { vector outData(dims.elements()); - outArray.host((void*)outData.data()); + outArray.host((void *)outData.data()); vector currGoldBar(tests[0].begin(), tests[0].end()); @@ -124,7 +125,7 @@ void meanDimTest(string pFileName, dim_t dim, bool isWeighted = false) { vector outData(dims.elements()); - outArray.host((void*)outData.data()); + outArray.host((void *)outData.data()); vector currGoldBar(tests[0].begin(), tests[0].end()); @@ -168,8 +169,32 @@ TYPED_TEST(Mean, Wtd_Dim1Matrix) { true); } +template +T random() { + return T(std::rand() % 10); +} + +template<> +half random() { + // create values from -0.5 to 0.5 to ensure sum does not deviate + // too far out of half's useful range + float r = static_cast(rand()) / static_cast(RAND_MAX) - 0.5f; + return half(r); +} + +template<> +cfloat random() { + return cfloat(float(std::rand() % 10), float(std::rand() % 10)); +} + +template<> +cdouble random() { + return cdouble(double(std::rand() % 10), double(std::rand() % 10)); +} + template void meanAllTest(T const_value, dim4 dims) { + UNUSED(const_value); typedef typename meanOutType::type outType; SUPPORTED_TYPE_CHECK(T); @@ -177,21 +202,48 @@ void meanAllTest(T const_value, dim4 dims) { using af::array; using af::mean; + std::srand(std::time(0)); - vector hundred(dims.elements(), const_value); + vector data(dims.elements()); + std::generate(data.begin(), data.end(), random); - outType gold = outType(0); - // for(auto i:hundred) gold += i; - for (int i = 0; i < (int)hundred.size(); i++) { gold = gold + hundred[i]; } - gold = gold / dims.elements(); + // Process vector by replacing each 2 elements by its weighted mean. + // Repeat until only 1 element remains + // Weighted means is necessary since 1 element is not processed in case of + // odd #elements. During the processing on the next round, the weight of + // this 1 element is different from the others. + vector> meanWeight(data.size()); + std::transform(data.cbegin(), data.cend(), meanWeight.begin(), + [](auto d) { return std::make_pair((outType)d, 1.); }); + + auto mwEnd = meanWeight.end(); + do { + auto saveIt = meanWeight.begin(); + // When odd #elements remain, skip the first. + if (std::distance(saveIt, mwEnd) % 2 == 1) ++saveIt; + for (auto mwIt = saveIt; mwIt != mwEnd; ++mwIt, ++mwIt, ++saveIt) { + auto nextIt = mwIt + 1; + saveIt->second = mwIt->second + nextIt->second; + if (saveIt->second != 0) { + // When the weight is 0 for both elements, we do not care about + // the value + const float nextScale = nextIt->second / saveIt->second; + const float currScale = 1.0 - nextScale; + saveIt->first = + mwIt->first * currScale + nextIt->first * nextScale; + } + } + // From now on, only process the saved elements (#elements/2) + mwEnd = saveIt; + } while (mwEnd != meanWeight.begin() + 1); + outType gold = meanWeight[0].first; - array a(dims, &(hundred.front())); + array a(dims, &(data.front())); outType output = mean(a); - ASSERT_NEAR(::real(output), ::real(gold), 1.0e-3); - ASSERT_NEAR(::imag(output), ::imag(gold), 1.0e-3); + ASSERT_NEAR(::real(output), ::real(gold), 1.e-5); + ASSERT_NEAR(::imag(output), ::imag(gold), 1.e-5); } - template<> void meanAllTest(half_float::half const_value, dim4 dims) { SUPPORTED_TYPE_CHECK(half_float::half); @@ -212,12 +264,11 @@ void meanAllTest(half_float::half const_value, dim4 dims) { // make sure output2 and output are binary equals. This is necessary // because af_half is not a complete type half output2_copy; - memcpy(static_cast(&output2_copy), &output2, sizeof(af_half)); + memcpy(static_cast(&output2_copy), &output2, sizeof(af_half)); ASSERT_EQ(output, output2_copy); - ASSERT_NEAR(output, gold, 1.0e-3); + ASSERT_NEAR(output, gold, 5e-4); } - TEST(MeanAll, f64) { meanAllTest(2.1, dim4(10, 10, 1, 1)); } TEST(MeanAll, f32) { meanAllTest(2.1f, dim4(10, 5, 2, 1)); } @@ -240,29 +291,6 @@ TEST(MeanAll, u16) { meanAllTest(2, dim4(100, 1, 1, 1)); } TEST(MeanAll, c64) { meanAllTest(cdouble(2.1), dim4(10, 10, 1, 1)); } -template -T random() { - return T(std::rand() % 10); -} - -template<> -half random() { - // create values from -0.5 to 0.5 to ensure sum does not deviate - // too far out of half's useful range - float r = static_cast(rand()) / static_cast(RAND_MAX) - 0.5f; - return half(r); -} - -template<> -cfloat random() { - return cfloat(float(std::rand() % 10), float(std::rand() % 10)); -} - -template<> -cdouble random() { - return cdouble(double(std::rand() % 10), double(std::rand() % 10)); -} - template class WeightedMean : public ::testing::Test { public: @@ -290,29 +318,53 @@ void weightedMeanAllTest(dim4 dims) { std::generate(data.begin(), data.end(), random); std::generate(wts.begin(), wts.end(), random); - outType wtdSum = outType(0); - wtsType wtsSum = wtsType(0); - - for (int i = 0; i < (int)data.size(); i++) { - wtdSum = wtdSum + data[i] * wts[i]; - wtsSum = wtsSum + wts[i]; - } - - outType gold = wtdSum / outType(wtsSum); + // Process vector by replacing each 2 elements by its weighted mean. + // Repeat until only 1 element remains + vector> meanWeight(data.size()); + std::transform( + data.cbegin(), data.cend(), wts.cbegin(), meanWeight.begin(), + [](auto d, auto w) { return std::make_pair((outType)d, w); }); + + auto mwEnd = meanWeight.end(); + do { + auto saveIt = meanWeight.begin(); + // When odd #elements remaining, skip the first. + if (std::distance(saveIt, mwEnd) % 2 == 1) ++saveIt; + for (auto mwIt = saveIt; mwIt != mwEnd; ++mwIt, ++mwIt, ++saveIt) { + auto nextIt = mwIt + 1; + saveIt->second = mwIt->second + nextIt->second; + if (saveIt->second != 0) { + // When the weight is 0 for both elements, we do not care about + // the value + const wtsType nextScale = nextIt->second / saveIt->second; + saveIt->first = + mwIt->first + (nextIt->first - mwIt->first) * nextScale; + } + } + // From now on, only process the saved elements (#elements/2) + mwEnd = saveIt; + } while (mwEnd != meanWeight.begin() + 1); + outType gold = meanWeight[0].first; array a(dims, &(data.front())); array w(dims, &(wts.front())); outType output = mean(a, w); - ASSERT_NEAR(::real(output), ::real(gold), 1.0e-2); - ASSERT_NEAR(::imag(output), ::imag(gold), 1.0e-2); + double tol = + ((af_dtype)af::dtype_traits::af_type == f16) ? 5.e-4 : 1.e-5; + ASSERT_NEAR(::real(output), ::real(gold), tol); + ASSERT_NEAR(::imag(output), ::imag(gold), tol); +} + +TYPED_TEST(WeightedMean, Small) { + weightedMeanAllTest(dim4(20, 2, 2, 2)); } TYPED_TEST(WeightedMean, Basic) { weightedMeanAllTest(dim4(32, 30, 33, 17)); } -TEST(WeightedMean, Broadacst) { +TEST(WeightedMean, Broadcast) { float val = 0.5f; array a = randu(4096, 32); array w = constant(val, a.dims()); @@ -344,15 +396,13 @@ TEST(Mean, Issue2093) { float outVal; out.host(&outVal); - float expected = 0.0; + double expected = 0.0; for (size_t i = 0; i < NELEMS; ++i) expected += hdata[i]; expected /= NELEMS; - - ASSERT_NEAR(outVal, expected, 0.001); + EXPECT_NEAR(outVal, expected, 1.0e-5); } -TEST(MeanAll, SubArray) { - // Fixes Issue 2636 +TEST(MeanAll, Issue2636) { using af::mean; using af::span; using af::sum; @@ -379,3 +429,212 @@ TEST(MeanHalf, dim0) { // 0.506836 ASSERT_ARRAYS_NEAR(m16.as(f32), m32, 0.001f); } + +#define TESTS_TEMP_FORMATS_ALL(form) \ + TEST(TEMP_FORMAT, form##_all) { \ + const dim4 dims(2, 2, 2, 2); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const float out = af::mean(toTempFormat(form, in)); \ + const float gold = af::mean(in); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5); \ + } \ + TEST(TEMP_FORMAT, form##_all_vector) { \ + const dim4 dims(20, 1, 1, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const float out = af::mean(toTempFormat(form, in)); \ + const float gold = af::mean(in); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5); \ + } \ + TEST(TEMP_FORMAT, form##_all_weighted) { \ + const dim4 dims(2, 2, 2, 2); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims) * 10.; \ + af::eval(in, weight); \ + \ + const float out = af::mean(toTempFormat(form, in), \ + toTempFormat(form, weight)); \ + const float out2 = af::mean(toTempFormat(form, in), weight); \ + const float out3 = af::mean(in, toTempFormat(form, weight)); \ + const float gold = af::mean(in, weight); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5) << "in & weight TempFormat"; \ + EXPECT_NEAR(out2, gold, 1.0e-5) << "in TempFormat & weight Linear"; \ + EXPECT_NEAR(out3, gold, 1.0e-5) << "in Linear & weight TempFormat"; \ + } \ + TEST(TEMP_FORMAT, form##_all_weighted_vector) { \ + const dim4 dims(20, 1, 1, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims) * 10.; \ + af::eval(in, weight); \ + \ + const float out = af::mean(toTempFormat(form, in), \ + toTempFormat(form, weight)); \ + const float out2 = af::mean(toTempFormat(form, in), weight); \ + const float out3 = af::mean(in, toTempFormat(form, weight)); \ + const float gold = af::mean(in, weight); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5) << "in & weight TempFormat"; \ + EXPECT_NEAR(out2, gold, 1.0e-5) << "in TempFormat & weight Linear"; \ + EXPECT_NEAR(out3, gold, 1.0e-5) << "in Linear & weight TempFormat"; \ + } \ + \ + TEST(TEMP_FORMAT, form##_all_large) { \ + const dim4 dims(2, 512, 60, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const float out = af::mean(toTempFormat(form, in)); \ + const float gold = af::mean(in); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5); \ + } \ + \ + TEST(TEMP_FORMAT, form##_all_weighted_large) { \ + const dim4 dims(2, 512, 60, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims) * 10.; \ + af::eval(in, weight); \ + \ + const float out = af::mean(toTempFormat(form, in), \ + toTempFormat(form, weight)); \ + const float out2 = af::mean(toTempFormat(form, in), weight); \ + const float out3 = af::mean(in, toTempFormat(form, weight)); \ + const float gold = af::mean(in, weight); \ + \ + EXPECT_NEAR(out, gold, 1.0e-5) << "in & weight TempFormat"; \ + EXPECT_NEAR(out2, gold, 1.0e-5) << "in TempFormat & weight Linear"; \ + EXPECT_NEAR(out3, gold, 1.0e-5) << "in Linear & weight TempFormat"; \ + } + +#define TESTS_TEMP_FORMAT_dim(form, dim) \ + TEST(TEMP_FORMAT, form##_##dim) { \ + const dim4 dims(2, 2, 2, 2); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const array out = af::mean(toTempFormat(form, in), dim); \ + const array gold = af::mean(in, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5); \ + } \ + TEST(TEMP_FORMAT, form##_##dim##_vector) { \ + const dim4 dims(20, 1, 1, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const array out = af::mean(toTempFormat(form, in), dim); \ + const array gold = af::mean(in, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5); \ + } \ + TEST(TEMP_FORMAT, form##_##dim##_weighted) { \ + const dim4 dims(2, 2, 2, 2); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims) * 10.; \ + af::eval(in, weight); \ + \ + const array out = \ + af::mean(toTempFormat(form, in), toTempFormat(form, weight), dim); \ + const array out2 = af::mean(toTempFormat(form, in), weight, dim); \ + const array out3 = af::mean(in, toTempFormat(form, weight), dim); \ + const array gold = af::mean(in, weight, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5) << "in & weight TempFormat"; \ + EXPECT_ARRAYS_NEAR(out2, gold, 1.0e-5) \ + << "in TempFormat & weight Linear"; \ + EXPECT_ARRAYS_NEAR(out3, gold, 1.0e-5) \ + << "in Linear & weight TempFormat"; \ + } \ + TEST(TEMP_FORMAT, form##_##dim##_weighted_vector) { \ + const dim4 dims(20, 1, 1, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims) * 10.; \ + af::eval(in, weight); \ + \ + const array out = \ + af::mean(toTempFormat(form, in), toTempFormat(form, weight), dim); \ + const array gold = af::mean(in, weight, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5); \ + } \ + \ + TEST(TEMP_FORMAT, form##_##dim##_large) { \ + const dim4 dims(2, 512, 60, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + in.eval(); \ + \ + const array out = af::mean(toTempFormat(form, in), dim); \ + const array gold = af::mean(in, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5); \ + } \ + \ + TEST(TEMP_FORMAT, form##_##dim##_weighted_large) { \ + const dim4 dims(2, 512, 60, 1); \ + /* Make sure that when a randum numbers are used, that they differ \ + * from the random numbers used in the creation of sub-arrays, so \ + * amplify the valid values by 10x */ \ + const array in = randu(dims) * 10.; \ + const array weight = randu(dims); \ + af::eval(in, weight); \ + \ + const array out = \ + af::mean(toTempFormat(form, in), toTempFormat(form, weight), dim); \ + const array out2 = af::mean(toTempFormat(form, in), weight, dim); \ + const array out3 = af::mean(in, toTempFormat(form, weight), dim); \ + const array gold = af::mean(in, weight, dim); \ + \ + EXPECT_ARRAYS_NEAR(out, gold, 1.0e-5) << "in & weight TempFormat"; \ + EXPECT_ARRAYS_NEAR(out2, gold, 1.0e-5) \ + << "in TempFormat & weight Linear"; \ + EXPECT_ARRAYS_NEAR(out3, gold, 1.0e-5) \ + << "in Linear & weight TempFormat"; \ + } + +#define TESTS_TEMP_FORMATS_dim(form) \ + TESTS_TEMP_FORMAT_dim(form, 0); \ + TESTS_TEMP_FORMAT_dim(form, 1); \ + TESTS_TEMP_FORMAT_dim(form, 2); \ + TESTS_TEMP_FORMAT_dim(form, 3); + +FOREACH_TEMP_FORMAT(TESTS_TEMP_FORMATS_ALL) +FOREACH_TEMP_FORMAT(TESTS_TEMP_FORMATS_dim) \ No newline at end of file