Wsrtjtyk

I need to paralelize one funcion without race conditions nor false sharing. I have tried many ways but I could not achieve that yet. The function is:

__inline static

void calculateClusterCentroIDs(int numCoords, int numObjs, int numClusters, float * dataSetMatrix, int * clusterAssignmentCurrent, float *clustersCentroID) {

    int * clusterMemberCount = (int *) calloc (numClusters,sizeof(float));



    // sum all points

    // for every point

    for (int i = 0; i < numObjs; ++i) {

        // which cluster is it in?

        int activeCluster = clusterAssignmentCurrent[i];



        // update count of members in that cluster

        ++clusterMemberCount[activeCluster];



        // sum point coordinates for finding centroid

        for (int j = 0; j < numCoords; ++j)

            clustersCentroID[activeCluster*numCoords + j] += dataSetMatrix[i*numCoords + j];

    }





    // now divide each coordinate sum by number of members to find mean/centroid

    // for each cluster

    for (int i = 0; i < numClusters; ++i) {

        if (clusterMemberCount[i] != 0)

            // for each coordinate

            for (int j = 0; j < numCoords; ++j)

                clustersCentroID[i*numCoords + j] /= clusterMemberCount[i];  /// XXXX will divide by zero here for any empty clusters!

    }

Any idea how could I achieve that?

Thank you.

This is pretty straight forward

// sum all points

// for every point

for (int i = 0; i < numObjs; ++i) {

    // which cluster is it in?

    int activeCluster = clusterAssignmentCurrent[i];



    // update count of members in that cluster

    ++clusterMemberCount[activeCluster];



    // sum point coordinates for finding centroid

    #pragma omp parallel for

    for (int j = 0; j < numCoords; ++j)

        clustersCentroID[activeCluster*numCoords + j] += dataSetMatrix[i*numCoords + j];

}

The inner loop is perfectly fine to parallelize as all writes happen to different elements of clustersCentroID. You can safely assume that the default schedule will not exhibit significant false sharing, it typically has large-enough chunks. Just don't try something like schedule(static,1).

The outer loop is not as easy to parallelize. You could either use a reduction on clusterMemberCount and clusterMemberCount, or do something like:

#pragma omp parallel // note NO for

for (int i = 0; i < numObjs; ++i) {

    int activeCluster = clusterAssignmentCurrent[i];

    // ensure that exactly one thread works on each cluster

    if (activeCluster % omp_num_threads() != omp_get_thread_num()) continue;

Only do this if the simple solution does not yield sufficient performance.

The other loop is simple as well

#pragma omp parallel for

for (int i = 0; i < numClusters; ++i) {

    if (clusterMemberCount[i] != 0)

        // for each coordinate

        for (int j = 0; j < numCoords; ++j)

            clustersCentroID[i*numCoords + j] /= clusterMemberCount[i];

}

Again, data access is perfectly isolated both in terms of correctness and, except for edge cases, false sharing.

You should give an order of magnitude for the expected values of numCoords, numObjs and numClusters as the optimal way to parallelize depends on that. Especially, numCoords is important to see if parallelizing/vectorizing the inner loop over coordinates makes sense; like, are you taking of 3D coordinates or 1000 dimensions?

Another attempt with the drawback of an if statement in the first loop (bad for performance), static schedule (possible load unbalance) but each thread incrementing contiguous parts of clusterMemberCount and clustersCentroID which limits the risk of false sharing.

#ifdef _OPENMP

   #include <omp.h>

#else

   #define omp_get_num_threads() 1

   #define omp_get_thread_num() 0

#endif





__inline static

void calculateClusterCentroIDs(int numCoords, int numObjs, int numClusters, float * dataSetMatrix, int * clusterAssignmentCurrent, float *clustersCentroID) {

    int * clusterMemberCount = (int *) calloc (numClusters,sizeof(float));

    // sum all points

    // for every point



    #pragma omp parallel

    {

        int nbOfThreads = omp_get_num_threads();

        int thisThread = omp_get_thread_num();

        // Schedule for the first step : process only cluster with ID in the [from , to[ range

        int clustFrom = (thisThread*numClusters)/nbOfThreads;

        int clustTo   = (thisThread+1 == nbOfThreads) ? numClusters : ((thisThread+1)*numClusters)/nbOfThreads;



        // Each thread will loop through all values of numObjs but only process them depending on activeCluster

        // The loop is skipped only if the thread was assigned no cluster

        if (clustTo>clustFrom){

            for (int i = 0; i < numObjs; ++i) {

                // which cluster is it in?

                int activeCluster = clusterAssignmentCurrent[i];



                if (activeCluster>=clustFrom && activeCluster<clustTo){

                    // update count of members in that cluster

                    ++clusterMemberCount[activeCluster];



                    // sum point coordinates for finding centroid

                    for (int j = 0; j < numCoords; ++j)

                        clustersCentroID[activeCluster*numCoords + j] += dataSetMatrix[i*numCoords + j];

                }

            }

        }



        #pragma omp barrier



        // now divide each coordinate sum by number of members to find mean/centroid

        // for each cluster

        #pragma omp for // straightforward

        for (int i = 0; i < numClusters; ++i) {

            if (clusterMemberCount[i] != 0)

                // for each coordinate

                for (int j = 0; j < numCoords; ++j)

                    clustersCentroID[i*numCoords + j] /= clusterMemberCount[i];  /// XXXX will divide by zero here for any empty clusters!

        }

    }

    free(clusterMemberCount);

}

Adding to my comment: ++clusterMemberCount[activeCluster] forms a histogram and is problematic when two threads try to update same item (or adjacent item sharing a cache line). This needs to be taken out of the loop either as a sequential part, or must be parallelized by having a separate copy of the histogram for each thread, then combined.

You can easily separate this part from the first parallel loop.

// Make the histogram

for (int i = 0; i < numObjs; ++i) {

    int activeCluster = clusterAssignmentCurrent[i];

    ++clusterMemberCount[activeCluster];

}

Then process everything exploiting parallelism

// parallel processing

#pragma omp parallel for

for (int i = 0; i < numObjs; ++i) {

    int activeCluster = clusterAssignmentCurrent[i];

    for (int j = 0; j < numCoords; ++j)

        clustersCentroID[activeCluster*numCoords + j] += dataSetMatrix[i*numCoords + j];

}

The second time there's a potential false sharing is when numCoords * sizeof(clustersCentroID[0]) % 64 != 0 assuming 64-byte cache line. This can be mitigated by overallocating clustersCentroID to full multiples of 64 bytes.

// Loop for numCoords, but index by numCoordsX

for (int j = 0; j < numCoords; ++j)

    clustersCentroID[activeCluster*numCoordsX + j] += dataSetMatrix[i*numCoords + j];