H1 optimization (!1) · Merge requests · Luca Tornatore / dADP

src/adp/adp.c

+394 −86

Changes for src/adp/adp.c: 394 added lines, 86 removed lines.
Original line number	Diff line number	Diff line
		@@ -118,8 +118,8 @@ idx_t get_j_ksel_idx(global_context_t* ctx, idx_t j, idx_t ksel, MPI_Win exposed
		}


		void compute_density_kstarnn_rma_v2(global_context_t* ctx, const float_t d, int verbose){

		void compute_density_kstarnn_rma_v2(global_context_t* ctx, const float_t d, int verbose)
		{
		/*
		* Point density computation:
		* args:
		@@ -510,8 +510,6 @@ datapoint_info_t find_possibly_halo_datapoint_rma(global_context_t* ctx, idx_t i
		else
		{
		datapoint_info_t tmp_dp;
		#pragma omp critical
		{
		idx_t i = idx - ctx -> rank_idx_start[owner];
		MPI_Request request;
		MPI_Status status;
		@@ -520,8 +518,6 @@ datapoint_info_t find_possibly_halo_datapoint_rma(global_context_t* ctx, idx_t i
		i * sizeof(datapoint_info_t), sizeof(datapoint_info_t), MPI_BYTE, win_datapoints, &request);
		MPI_Wait(&request, MPI_STATUS_IGNORE);

		}

		return tmp_dp;
		}
		}
		@@ -648,21 +644,43 @@ lock_t h1_lock_free(global_context_t* ctx, MPI_Win lock_window, int owner, idx_t

		}

		void enqueue_center_removal(center_removal_queue_t* queue,
		center_removal_t element)
		{
		int default_incr_size = 100;
		if(queue -> count >= queue -> size)
		{
		queue -> size += default_incr_size;
		queue -> data = realloc(queue -> data, queue -> size * sizeof(center_removal_t));
		}
		queue -> data[queue -> count] = element;
		queue -> count++;
		}

		int compare_removal_by_target(const void* a, const void* b)
		{
		idx_t arg1 = ((const center_removal_t )a).target_id;
		idx_t arg2 = ((const center_removal_t )b).target_id;
		return (arg1 > arg2) - (arg1 < arg2);
		}

		clusters_t Heuristic1(global_context_t *ctx)
		{
		/*
		/**
		* Heurisitc 1, from paper of Errico, Facco, Laio & Rodriguez
		* ( https://doi.org/10.1016/j.ins.2021.01.010 )
		*/
		**/

		datapoint_info_t* dp_info = ctx -> local_datapoints;
		idx_t n = ctx -> local_n_points;

		struct timespec start_tot, finish_tot;
		double elapsed_tot;
		double elapsed_time;

		TIME_DEF;

		TIME_START;
		lu_dynamic_array_t all_centers, removed_centers, actual_centers, max_rho;

		lu_dynamic_array_allocate(&all_centers);
		@@ -679,7 +697,7 @@ clusters_t Heuristic1(global_context_t *ctx)
		1, MPI_INFO_NULL, ctx -> mpi_communicator, &win_datapoints);
		MPI_Win_fence(0, win_datapoints);
		MPI_Win_lock_all(0, win_datapoints);

		//
		#if !defined(THREAD_FUNNELED)
		#pragma omp parallel for
		#endif
		@@ -715,45 +733,166 @@ clusters_t Heuristic1(global_context_t *ctx)
		}
		}

		/*
		/**
		* optimized version
		*
		* Generate a mask that keeps track of the point has been eliminating the
		* point considered. Each thread updates this mask, then after the procedure
		* ends, center, removed centers, and max_rho arrays are populated
		*/

		lock_t* lock_array = (lock_t)MY_MALLOC(n sizeof(lock_t));
		*
		* optimized v2 use a queue of center removal and then exchange them
		**/
		heap_node* to_remove_mask = (heap_node)MY_MALLOC(nsizeof(heap_node));

		for(idx_t p = 0; p < n; ++p)
		{
		to_remove_mask[p].array_idx = MY_SIZE_MAX;
		to_remove_mask[p].value = 9999999;
		lock_array[p] = LOCK_FREE;
		to_remove_mask[p].value = -9999999;
		}

		// sort by density

		qsort(dp_info_ptrs, n, sizeof(datapoint_info_t*), cmpPP);

		/**
		* workflow
		* find the elements we need to eliminate
		* compact them to a single array
		* send/recieve
		**/

		MPI_Win win_to_remove_mask;
		MPI_Win_create(to_remove_mask, n * sizeof(heap_node), 1, MPI_INFO_NULL, ctx -> mpi_communicator, &win_to_remove_mask);
		MPI_Win_fence(0, win_to_remove_mask);
		center_removal_queue_t* removal_queues = (center_removal_queue_t)MY_MALLOC(n sizeof(center_removal_queue_t));
		omp_lock_t* lock_array = (omp_lock_t)MY_MALLOC(n sizeof(omp_lock_t));

		MPI_Win win_locks;
		MPI_Win_create(lock_array, n * sizeof(lock_t), sizeof(lock_t), MPI_INFO_NULL, ctx -> mpi_communicator, &win_locks);
		MPI_Win_fence(0, win_locks);
		//zero out all queues

		#ifdef EXP_H1
		MPI_Win_lock_all(0, win_to_remove_mask);
		MPI_Win_lock_all(0, win_locks);
		#endif
		for(idx_t i = 0; i < n; ++i)
		{
		removal_queues[i].count = 0;
		removal_queues[i].size = 0;
		removal_queues[i].data = NULL;
		omp_init_lock(lock_array + i);
		}

		#ifdef EXP_H1
		printf("Using experimental h1\n");
		#endif
		elapsed_time = TIME_STOP;
		LOG_WRITE("Putative centers", elapsed_time);

		#if !defined(THREAD_FUNNELED)
		TIME_START;



		//#define EXP_CENTER_PRUNING
		#if defined(EXP_CENTER_PRUNING)
		int all_have_finished = 0;
		int i_have_finished = 0;

		datapoint_info_t* tmp_datapoints = (datapoint_info_t)MY_MALLOC(ctx -> local_n_points sizeof(datapoint_info_t));

		// use strategy similar to the one employed in density estimation
		// vertical loop, first on k the on n
		// request points
		// then process the points


		MPI_Win_fence(MPI_MODE_NOPUT, win_datapoints);
		//MPI_Win_lock_all(0, win_datapoints);
		MPI_DB_PRINT("Using vertical loop on centers pruning\n");
		int kmax = ctx -> k;
		for(idx_t j = 1; j < kmax; ++j)
		{
		i_have_finished = 1;

		// request points
		#pragma omp parallel for schedule(dynamic)
		for(idx_t p = 0; p < n; ++p)
		{

		datapoint_info_t i_point = *(dp_info_ptrs[p]);
		if(j < i_point.kstar + 1)
		{
		#pragma omp atomic write
		i_have_finished = 0;

		idx_t jidx = i_point.ngbh[j].array_idx;
		int owner = foreign_owner(ctx, jidx);
		/* find if datapoint is halo or not */
		idx_t jpos = jidx - ctx -> rank_idx_start[owner];
		if(owner == ctx -> mpi_rank)
		{
		tmp_datapoints[p] = ctx -> local_datapoints[jpos];
		}
		else
		{
		MPI_Get(tmp_datapoints + p, sizeof(datapoint_info_t), MPI_BYTE, owner,
		jpos * sizeof(datapoint_info_t), sizeof(datapoint_info_t), MPI_BYTE, win_datapoints);
		}
		}

		}

		//MPI_Win_flush_all(win_datapoints);
		MPI_Win_fence(MPI_MODE_NOPUT ^ MPI_MODE_NOSTORE, win_datapoints);
		// wait for all comms to happen

		#pragma omp parallel for schedule(dynamic)
		for(idx_t p = 0; p < n; ++p)
		{
		datapoint_info_t i_point = *(dp_info_ptrs[p]);
		if(j < i_point.kstar + 1)
		{
		datapoint_info_t j_point = tmp_datapoints[p];
		idx_t jidx = j_point.array_idx;

		if(j_point.is_center && i_point.g > j_point.g)
		{
		// set the lock at position i
		// actually is the p-th point
		int owner = foreign_owner(ctx, jidx);
		//if local process it
		if(owner == ctx -> mpi_rank)
		{
		//acquire the lock
		idx_t jpos = jidx - ctx -> idx_start;
		omp_set_lock(lock_array + jpos);
		if(i_point.g > to_remove_mask[jpos].value)
		{
		to_remove_mask[jpos].array_idx = i_point.array_idx;
		to_remove_mask[jpos].value = i_point.g;
		}
		omp_unset_lock(lock_array + jpos);
		}
		////otherwise enqueue for sending
		else
		{
		center_removal_t element = {.rank = owner, .source_id = i_point.array_idx,
		.target_id = j_point.array_idx, .source_density = i_point.g};
		enqueue_center_removal(removal_queues + p, element);
		}
		}

		}

		}

		//MPI_Allreduce(&i_have_finished, &all_have_finished, 1, MPI_INT, MPI_LAND, ctx -> mpi_communicator);
		//if(all_have_finished) break;
		if(i_have_finished) break;
		}

		MPI_Win_unlock_all(win_datapoints);
		MPI_Win_fence(0, win_datapoints);
		free(tmp_datapoints);

		#else
		//MPI_Win_unlock_all(win_datapoints);
		//MPI_Win_fence(MPI_MODE_NOPUT, win_datapoints);

		idx_t n_foreign_req = 0;
		idx_t n_local_req = 0;
		float_t elapsed_proc = 0.;;

		#if !defined(THREAD_FUNNELED)
		#pragma omp parallel for schedule(dynamic) reduction(+:elapsed_proc, n_local_req, n_foreign_req)
		#endif
		for(idx_t p = 0; p < n; ++p)
		{
		@@ -762,87 +901,234 @@ clusters_t Heuristic1(global_context_t *ctx)
		for(idx_t j = 1; j < i_point.kstar + 1; ++j)
		{
		idx_t jidx = i_point.ngbh[j].array_idx;

		struct timespec __start, __end;
		clock_gettime(CLOCK_MONOTONIC,&__start);
		datapoint_info_t j_point = find_possibly_halo_datapoint_rma(ctx, jidx, win_datapoints);
		clock_gettime(CLOCK_MONOTONIC,&__end);
		elapsed_proc += (double)(__end.tv_sec - __start.tv_sec) + (__end.tv_nsec - __start.tv_nsec)/1e9;

		if(j_point.is_center && i_point.g > j_point.g)
		int owner = foreign_owner(ctx, jidx);
		if(owner == ctx -> mpi_rank)
		{
		/*
		*
		* TODO: Implement it without this but using private locks
		* use an array of locks, and compare and swap to actually gain control of the thing
		*
		* */
		n_local_req++;
		}
		else
		{
		n_foreign_req++;
		}

		#ifdef EXP_H1
		#pragma omp critical (h1_exp)
		if(j_point.is_center && i_point.g > j_point.g)
		{
		int owner = foreign_owner(ctx, jidx);
		idx_t jpos = jidx - ctx -> rank_idx_start[owner];
		// set the lock at position i
		// actually is the p-th point
		// if local process it
		if(owner == ctx -> mpi_rank)
		{
		// acquire the lock
		idx_t jpos = jidx - ctx -> idx_start;
		omp_set_lock(lock_array + jpos);
		if(i_point.g > to_remove_mask[jpos].value)
		{
		to_remove_mask[jpos].array_idx = i_point.array_idx;
		to_remove_mask[jpos].value = i_point.g;
		}
		omp_unset_lock(lock_array + jpos);
		}
		// otherwise enqueue for sending
		else
		{
		center_removal_t element = {.rank = owner, .source_id = i_point.array_idx,
		.target_id = j_point.array_idx, .source_density = i_point.g};
		enqueue_center_removal(removal_queues + p, element);
		}
		}
		}
		}
		//MPI_Win_fence(MPI_MODE_NOPUT, win_datapoints);

		lock_t state = LOCK_FREE;
		// MPI_Barrier(ctx -> mpi_communicator);
		// DB_PRINT("Rank %d: foreign requests points %lu out of %lu -> fraction %.2lf time %.2lfs\n",
		// ctx -> mpi_rank, n_foreign_req, n_local_req + n_foreign_req, (float)n_foreign_req/(float)n_local_req, elapsed_proc);
		// MPI_Barrier(ctx -> mpi_communicator);
		#endif

		state = h1_lock_acquire(ctx, win_locks, owner, jpos, state);
		//assemble arrays into a single buffer

		heap_node mask_element;
		MPI_Request request;

		MPI_Rget(&mask_element, sizeof(heap_node), MPI_BYTE,
		owner, jpos * sizeof(heap_node), sizeof(heap_node), MPI_BYTE, win_to_remove_mask, &request);
		MPI_Wait(&request, MPI_STATUS_IGNORE);
		elapsed_time = TIME_STOP;
		LOG_WRITE("Finding centers to prune", elapsed_time);
		TIME_START;

		int flag = mask_element.array_idx == MY_SIZE_MAX;
		if(flag \|\| i_point.g > mask_element.value )
		idx_t tot_removal = 0;
		for(idx_t p = 0; p < n; ++p)
		{
		heap_node tmp_mask_element = {.array_idx = i_point.array_idx, .value = i_point.g};
		MPI_Request request;
		MPI_Rput(&tmp_mask_element, sizeof(heap_node), MPI_BYTE, owner,
		jpos*sizeof(heap_node), sizeof(heap_node), MPI_BYTE, win_to_remove_mask, &request);
		MPI_Wait(&request, MPI_STATUS_IGNORE);
		tot_removal += removal_queues[p].count;
		}

		idx_t buffer_idx = 0;
		center_removal_t* removal_buffer = (center_removal_t)MY_MALLOC(tot_removal sizeof(center_removal_t));

		for(idx_t p = 0; p < n; ++p)
		{
		if(removal_queues[p].count > 0)
		{
		memcpy(removal_buffer + buffer_idx, removal_queues[p].data, removal_queues[p].count * sizeof(center_removal_t));
		buffer_idx += removal_queues[p].count;
		}
		}

		//sort by array idx (it sorts also by rank)

		qsort(removal_buffer, tot_removal, sizeof(center_removal_t), compare_removal_by_target);

		state = h1_lock_free(ctx, win_locks, owner, jpos, state);
		// remove
		for(int i = 1; i < tot_removal; ++i)
		{
		if(removal_buffer[i - 1].rank > removal_buffer[i].rank)
		printf("Unsorted removal buffer");
		}
		#else
		#pragma omp critical (centers_elimination)

		//prepare for the sendrcv

		int* recv_counts = (int)MY_MALLOC(ctx -> world_size sizeof(int));
		int* send_counts = (int)MY_MALLOC(ctx -> world_size sizeof(int));

		int* recv_displs = (int)MY_MALLOC(ctx -> world_size sizeof(int));
		int* send_displs = (int)MY_MALLOC(ctx -> world_size sizeof(int));

		//zero_out_all

		for(int i = 0; i < ctx -> mpi_rank; ++i)
		{
		int owner = foreign_owner(ctx, jidx);
		idx_t jpos = jidx - ctx -> rank_idx_start[owner];
		recv_counts[i] = 0;
		send_counts[i] = 0;
		recv_displs[i] = 0;
		send_displs[i] = 0;
		}

		MPI_Win_lock(MPI_LOCK_EXCLUSIVE, owner, 0, win_to_remove_mask);
		heap_node mask_element;
		MPI_Request request;
		MPI_Rget(&mask_element, sizeof(heap_node), MPI_BYTE,
		owner, jpos * sizeof(heap_node), sizeof(heap_node), MPI_BYTE, win_to_remove_mask, &request);
		MPI_Wait(&request, MPI_STATUS_IGNORE);
		for(idx_t i = 0; i < tot_removal; ++i)
		{
		int rank_idx = removal_buffer[i].rank;
		send_counts[rank_idx]++;
		}

		// all to all to recieve counts

		MPI_Alltoall(send_counts, 1, MPI_INT,
		recv_counts, 1, MPI_INT, ctx -> mpi_communicator);

		int flag = mask_element.array_idx == MY_SIZE_MAX;
		if(flag \|\| i_point.g > mask_element.value )
		// compute displacements

		for(int i = 1; i < ctx -> world_size; ++i)
		{
		heap_node tmp_mask_element = {.array_idx = i_point.array_idx, .value = i_point.g};
		MPI_Request request;
		MPI_Rput(&tmp_mask_element, sizeof(heap_node), MPI_BYTE, owner,
		jpos*sizeof(heap_node), sizeof(heap_node), MPI_BYTE, win_to_remove_mask, &request);
		MPI_Wait(&request, MPI_STATUS_IGNORE);
		recv_displs[i] = recv_displs[i - 1] + recv_counts[i - 1];
		send_displs[i] = send_displs[i - 1] + send_counts[i - 1];
		}

		idx_t tot_recv_counts = 0;

		// count how many elements to recieve
		// MPI_DB_PRINT("Using centers elimination queue experiment\n");

		for(int i = 0; i < ctx -> world_size; ++i) tot_recv_counts += recv_counts[i];
		/*
		if(ctx -> mpi_rank == 0){
		for(int i = 0; i < ctx -> world_size; ++i){
		DB_PRINT("%d mpi rank recv_count %d to %d\n", ctx -> mpi_rank, recv_counts[i], i);
		DB_PRINT("%d mpi rank send_count %d to %d\n", ctx -> mpi_rank, send_counts[i], i);
		}
		}
		*/
		/DB_PRINT("rank %d: %lu recv counts\n", ctx -> mpi_rank, tot_recv_counts);/

		// change dimensions to bytes
		//

		//comm matrices

		/* DEBUG PRINT
		* COMMUNICATION MATRICES */

		/*
		int* all_send = (int)MY_MALLOC(ctx -> world_size ctx -> world_size * sizeof(int));
		int* all_recv = (int)MY_MALLOC(ctx -> world_size ctx -> world_size * sizeof(int));

		MPI_Gather(send_counts, ctx -> world_size, MPI_INT,
		all_send, ctx -> world_size , MPI_INT, 0, ctx -> mpi_communicator);

		MPI_Win_unlock(owner, win_to_remove_mask);
		MPI_Gather(recv_counts, ctx -> world_size, MPI_INT,
		all_recv, ctx -> world_size , MPI_INT, 0, ctx -> mpi_communicator);

		if(I_AM_MASTER)
		{
		printf("Send matrix\n");
		for(int i = 0; i < ctx -> world_size; ++i)
		{
		for(int j = 0; j < ctx -> world_size; ++j)
		{
		printf("%4d ", all_send[i*ctx -> world_size + j]);
		}
		#endif
		printf("\n");
		}

		printf("Recv matrix\n");
		for(int i = 0; i < ctx -> world_size; ++i)
		{
		for(int j = 0; j < ctx -> world_size; ++j)
		{
		printf("%4d ", all_recv[i*ctx -> world_size + j]);
		}
		printf("\n");
		}
		}

		#ifdef EXP_H1
		MPI_Win_unlock_all(win_to_remove_mask);
		MPI_Win_unlock_all(win_locks);
		#endif
		free(all_recv);
		free(all_send);

		*/


		for(int i = 0; i < ctx -> world_size; ++i)
		{
		recv_counts[i] = recv_counts[i] * sizeof(center_removal_t);
		send_counts[i] = send_counts[i] * sizeof(center_removal_t);
		recv_displs[i] = recv_displs[i] * sizeof(center_removal_t);
		send_displs[i] = send_displs[i] * sizeof(center_removal_t);
		}

		//allocate buffer to recieve center elminiations


		center_removal_t* recv_removals = (center_removal_t)MY_MALLOC(tot_recv_counts sizeof(center_removal_t));

		// all to all
		MPI_Alltoallv(removal_buffer, send_counts, send_displs, MPI_BYTE,
		recv_removals, recv_counts, recv_displs, MPI_BYTE, ctx -> mpi_communicator);

		// merge into the mask

		elapsed_time = TIME_STOP;
		LOG_WRITE("Communicating eliminations", elapsed_time);
		TIME_START;

		#pragma omp parallel for schedule(dynamic)
		for(idx_t i = 0; i < tot_recv_counts; ++i)
		{
		idx_t el_pos = recv_removals[i].target_id - ctx -> idx_start;
		int owner = foreign_owner(ctx, recv_removals[i].target_id);
		if(owner != ctx -> mpi_rank){
		printf("Error here\n");
		exit(1);
		}
		omp_set_lock(lock_array + el_pos);
		if(recv_removals[i].source_density > to_remove_mask[el_pos].value)
		{
		to_remove_mask[el_pos].array_idx = recv_removals[i].source_id;
		to_remove_mask[el_pos].value = recv_removals[i].source_density;
		}
		omp_unset_lock(lock_array + el_pos);
		}

		MPI_Win_fence(0, win_to_remove_mask);
		MPI_Win_fence(0, win_locks);
		MPI_Barrier(ctx -> mpi_communicator);

		/* populate the usual arrays */
		for(idx_t p = 0; p < all_centers.count; ++p)
		@@ -888,11 +1174,29 @@ clusters_t Heuristic1(global_context_t *ctx)
		}
		}

		MPI_Win_free(&win_to_remove_mask);
		free(to_remove_mask);

		MPI_Win_free(&win_locks);
		for(idx_t i = 0; i < n; ++i)
		{
		removal_queues[i].count = 0;
		removal_queues[i].size = 0;
		free(removal_queues[i].data);
		omp_destroy_lock(lock_array+ i);
		}

		free(removal_queues);
		free(removal_buffer);
		free(send_counts);
		free(send_displs);
		free(recv_counts);
		free(recv_displs);
		free(lock_array);
		free(recv_removals);

		elapsed_time = TIME_STOP;
		LOG_WRITE("Merging", elapsed_time);

		TIME_START;

		int n_centers = (int)actual_centers.count;
		int tot_centers;
		@@ -900,9 +1204,11 @@ clusters_t Heuristic1(global_context_t *ctx)

		MPI_DB_PRINT("Found %d temporary centers\n", tot_centers);

		/* bring on master all centers
		/**
		* bring on master all centers
		* order them in ascending order of density,
		* then re-scatter them around to get unique cluster labels */
		* then re-scatter them around to get unique cluster labels
		**/

		center_t* private_centers_buffer = (center_t)MY_MALLOC(actual_centers.count sizeof(center_t));
		center_t* global_centers_buffer = (center_t)MY_MALLOC(tot_centers sizeof(center_t));
		@@ -1055,6 +1361,8 @@ clusters_t Heuristic1(global_context_t *ctx)
		free(ks);
		#endif

		elapsed_time = TIME_STOP;
		LOG_WRITE("Cluster assign", elapsed_time);

		free(actual_centers.data);
		actual_centers.size = tot_centers;

