Commit f2b0b20b authored by lykos98's avatar lykos98
Browse files

added openmp support, still working on optimization

parent a01f1f76

Makefile

+1 −1
Original line number Diff line number Diff line 
CC=mpicc

CFLAGS=-O3 -g

LDFLAGS=-lm

LDFLAGS=-lm -fopenmp



all: main

bb_plot.ipynb

+0 −0
Original line number Diff line number Diff line

check.py

0 → 100644
+115 −0
Original line number Diff line number Diff line 
#!/usr/bin/env python

# coding: utf-8



import matplotlib.pyplot as plt

import numpy as np

from sklearn.neighbors import NearestNeighbors



ndims = 2

k     = 100 

p     = 12



with open("bb/nodes_50_blobs_more_var.csv","r") as f:

    l = f.readlines() 



def parse_lines(l,n_dims):

    ll = [line.split(",") for line in l]

    level = np.array([ int(line[0]) for line in ll])

    owner = np.array([ int(line[1]) for line in ll])

    split_dim = np.array([ int(line[2]) for line in ll])

    split_val = np.array([ float(line[3]) for line in ll])

    box_lb = np.array([ [float(el) for el in line[4:(4+n_dims)]] for line in ll])

    box_ub = np.array([ [float(el) for el in line[4 + n_dims:]] for line in ll])

    return level, owner, split_dim, split_val, box_lb, box_ub



def plot_boxes(x,d0,d1,owner, split_dim, split_val, box_lb, box_ub, ratio = 0.7):

    from matplotlib.patches import Rectangle

    fig, ax = plt.subplots(figsize = (12 * ratio,10 * ratio))

    ax.scatter(x[:,d0],x[:,d1], s = 0.1)

    procs = np.where(owner != -1)

    for p in procs[0]:

        lbx = box_lb[p,d0]

        ubx = box_ub[p,d0]

        lby = box_lb[p,d1]

        uby = box_ub[p,d1]

        bw  = ubx - lbx

        bh  = uby - lby

        col = (np.random.rand(),np.random.rand(),np.random.rand(),0.5)

        ax.add_patch(Rectangle((lbx,lby),bw,bh, facecolor = col, label = owner[p]))

    plt.legend(loc = "lower left")

        #ax.add_patch(Rectangle((lbx,lby),2,2, facecolor = (np.random.rand(),np.random.rand(),np.random.rand(),0.3)))



def plot_planes(x,d0,d1,owner, split_dim, split_val, box_lb, box_ub, ratio=0.7):

    from matplotlib.patches import Rectangle

    fig, ax = plt.subplots(figsize = (12 * ratio,10 * ratio))

    ax.scatter(x[:,d0],x[:,d1], s = 0.1)

    procs = np.where(owner == -1)[0]

    for p in procs:

        if split_dim[p] == d0:

            line_bounds = [box_lb[p,d1],box_ub[p,d1]]

            line_coord  = split_val[p] 

            #print("vline",split_dim[p],split_dim[p], line_bounds, line_coord)

            plt.vlines(line_coord, line_bounds[0], line_bounds[1], color = "y")

        elif split_dim[p] == d1:

            line_bounds = [box_lb[p,d0],box_ub[p,d0]]

            line_coord  = split_val[p] 

            #print("hline",split_dim[p],split_dim[p], line_bounds, line_coord)

            plt.hlines(line_coord, line_bounds[0], box_ub[p,d0], color = "y")

    plt.show()





