interatomicDist-FW-ParallelExample.cu

#include <stdio.h>
#include <cuda/std/cstdlib>

#define ATOMS_IN_MOLECULE 10
#define BLOCK_SIZE_LIMIT 32

//These bounds breaks the physical realism of molecule sizes
#define UNKNOWN_UPPER_BOUND 1.0f
#define UNKNOWN_LOWER_BOUND 0.01f
#define UNKNOWN_ACTUAL_VALUE 0.0f

//Struct representing the weight of an edge in the graph
typedef struct __align__(16)
{
    float upper;
    float lower;
    float actual;
} weight;

float generateRandFloatInWeightRange();
int generateRandInt(int maxRageFromZero);
void generateGraph(weight *adjMatrix, int atomsRequested);
void printGraphMatrix(weight *adjMatrix, int atomsHeightWidth, bool fwPass);
__global__ void floydWarshallTheChemist(weight *adjMatrix, int atomsHeightWidth);

int main(){
    //Initializing the RNG seed:
    srand(time(NULL));

    //Create matrix of size Atoms x Atoms, for a molecule made of a Atoms
    weight* moleculeGraph = (weight*)malloc((ATOMS_IN_MOLECULE * ATOMS_IN_MOLECULE) * sizeof(weight));
    generateGraph(moleculeGraph, ATOMS_IN_MOLECULE);

    printGraphMatrix(moleculeGraph, ATOMS_IN_MOLECULE, false);
    
    //Transfering the molecule/graph to the GPU
    weight* d_moleculeGraph;
    size_t moleculeSize = ATOMS_IN_MOLECULE * ATOMS_IN_MOLECULE * sizeof(weight);
    cudaMalloc(&d_moleculeGraph, moleculeSize);
    cudaMemcpy(d_moleculeGraph, moleculeGraph, moleculeSize, cudaMemcpyHostToDevice);

    //Calculating the amount of blocks, amount of threads and then launching kernel
    dim3 amountOfBlocks((ATOMS_IN_MOLECULE + BLOCK_SIZE_LIMIT - 1)/BLOCK_SIZE_LIMIT, (ATOMS_IN_MOLECULE + BLOCK_SIZE_LIMIT - 1) / BLOCK_SIZE_LIMIT);
    dim3 blockShape(BLOCK_SIZE_LIMIT, BLOCK_SIZE_LIMIT);
    floydWarshallTheChemist<<<amountOfBlocks, blockShape>>>(d_moleculeGraph, ATOMS_IN_MOLECULE);

    //Waiting for the GPU to finish executing the kernel
    cudaDeviceSynchronize();

    //Transfers the result from the GPU to the Host
    weight* processedMoleculeGraph = (weight*)malloc((ATOMS_IN_MOLECULE * ATOMS_IN_MOLECULE) * sizeof(weight));
    cudaMemcpy(processedMoleculeGraph, d_moleculeGraph, moleculeSize, cudaMemcpyDeviceToHost);
    cudaFree(d_moleculeGraph);

    printf("\n\nNewly Computed Graph: \n\n");
    printGraphMatrix(processedMoleculeGraph, ATOMS_IN_MOLECULE, true);

    free(moleculeGraph);
    free(processedMoleculeGraph);

    return 0;
}


float generateRandFloatInWeightRange(){
    int randNumPull = rand();
    while (randNumPull == 0){
        randNumPull = rand();
    }

    float randNum = ((float)(randNumPull % 101)) / 100;
    return randNum;
}

int generateRandInt(int maxRageFromZero){
    return (rand() % maxRageFromZero);
}