src/adp/adp.h

+15 −0

Changes for src/adp/adp.h: 15 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -45,6 +45,21 @@ typedef struct merge_t
		float_t density;
		} merge_t;

		typedef struct center_removal_t
		{
		int rank;
		idx_t source_id;
		idx_t target_id;
		float_t source_density;
		} center_removal_t;

		typedef struct center_removal_queue_t
		{
		center_removal_t* data;
		idx_t count;
		idx_t size;
		} center_removal_queue_t;



		void compute_density_kstarnn_rma(global_context_t* ctx, const float_t d, int verbose);

src/common/common.c

+4 −2

Changes for src/common/common.c: 4 added lines, 2 removed lines.
Original line number	Diff line number	Diff line
		@@ -21,6 +21,8 @@ void get_context(global_context_t* ctx)
		ctx -> __local_heap_buffers = NULL;
		ctx -> input_data_in_float32 = -1;
		ctx -> dims = 0;
		ctx -> k = 300;
		ctx -> z = 3;
		}

		void get_dataset_diagnostics(global_context_t* ctx, float_t* data)
		@@ -42,7 +44,7 @@ void get_dataset_diagnostics(global_context_t* ctx, float_t* data)
		for(int j = 0; j < ctx -> dims; ++j) pvt_mean[j] = 0.;

		#pragma omp for
		for(int i = 0; i < ctx -> n_points; ++i)
		for(idx_t i = 0; i < ctx -> n_points; ++i)
		{
		int j = 0;
		for(j = 0; j < jmax; j+=4)
		@@ -79,7 +81,7 @@ void get_dataset_diagnostics(global_context_t* ctx, float_t* data)
		for(int j = 0; j < ctx -> dims; ++j) pvt_var[j] = 0.;

		#pragma omp for
		for(int i = 0; i < ctx -> n_points; ++i)
		for(idx_t i = 0; i < ctx -> n_points; ++i)
		{
		int j = 0;
		for(j = 0; j < jmax; j+=4)

src/common/common.h

+5 −2

Changes for src/common/common.h: 5 added lines, 2 removed lines.
Original line number	Diff line number	Diff line
		@@ -45,9 +45,11 @@
		#define MY_TRUE 1
		#define MY_FALSE 0

		#define CHECK_ALLOCATION(x) if(!x){printf("[!!!] %d rank encountered failed allocation at line %s \n", ctx -> mpi_rank, __LINE__ ); exit(1);};
		#define HERE printf("%d in file %s reached line %d\n", ctx -> mpi_rank, __FILE__, __LINE__); MPI_Barrier(ctx -> mpi_communicator);

		#define CHECK_ALLOCATION_NO_CTX(x) if(!x){printf("[!!!] Failed allocation at line %d \n", __LINE__ ); exit(1);}
		#define CHECK_ALLOCATION(x) if(!x){printf("[!!!] %d rank encountered failed allocation: %s at line %s \n", ctx -> mpi_rank, __FILE__, __LINE__ ); exit(1);};

		#define CHECK_ALLOCATION_NO_CTX(x) if(!x){printf("[!!!] Failed allocation: %s at line %d \n", __FILE__, __LINE__ ); exit(1);}
		#define MY_MALLOC(n) ({void* p = calloc(n,1); CHECK_ALLOCATION_NO_CTX(p); p; })

		#define DB_PRINT(...) printf(__VA_ARGS__)
		@@ -149,6 +151,7 @@ typedef struct global_context_t
		int mpi_rank; //rank of the processor
		uint32_t dims; //number of dimensions of the dataset
		idx_t k; //number of neighbors
		float_t z; //z parameter
		float_t* local_data; //pointer to the dataset stored into the node
		float_t* lb_box; //bounding box of the dataset
		float_t* ub_box; //bounding box of the dataset

src/main/main.c

+60 −21

Changes for src/main/main.c: 60 added lines, 21 removed lines.
Original line number	Diff line number	Diff line
		@@ -44,6 +44,8 @@ struct option long_options[] =
		{"in-dims" , required_argument, 0, 'd'},
		{"out-data", optional_argument, 0, 'o'},
		{"out-assignment", optional_argument, 0, 'a'},
		{"kngbh", optional_argument, 0, 'k'},
		{"z", optional_argument, 0, 'z'},
		{"help", optional_argument, 0, 'h'},
		{0, 0, 0, 0}
		};
		@@ -59,7 +61,9 @@ const char* help = "Distributed Advanced Density Peak\n"\
		" mpi tasks and datapoints are ordered default is `out_data` \n"\
		" -a --out-assignment (optional) output path for the cluster assignment output ranges [0 ... Nc - 1]\n"\
		" for core points halo points have indices [-Nc ... -1] conversion\n"\
		" of idx for an halo point is cluster_idx = -halo_idx - 1, default is `out_assignment`\n";
		" of idx for an halo point is cluster_idx = -halo_idx - 1, default is `out_assignment`\n"\
		" -k --kngbh (optional) number of nearest neighbors to compute\n"\
		" -z --z (optional) number of nearest neighbors to compute\n";

		void parse_args(global_context_t* ctx, int argc, char** argv)
		{
		@@ -70,7 +74,7 @@ void parse_args(global_context_t* ctx, int argc, char** argv)
		snprintf(ctx -> output_assignment_file, DEFAULT_STR_LEN, "%s", OUT_CLUSTER_ASSIGN);
		snprintf(ctx -> output_data_file, DEFAULT_STR_LEN, "%s", OUT_DATA);

		while((opt = getopt_long(argc, argv, "i:t:d:o:a:h", long_options, NULL)) != -1)
		while((opt = getopt_long(argc, argv, "i:t:d:o:a:k:z:h", long_options, NULL)) != -1)
		{
		switch(opt)
		{
		@@ -103,10 +107,17 @@ void parse_args(global_context_t* ctx, int argc, char** argv)
		case 'a':
		strncpy(ctx -> output_assignment_file, optarg, DEFAULT_STR_LEN);
		break;
		case 'k':
		ctx -> k = atoi(optarg);
		break;
		case 'z':
		ctx -> z = atof(optarg);
		break;
		case 'h':
		mpi_printf(ctx, "%s\n", help);
		MPI_Finalize();
		exit(0);

		default:
		mpi_printf(ctx, "%s\n", help);
		MPI_Finalize();
		@@ -121,6 +132,42 @@ void parse_args(global_context_t* ctx, int argc, char** argv)
		if(err){MPI_Finalize(); exit(1);};
		}

		void print_hello(global_context_t* ctx)
		{
		char * hello = ""
		" _ _ \n"
		" __\| \| __ _ __\| \|_ __ \n"
		" / _` \|/ _` \|/ _` \| '_ \\ \n"
		"\| (_\| \| (_\| \| (_\| \| \|_) \| \n"
		" \\__,_\|\\__,_\|\\__,_\| .__/ \n"
		" \|_\| \n"
		" Distributed Advanced Density Peak";
		mpi_printf(ctx, "%s\n\n", hello);

		#if defined (_OPENMP)
		mpi_printf(ctx,"Running Hybrid (Openmp + MPI) code\n");
		#else
		mpi_printf(ctx,"Running pure MPI code\n");
		#endif

		#if defined (THREAD_FUNNELED)
		mpi_printf(ctx,"/!\\ Code built with MPI_THREAD_FUNNELED level\n");
		#else
		mpi_printf(ctx,"/!\\ Code built with MPI_THREAD_MULTIPLE level\n");
		#endif

		mpi_printf(ctx, "\nConfigs: \n");
		mpi_printf(ctx, "Data file .............> %s\n", ctx -> input_data_file);
		mpi_printf(ctx, "Input Type .............> float%d\n", ctx -> input_data_in_float32 ? 32 : 64);
		mpi_printf(ctx, "Dimensions .............> %d\n", ctx -> dims);
		mpi_printf(ctx, "Output data file .......> %s\n", ctx -> output_data_file);
		mpi_printf(ctx, "Output assignment file -> %s\n", ctx -> output_assignment_file);
		mpi_printf(ctx, "k ......................> %lu\n", ctx -> k);
		mpi_printf(ctx, "Z ......................> %.2lf\n", ctx -> z);
		mpi_printf(ctx, "\nRUNNING!\n");

		}

		int main(int argc, char** argv) {
		#if defined (_OPENMP)
		int mpi_provided_thread_level;
		@@ -160,32 +207,25 @@ int main(int argc, char** argv) {
		ctx.mpi_communicator = MPI_COMM_WORLD;
		get_context(&ctx);

		#if defined (_OPENMP)
		mpi_printf(&ctx,"Running Hybrid (Openmp + MPI) code\n");
		#else
		mpi_printf(&ctx,"Running pure MPI code\n");
		#endif

		#if defined (THREAD_FUNNELED)
		mpi_printf(&ctx,"/!\\ Code built with MPI_THREAD_FUNNELED level\n");
		#else
		mpi_printf(&ctx,"/!\\ Code built with MPI_THREAD_MULTIPLE level\n");
		#endif

		/*
		* Mock reading some files, one for each processor
		*/

		int d = 5;
		int k = 300;

		float_t* data;

		//parse command line

		parse_args(&ctx, argc, argv);
		print_hello(&ctx);
		/*
		* Generate a random matrix of lenght of some kind
		*/


		if(ctx.mpi_rank == 0)
		{
		simulate_master_read_and_scatter(5, 1000000, &ctx);
		@@ -211,10 +251,8 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
		TIME_DEF
		double elapsed_time;

		float_t z = 3;
		int halo = MY_TRUE;
		int halo = MY_FALSE;
		float_t tol = 0.002;
		int k = 300;

		if(I_AM_MASTER && ctx -> world_size <= 6)
		{
		@@ -343,7 +381,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
		dp_info[i].array_idx = -1;
		dp_info[i].kstar = -1;
		dp_info[i].is_center = -1;
		dp_info[i].cluster_idx = -1;
		dp_info[i].cluster_idx = -1717171717;
		//dp_info[i].halo_flag = 0;
		}
		ctx -> local_datapoints = dp_info;
		@@ -356,7 +394,8 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
		MPI_DB_PRINT("----- Performing ngbh search -----\n");
		MPI_Barrier(ctx -> mpi_communicator);

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

		MPI_Barrier(ctx -> mpi_communicator);
		elapsed_time = TIME_STOP;
		@@ -376,7 +415,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)


		compute_density_kstarnn_rma_v2(ctx, id, MY_FALSE);
		compute_correction(ctx, z);
		compute_correction(ctx, ctx -> z);
		elapsed_time = TIME_STOP;
		LOG_WRITE("Density estimate", elapsed_time)

		@@ -393,7 +432,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)


		TIME_START;
		Heuristic3(ctx, &clusters, z, halo);
		Heuristic3(ctx, &clusters, ctx -> z, halo);
		elapsed_time = TIME_STOP;
		LOG_WRITE("H3", elapsed_time)

		@@ -402,7 +441,7 @@ void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx)
		int* cl = (int)MY_MALLOC(ctx -> local_n_points sizeof(int));
		for(int i = 0; i < ctx -> local_n_points; ++i) cl[i] = ctx -> local_datapoints[i].cluster_idx;

		if(ctx -> world_size <= 6)
		if(ctx -> world_size <= 32)
		{
		big_ordered_buffer_to_file(ctx, cl, sizeof(int), ctx -> local_n_points, ctx -> output_assignment_file);
		big_ordered_buffer_to_file(ctx, ctx -> local_data, sizeof(double), ctx -> local_n_points * ctx -> dims, ctx -> output_data_file);

H1 optimization

Merge request reports