if __name__ == "__main__":

    level, owner, split_dim, split_val, box_lb, box_ub = parse_lines(l,ndims)



    #x = np.fromfile("../../robavaria/50_blobs_more_var.npy", np.float32)

    print("Loading data file")

    x = np.fromfile("./bb/ordered_data.npy", np.float64)

    x = x.reshape((x.shape[0]//ndims,ndims))



    #plot_boxes(x,0,1,owner,split_dim,split_val,box_lb,box_ub)

    #plot_planes(x,0,1,owner,split_dim,split_val,box_lb,box_ub)



    print("Loading ngbh results")

    ngbh = []

    for pp in range(p):

        ngbh.append(np.fromfile(f"./bb/rank_{pp}.ngbh", dtype = [("value","f8"),("array_idx","u8")]))

    ngbh = np.concatenate(ngbh)



    print("Searching for neighbors")

    nn = NearestNeighbors(n_jobs=-1,n_neighbors=k)



    nn.fit(x)

    dist, idx = nn.kneighbors(x)



    idx_c = ngbh["array_idx"]

    idx_c.shape

    dist_c = ngbh["value"]





    idx_c = idx_c.reshape((len(idx_c)//k,k))

    dist_c = dist_c.reshape((len(dist_c)//k,k))



    same_dist = 0

    sd_el = []

    abs_errors = 0



    print("Check")

    for i in range(len(idx_c)):

        r1 = idx[i]

        r2 = idx_c[i]

        w = np.where(r1 != r2)

        if len(w[0]) > 0:

            d1 = dist[i,w[0][0]]

            d2 = dist[i,w[0][1]]

            #print(i, w[0])

            if not np.isclose(d1,d2):

                abs_errors += 1

                same_dist += 1

                print("   Found error in ", w[0], d1, d2)

    print(f"Found {abs_errors} errors")

src/main/main.c

+15 −1
Original line number Diff line number Diff line 
#include <mpi.h>

#include <omp.h>

#include <stdio.h>

#include "../common/common.h"

#include "../tree/tree.h"



int main(int argc, char** argv) {

    #if defined (_OPENMP)

        printf("Running Hybrid (Openmp + MPI) code\n");

        int mpi_provided_thread_level;

        MPI_Init_thread( &argc, &argv, MPI_THREAD_FUNNELED, &mpi_provided_thread_level);

        if ( mpi_provided_thread_level < MPI_THREAD_FUNNELED ) 

        {

            printf("a problem arise when asking for MPI_THREAD_FUNNELED level\n");

            MPI_Finalize();

            exit( 1 );

        }

    #else

        printf("Running pure MPI code\n");

        MPI_Init(NULL, NULL);

    #endif



    char processor_name[MPI_MAX_PROCESSOR_NAME];

    int name_len;

src/tree/tree.c

+44 −47
Original line number Diff line number Diff line
@@ -16,6 +16,10 @@ 
#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <omp.h>



#define WRITE_NGBH

#define WRITE_TOP_NODES



#ifdef USE_FLOAT32

#define MPI_MY_FLOAT MPI_FLOAT
@@ -326,6 +330,9 @@ void compute_bounding_box_pointset(global_context_t *ctx, pointset_t *ps) { 


	/* compute minimum and maximum for each dimensions, store them in local bb */

	/* each processor on its own */



    /* TODO: reduction using omp directive */



	for (size_t i = 0; i < ps->n_points; ++i) 

	{

		for (size_t d = 0; d < ps->dims; ++d)
@@ -473,8 +480,6 @@ void compute_pure_global_binning(global_context_t *ctx, pointset_t *ps, 
{

	/* compute binning of data along dimension d */

	uint64_t *local_bin_count = (uint64_t *)malloc(k_global * sizeof(uint64_t));

	//MPI_DB_PRINT("%p %p %p %p %p\n", local_bin_count, global_bin_counts, ps -> data, ps -> lb_box, ps -> ub_box);

	//DB_PRINT("rank %d npoints %lu %p %p %p %p %p\n",ctx -> mpi_rank, ps -> n_points, local_bin_count, global_bin_counts, ps -> data, ps -> lb_box, ps -> ub_box);

	for (size_t k = 0; k < k_global; ++k) 

	{

		local_bin_count[k] = 0;
@@ -487,25 +492,16 @@ void compute_pure_global_binning(global_context_t *ctx, pointset_t *ps, 
	MPI_DB_PRINT("\n");

	*/







	float_t bin_w = (ps-> ub_box[d] - ps->lb_box[d]) / (float_t)k_global;



    #pragma omp parallel for

	for (size_t i = 0; i < ps->n_points; ++i) 

	{

		float_t p = ps->data[i * ps->dims + d];

		/* to prevent the border point in the box to have bin_idx == k_global causing invalid memory access */

		int bin_idx = MIN((int)((p - ps->lb_box[d]) / bin_w), k_global - 1);

		//int bin_idx = (int)((p - ps->lb_box[d]) / bin_w), k_global - 1;

		/*

		if(bin_idx < 0) 

		{

			DB_PRINT("rank %d qua %lf %lf %d %lf\n",ctx -> mpi_rank, (p - ps->lb_box[d]), (p - ps->lb_box[d]) / bin_w, bin_idx, bin_w);

			DB_PRINT("[PS BOUNDING BOX %d i have %d]: ", ctx -> mpi_rank,d);

			for(size_t d = 0; d < ps -> dims; ++d) DB_PRINT("d%d:[%lf, %lf] ",(int)d, ps -> lb_box[d], ps -> ub_box[d]); MPI_DB_PRINT("\n");

			DB_PRINT("\n");

		}

		*/

        

        #pragma omp atomic update

		local_bin_count[bin_idx]++;

	}
@@ -961,12 +957,14 @@ void build_top_kdtree(global_context_t *ctx, pointset_t *og_pointset, top_kdtree 
	}

	tree -> root = tree -> _nodes;



    #if defined(WRITE_TOP_NODES)

	MPI_DB_PRINT("Root is %p\n", tree -> root);

        if(I_AM_MASTER)

        {

            //tree_print(ctx, tree -> root);

            write_nodes_to_file(ctx, tree, "bb/nodes_50_blobs_more_var.csv");

        }

    #endif



	

	free(current_pointset.lb_box);
@@ -1042,17 +1040,9 @@ void exchange_points(global_context_t* ctx, top_kdtree_t* tree) 
	int* points_per_proc = (int*)malloc(ctx -> world_size * sizeof(int));	

	int* points_owners 	 = (int*)malloc(ctx -> dims * ctx -> local_n_points * sizeof(float_t));

	int* partition_offset = (int*)malloc(ctx -> world_size * sizeof(int));	

    for(int i = 0; i < ctx -> local_n_points; ++i)

    {

        float_t d1 = ctx -> local_data[i * ctx -> dims] + 8.33416939;

        float_t d2 = ctx -> local_data[i * ctx -> dims + 1] + 8.22858047;

        if (sqrt(d1 * d1 + d2 * d2) < 1e-5)

        {

            DB_PRINT("Rank %d found it!!! idx %d\n", ctx -> mpi_rank, i);

        }

    }



	/* compute owner */

    #pragma omp parallel for

	for(size_t i = 0; i < ctx -> local_n_points; ++i)

	{

		/* tree walk */
@@ -1130,8 +1120,7 @@ void exchange_points(global_context_t* ctx, top_kdtree_t* tree) 
	{

		ctx -> idx_start += ppp[i];

	}

    DB_PRINT("rank %d start %lu\n", ctx -> mpi_rank, ctx -> idx_start);



    //DB_PRINT("rank %d start %lu\n", ctx -> mpi_rank, ctx -> idx_start);



	/* free prv pointer */

    free(ppp);
@@ -1304,6 +1293,7 @@ void mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtre 
{

	/* local search */

	MPI_DB_PRINT("Ngbh search\n");

    #pragma omp parallel for

	for(int p = 0; p < ctx -> local_n_points; ++p)

	{

		idx_t idx = local_tree -> _nodes[p].array_idx;
@@ -1395,9 +1385,6 @@ void mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtre 


    MPI_Alltoallv(__snd_points, snd_count, snd_displ, MPI_MY_FLOAT, 

                  __rcv_points, rcv_count, rcv_displ, MPI_MY_FLOAT, ctx -> mpi_communicator); 

    HERE;



   



    float_t** rcv_work_batches = (float_t**)malloc(ctx -> world_size * sizeof(float_t*));

    for(int i = 0; i < ctx -> world_size; ++i)
@@ -1443,6 +1430,7 @@ void mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtre 


    /* compute everything */



    #pragma omp parallel for

    for(int p = 0; p < ctx -> world_size; ++p)

    {

        if(point_to_rcv_count[p] > 0)
@@ -1486,6 +1474,7 @@ void mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtre 


	for(int i = 0; i < ctx -> world_size; ++i)

    {

        #pragma omp paralell for

        for(int b = 0; b < point_to_snd_count[i]; ++b)

        {

            int idx = local_idx_of_the_point[i][b];
@@ -1506,20 +1495,30 @@ void mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtre 
    

    /* heapsort them */



    #pragma omp parallel for

    for(int i = 0; i < ctx -> local_n_points; ++i)

    {

        heap_sort(&(dp_info[i].ngbh));

    }



    #if defined(WRITE_NGBH)

    MPI_DB_PRINT("Writing ngbh to files\n");

        char ngbh_out[80];

        sprintf(ngbh_out, "./bb/rank_%d.ngbh",ctx -> mpi_rank);

        FILE* file = fopen(ngbh_out,"w");

        if(!file) 

        {

            printf("Cannot open file %s\n",ngbh_out);

        }

        else

        {

            for(int i = 0; i < ctx -> local_n_points; ++i)

            {

                fwrite(dp_info[i].ngbh.data, sizeof(heap_node), k, file);

            }

            fclose(file);

        }

    #endif



    MPI_Barrier(ctx -> mpi_communicator);
@@ -1652,12 +1651,12 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) 


	if (ctx->mpi_rank == 0) 

	{

		data = read_data_file(ctx, "../norm_data/50_blobs_more_var.npy", MY_TRUE);

        ctx->dims = 2;

		//data = read_data_file(ctx, "../norm_data/50_blobs_more_var.npy", MY_TRUE);

        //ctx->dims = 2;

		//data = read_data_file(ctx, "../norm_data/50_blobs.npy", MY_TRUE);

		// std_g0163178_Me14_091_0000

		//data = read_data_file(ctx,"../norm_data/std_LR_091_0001",MY_TRUE);

        //ctx->dims = 5;

		data = read_data_file(ctx,"../norm_data/std_LR_091_0001",MY_TRUE);

        ctx->dims = 5;



		// ctx -> n_points = 48*5*2000;

		ctx->n_points = ctx->n_points / ctx->dims;
@@ -1732,11 +1731,9 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) 
	kdtree_v2_init( &local_tree, ctx -> local_data, ctx -> local_n_points, (unsigned int)ctx -> dims);

	int k = 100;



	MPI_DB_PRINT("uu\n");

	datapoint_info_t* dp_info = (datapoint_info_t*)malloc(ctx -> local_n_points * sizeof(datapoint_info_t));			

	build_local_tree(ctx, &local_tree);



	MPI_DB_PRINT("Mi pianto qua\n");

	mpi_ngbh_search(ctx, dp_info, &tree, &local_tree, ctx -> local_data, k);