//Function to generate graph of molecule
void generateGraph(weight *adjMatrix, int atomsRequested){
    weight noEdgeWeight;
    noEdgeWeight.upper = 0.0f;
    noEdgeWeight.lower = 0.0f;
    noEdgeWeight.actual = 0.0f;

    weight unknownMessurementWeight;
    unknownMessurementWeight.upper = UNKNOWN_UPPER_BOUND;
    unknownMessurementWeight.lower = UNKNOWN_LOWER_BOUND;
    unknownMessurementWeight.actual = UNKNOWN_ACTUAL_VALUE;

    //Initializing the molecule/graph with unknown messurements/weights
    for(int row = 0; row < atomsRequested; row++){
        for (int column = 0; column < atomsRequested; column++)
        {
            //Since the graph is undirected, there should be no edges from a vertex to itself
            if (row == column)
            {
                //Performing 2D arr pointer arithmetic
                *(adjMatrix + (row * atomsRequested + column)) = noEdgeWeight;
            }
            else
            {
                //Setting the rest of the edges to their unknow messurement boundaries
                *(adjMatrix + (row * atomsRequested + column)) = unknownMessurementWeight;
            }
            
        }   
    }

    //For demo's sake lets say we know the actual distance between 20% of the atom pairs
    int atomPairsKnown = ((float)(atomsRequested * atomsRequested)) * 0.2f;
    printf("Atom Pairs Known: %d\n", atomPairsKnown);
    for (int i = 0; i < atomPairsKnown; i++)
    {
        int randRow = generateRandInt(atomsRequested);
        int randColumn = generateRandInt(atomsRequested);
        
        if (randRow != randColumn)
        {
            (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual = generateRandFloatInWeightRange();
            (*(adjMatrix + (randRow * atomsRequested + randColumn))).upper = (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual;
            (*(adjMatrix + (randRow * atomsRequested + randColumn))).lower = (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual;

            //Enforcing Symmetry
            (*(adjMatrix + (randColumn * atomsRequested + randRow))).actual = (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual;
            (*(adjMatrix + (randColumn * atomsRequested + randRow))).upper = (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual;
            (*(adjMatrix + (randColumn * atomsRequested + randRow))).lower = (*(adjMatrix + (randRow * atomsRequested + randColumn))).actual;
        }
        
    }

}

void printGraphMatrix(weight *adjMatrix, int atomsHeightWidth, bool fwPass){
    int atomsInTotal = 0;
    int atomsPossiblyUnknown = 0;
    int atomsNotConnected = 0;

    for(int row = 0; row < atomsHeightWidth; row++){
        for (int column = 0; column < atomsHeightWidth; column++)
        {
            float currUpper = (*(adjMatrix + (row * atomsHeightWidth + column))).upper;
            float currLower = (*(adjMatrix + (row * atomsHeightWidth + column))).lower;
            float currActual = (*(adjMatrix + (row * atomsHeightWidth + column))).actual;

            float symCurrUpper = (*(adjMatrix + (column * atomsHeightWidth + row))).upper;
            float symCurrLower = (*(adjMatrix + (column * atomsHeightWidth + row))).lower;
            float symCurrActual = (*(adjMatrix + (column * atomsHeightWidth + row))).actual;

            printf("adjMatrix[%d][%d] => (Upper: %f, Lower: %f, Actual: %f)", row, column, currUpper, currLower, currActual);

            //Checks if we by accidently break physical bounds
            if (currLower > currUpper)
            {
                printf("[BREAKS PHYSICS]");
            }
            //Checks if symmetry in the undirected graph is upheld
            if ((currUpper != symCurrUpper) || (currLower != symCurrLower) || (currActual != symCurrActual))
            {
                printf("[BREAKS SYMMETRY]");
            }
            if (currUpper == 1.0f && currLower == 0.01f && currActual == 0.0f)
            {
                //Checks wheter the graph has has a pass in the bound-smoothing algorithm
                if (fwPass == false)
                {
                    printf("[POSSIBLY UNKNOWN]");
                    atomsPossiblyUnknown++;
                }
                else
                {
                    //Means that no direct or transitive data including this edge existed
                    printf("[UNABLE TO INFER FROM AVAILABLE DATA]");
                    atomsPossiblyUnknown++;
                }
                
                
            }
            if (currUpper == 0.0f && currLower == 0.0f && currActual == 0.0f)
            {
                atomsNotConnected++;
            }
            
            
            printf("\n");

            atomsInTotal++;
        }
        
    }
    printf("Atoms in Total: %d\n", atomsInTotal);
    printf("Atoms Possibly Unknown in Total: %d\n", atomsPossiblyUnknown);
    printf("Atom Pairs Not Connected: %d\n", atomsNotConnected);
}

__global__ void floydWarshallTheChemist(weight *adjMatrix, int atomsHeightWidth){
    //We are making sure that the row (atom_i) and column (atom_j) indicies are thread-block agnostic, a.k.a making sure they work across thread blocks
    int atom_i = blockIdx.y * blockDim.y + threadIdx.y;
    int atom_j = blockIdx.x * blockDim.x + threadIdx.x;

    //Edge case & limiting threads acting out of bounds
    if (atom_i >= atomsHeightWidth || atom_j >= atomsHeightWidth)
    {
        return;
    }
    
    
    for (int atom_k = 0; atom_k < atomsHeightWidth; atom_k++)
    {   
        //Gotta wait for all threads to finish their individual computation. Otherwise we might read a false value in the next iteration
        __syncthreads();

        //Making sure that the threads are only working on one diagonal of the matrix, and then mirrors it to the other side.
        //This is to avoid two threads working on the same atom pair/edge
        if (atom_i >= atom_j)
        {
            //Skipping the nodes attempting to compute distance bounds to themselves, or those who have true data
            if ((atom_i == atom_j) || adjMatrix[atom_i * atomsHeightWidth + atom_j].actual != 0.0f)
            {
                continue;
            }


            //Computing upper bound
            float currUpperDistanceFromItoJ = adjMatrix[atom_i * atomsHeightWidth + atom_j].upper;
            float upperDistanceFromItoK = adjMatrix[atom_i * atomsHeightWidth + atom_k].upper;
            float upperDistanceFromKtoJ = adjMatrix[atom_k * atomsHeightWidth + atom_j].upper;
            float newUpper = (currUpperDistanceFromItoJ <= (upperDistanceFromItoK + upperDistanceFromKtoJ)) ? currUpperDistanceFromItoJ : (upperDistanceFromItoK + upperDistanceFromKtoJ);


            //New upper bound for the atom pair
            //Seeing wheter the current route is shorther or if the route through other atoms is shorter
            adjMatrix[atom_i * atomsHeightWidth + atom_j].upper = newUpper;

            //Upholding symmetry
            adjMatrix[atom_j * atomsHeightWidth + atom_i].upper = adjMatrix[atom_i * atomsHeightWidth + atom_j].upper;
            

        }
        


        
    }

    //Makes sure that all threads have finished updating upper bounds
    __syncthreads();

    for (int atom_k = 0; atom_k < atomsHeightWidth; atom_k++)
    {
        //Sync the threads before doing bound-smoothing, so that conflicts dont arise.
        __syncthreads();

        //Making sure that the threads are only working on one diagonal of the matrix, and then mirrors it to the other side.
        //This is to avoid two threads working on the same atom pair/edge
        if (atom_i >= atom_j)
        {
            //Skipping the nodes attempting to compute distance bounds to themselves, or those who have true data
            if ((atom_i == atom_j) || adjMatrix[atom_i * atomsHeightWidth + atom_j].actual != 0.0f)
            {
                continue;
            }

            //Computing lower bound
            float currLowerDistanceFromItoJ = (*(adjMatrix + (atom_i * atomsHeightWidth + atom_j))).lower;
            float newLower = max(currLowerDistanceFromItoJ, max((adjMatrix[atom_i * atomsHeightWidth + atom_k].lower - adjMatrix[atom_k * atomsHeightWidth + atom_j].upper), (adjMatrix[atom_k * atomsHeightWidth + atom_j].lower - adjMatrix[atom_i * atomsHeightWidth + atom_k].upper)));


            //New lower bound for the atom pair
            //Seeing wheter the current route is longest or if the route through other atoms is longer
            adjMatrix[atom_i * atomsHeightWidth + atom_j].lower = newLower;

            //Ensure that we dont break physics by surpassing the upper bounds
            if (adjMatrix[atom_i * atomsHeightWidth + atom_j].lower > adjMatrix[atom_i * atomsHeightWidth + atom_j].upper)
            {
                adjMatrix[atom_i * atomsHeightWidth + atom_j].lower = adjMatrix[atom_i * atomsHeightWidth + atom_j].upper;
            }

            //Upholding symmetry
            adjMatrix[atom_j * atomsHeightWidth + atom_i].lower = adjMatrix[atom_i * atomsHeightWidth + atom_j].lower;

        }
        

        
        
    }
    
}