Commit 39e6a3e4 authored by Vincent Hourdin's avatar Vincent Hourdin
Browse files

fix threading of histogram_median (for MAD and median rejection)

parent 5504ca49
......@@ -571,13 +571,18 @@ void sortnet (WORD *a, size_t n) {
* @param a array of unsigned short to search
* @param n size of the array
* @return median as a double (for n odd)
* Use temp storage h to build the histogram. Complexity O(2*N)
* Use temp storage h to build the histogram. Complexity O(2*N) in time and
* O(N) in memory if single threaded or O(N+N*nb_threads) if multi-threaded
*/
double histogram_median(WORD *a, size_t n, threading_type threading) {
// For arrays n < 10 histogram is use fast and simple sortnet_median
if (n < 10)
return sortnet_median(a, n);
threading = limit_threading(&threading, 2000000, n);
#ifndef _OPENMP
threading = SINGLE_THREADED;
#endif
const size_t s = sizeof(unsigned int);
unsigned int *h = (unsigned int*) calloc(USHRT_MAX + 1, s);
if (!h) {
......@@ -585,39 +590,43 @@ double histogram_median(WORD *a, size_t n, threading_type threading) {
return -1.0;
}
#ifdef _OPENMP
threading = limit_threading(&threading, 400000, n);
#pragma omp parallel num_threads(threading) if (threading > 1)
#endif
{
unsigned int *hthr = (unsigned int*) calloc(USHRT_MAX + 1, s);
if (!hthr) {
PRINT_ALLOC_ERR;
if (threading == SINGLE_THREADED) {
for (size_t i = 0; i < n; i++) {
h[a[i]]++;
}
else {
}
#ifdef _OPENMP
#pragma omp for nowait
#endif
for (size_t i = 0; i < n; i++) {
hthr[a[i]]++;
else {
/* parallel version */
#pragma omp parallel num_threads(threading) if (threading > 1)
{
unsigned int *hthr = (unsigned int*) calloc(USHRT_MAX + 1, s);
if (!hthr) {
PRINT_ALLOC_ERR;
}
#ifdef _OPENMP
else {
#pragma omp for schedule(static) nowait
for (size_t i = 0; i < n; i++) {
hthr[a[i]]++;
}
/* this critical block doesn't only apply to this parallel group,
* but to the callers as well, so it prevents this function from
* running in parallel with one thread for each call */
#pragma omp critical
#endif
{
// add per thread histogram to main histogram
#ifdef _OPENMP
{
// add per-thread histogram to main histogram
#pragma omp simd
#endif
for (int ii = 0; ii <= USHRT_MAX; ++ii) {
h[ii] += hthr[ii];
for (int ii = 0; ii <= USHRT_MAX; ++ii) {
h[ii] += hthr[ii];
}
}
free(hthr);
}
free(hthr);
}
}
unsigned int i= 0, j = 0, k = n / 2;
#endif
unsigned int i = 0, j = 0, k = n / 2;
unsigned int sum = 0;
if (n % 2 == 0) {
for (; sum <= k - 1; j++)
......
......@@ -30,7 +30,6 @@
* It compares the results with the quicksort*/
#define NBTRIES 200 // for result checking, unit test of implementations
#define THREADING_TYPE MULTI_THREADED
cominfo com; // the core data struct
guiinfo gui; // the gui data struct
......@@ -44,7 +43,7 @@ double median_from_sorted_array(WORD *arr, int size)
return (double)sum/2.0;
}
int compare_median_algos(int datasize)
int compare_median_algos(int datasize, int threads)
{
WORD *data, *data_backup;
double result_qsel1, result_qsel2, result_qsort;
......@@ -65,7 +64,7 @@ int compare_median_algos(int datasize)
result_qsel1 = quickmedian(data, datasize);
memcpy(data_backup, data, datasize * sizeof(WORD));
result_qsel2 = histogram_median(data, datasize, THREADING_TYPE);
result_qsel2 = histogram_median(data, datasize, threads);
memcpy(data_backup, data, datasize * sizeof(WORD));
if (result_qsel1 != result_qsort || result_qsel2 != result_qsort) {
......@@ -91,6 +90,10 @@ Test(Sorting, Median)
{
int size = 1;
for (int i = 0; i < NBTRIES; i++, size++) {
cr_assert(compare_median_algos(size) == 0, "Failed at size=%u", size);
cr_assert(compare_median_algos(size, 1) == 0, "Failed at size=%u", size);
}
for (int i = 0; i < NBTRIES; i++, size++) {
cr_assert(compare_median_algos(size, 2) == 0, "Failed at size=%u", size);
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment