Radix Sort

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <time.h>
#include <omp.h>

#define NO_MKL
//#define TEST
#define TEST2

//Ignore the test functions, I didn't bother removing them

void checkNullPtr(void *ptr){
    if(ptr == NULL){
        printf("Allocation Failed\n");
        exit(1);

    }
}

#ifdef TEST2
void printArr(int *arr, size_t length){
    for(int i=0; i<length; i++){
        printf("%d ", arr[i]);
    }
    printf("\n");
}

void checkAscending(int *arr, size_t length){
    bool ascending = true;
    if(length <= 1){
        printf("%zd element in array\n", length);
    }else{
        for(int i=1; i<length; i++){
            if(arr[i] < arr[i - 1]){
                ascending = false;
            }
        }
        if(ascending){
            printf("Ascending\n");
        }else{
            printf("Not Ascending - Test Failed\n");
        }
    }
}

#endif

//Was gonna write a version using mkl - hasn't happened yet
#ifdef NO_MKL
int getMax(int *arr, size_t length){
    int max = 0;
    for(int i=0; i<length; i++){
        if(arr[i] > max){
            max = arr[i];
        }
    }
    return max;
}
#endif

inline int tgtPlace(int value, int base, int power){
    return  (int)(value / pow(base, power)) % base;
}

void countingSort(int *arr, int *out, size_t length, int base, int power, int *count, int *partitions, int pSize){
    int digit;
    memset(count, 0, sizeof(int) * base * pSize);
    #ifdef TEST
    printf("Counting\n");
    #endif
    #pragma omp parallel
    {
        #pragma omp for private(digit)
        for(int i=0; i<pSize; i++){
            if(i != pSize - 1){
                for(int j=partitions[i]; j<partitions[i + 1]; j++){
                    digit = tgtPlace(arr[j], base, power);
                    count[i + digit * pSize + 1]++;
                }
            }else{
                for(int j=partitions[i]; j<length; j++){
                    digit = tgtPlace(arr[j], base, power);
                    if(digit < base - 1){
                        count[i + digit * pSize + 1]++;
                    }
                }
            }
        }

        #pragma omp barrier
        if(omp_get_thread_num() == 0){
            #ifdef TEST
            printf("Histogram\n");
            #endif
            for(int i = 1; i < base * pSize; i++){
                count[i] += count[i - 1];
            }

            #ifdef TEST
            printf("Rearranging\n");
            #endif
        }
        #pragma omp barrier

        #pragma omp for
        for(int i=0; i<pSize; i++){
            if(i != pSize - 1){
                for(int j=partitions[i]; j<partitions[i + 1]; j++){
                    digit = tgtPlace(arr[j], base, power);
                    out[count[i + digit * pSize]++] = arr[j];
                }
            }else{
                for(int j=partitions[i]; j<length; j++){
                    digit = tgtPlace(arr[j], base, power);
                    out[count[i + digit * pSize]++] = arr[j];
                }
            }
        }
    }
}

void radixSort(int *arr, size_t length, int max, unsigned int maxPartitions){
    int base = 0;
    int lL2 = (int) log2(length);
    if(lL2 < 15){
        //base = pow(2, lL2);
        base = 2<<(lL2 - 1);

    }else{
        base = 32768;
    }

    //printf("base %d\n", base);
    int maxI = (int) log(max) / log(base);


    int *count = NULL;
    int *partitions = NULL;
    int pSize;
    if(length > maxPartitions){
        count = malloc(base * maxPartitions * sizeof(int));
        partitions = malloc(maxPartitions * sizeof(int));
        pSize = maxPartitions;
    }else{
        count = malloc(base * length * sizeof(int));
        partitions = malloc(length * sizeof(int));
        pSize = length;
    }
    checkNullPtr(count);
    checkNullPtr(partitions);


    partitions[0] = 0;
    #pragma omp parallel for
    for(int i=1; i<pSize; i++){
        partitions[i] = (int) length / pSize * i;
    }

    int *arrT = NULL;
    arrT = malloc(length * sizeof(int));
    int *arrT2 = NULL;
    arrT2 = malloc(length * sizeof(int));
    checkNullPtr(arrT);
    checkNullPtr(arrT2);

    int *temp;

    #ifdef TEST
    printf("Running\n");
    #endif

    countingSort(arr, arrT, length, base, 0, count, partitions, pSize);
    if(maxI > 0){
        for(int i=1; i<maxI; i++){

            countingSort(arrT, arrT2, length, base, i, count, partitions, pSize);

            #ifdef TEST
            printArr(arrT, length);
            #endif

            temp = arrT;
            arrT = arrT2;
            arrT2 = temp;
        }
        countingSort(arrT, arr, length, base, maxI, count, partitions, pSize);
    }else{
        #ifdef TEST
        printf("Copying\n");
        #endif
        memcpy(arr, arrT, sizeof(int) * length);
    }


    //memcpy(arr, arrT, sizeof(int) * length);

    free(arrT);
    free(arrT2);
    free(count);
    free(partitions);
    //return arr;
}


//For benchmarking purposes
int compFunc (const void * a, const void * b) {
   return ( *(int*)a - *(int*)b );
}

#define ARR_LENGTH 200000
int main(){
    //int arr[] = {0, 1, 2, 4, 5, 2, 5, 3, 8, 9};
    int length = ARR_LENGTH;
    int arr[ARR_LENGTH];

    double rSSum = 0, qSSum = 0;
    int repeats=100;

    clock_t start, end;

    for(int count=0; count<repeats; count++){
        srand(time(NULL));
        for(int i=0; i<length; i++){
            arr[i] = rand();
            //printf("%d ", arr[i]);
        }
        //printf("\n");
        omp_set_num_threads(6);
        start = clock();
        radixSort(arr, length, getMax(arr, length), 6);
        end = clock();
        rSSum += (double) ((double)(end - start) / CLOCKS_PER_SEC);
        //printf("%f\n", (double) ((double)(end - start) / CLOCKS_PER_SEC) );
        checkAscending(arr, length);

        srand(time(NULL));
        for(int i=0; i<length; i++){
            arr[i] = rand();
            //printf("%d ", arr[i]);
        }

        start = clock();
        qsort(arr, length, sizeof(int), compFunc);
        end = clock();
        qSSum += (double) ((double)(end - start) / CLOCKS_PER_SEC);
        //printf("%f\n", (double) ((double)(end - start) / CLOCKS_PER_SEC) );
        checkAscending(arr, length);
    }
    printf("Radix Average: %f \nqSort Average: %f \nqSort / Radix: %f \n", rSSum / repeats, qSSum / repeats, qSSum/rSSum);



    //printArr(arr, length);


}