pddp_2means.c

	/*
	 * All data points belong to the root node at the beginning.
	 */
	#pragma vector aligned
	#pragma ivdep
	for (i = 0; i < data_points; i++) {
		indices[i] = i;
	}

	/*
	 * Initialize centroid.
	 */
	memset(centroid, 0, attributes * sizeof(double));

	#pragma omp parallel default(none) private(i, j, M_line) shared(M, indices, centroid, data_points, attributes, padded_attributes)
	{
	__attribute__((aligned(64))) double	centroid_local[attributes];

	memset(centroid_local, 0, attributes * sizeof(double));

	#pragma omp for
	for (i = 0; i < data_points; i++) {
		M_line = &M[indices[i] * padded_attributes];
		#pragma vector aligned
		#pragma ivdep
		for (j = 0; j < attributes; j++) {
			centroid_local[j] += M_line[j];
		}
	}

	#pragma vector aligned
	#pragma ivdep
	for (j = 0; j < attributes; j++) {
		#pragma omp atomic
		centroid[j] += centroid_local[j];
	}
	}

	#pragma vector aligned
	#pragma ivdep
	for (i = 0; i < attributes; i++) {
		centroid[i] /= data_points;
	}

	init_node(root, root, M, data_points, indices, centroid, attributes, padded_attributes, num_of_cores);

	current_level = 0;
	stop = 0;

#pragma omp parallel private(i) num_threads(num_of_cores < ((1 << max_level) / 4) ? num_of_cores : ((1 << max_level) / 4))
	{
#pragma omp single
	{
	do {
		cluster_num = 0;
#pragma omp taskgroup
		{
		error = posix_memalign((void **)(&cluster_nodes), 64, (1 << current_level) * sizeof(Node *));

		if (error != 0) {
			printf("ERROR: Could not allocate memory for vector of cluster nodes.\n");
			exit(0);
		}

		find_all_splittable_leaves(root, cluster_nodes, &cluster_num);

		for (i = 0; i < cluster_num; i++) {
			#pragma omp task
			process_node(cluster_nodes[i], M, attributes, padded_attributes, current_level, max_level);
		}
		}	// #pragma omp taskgroup

		free(cluster_nodes);

		/*
		 * If there have not been created as many clusters as requested,
		 * then we have to continue creating new clusters.
		 * 
		 * Except if there are no more clusters that can be split.
		 */
		if (leaves < clusters) {
			if (cluster_num == 0) {
				stop = 1;
				clusters = leaves;
				keep_all_nodes(root);
			}
		} else {
			/*
			 * The number of clusters created is equal or larger to the number of clusters requested.
			 * Check which ones have to be kept and if they are the ones required for the proper solution.
			 */
			stop = find_nodes_to_keep(root, clusters);
		}

		current_level = max_level;
		max_level++;
	} while (stop == 0);
	}	// #pragma omp single
	}	// #pragma omp parallel

	printf("Created a total of %ld clusters.\n", clusters);

	gettimeofday(&end, NULL);

	print_elapsed_time(start, end, "Time for calculations");

	gettimeofday(&start, NULL);

	assign_cluster_numbers(root, result, data_points);

	gettimeofday(&end, NULL);

	print_elapsed_time(start, end, "Time to assign cluster numbers");

	gettimeofday(&start, NULL);

	max_diam = sqrt(max_cluster_diameter(root, M, attributes, padded_attributes));

	gettimeofday(&end, NULL);

	print_elapsed_time(start, end, "Time to calculate maximum cluster diameter");

	printf("Maximum cluster diameter is %f\n", max_diam);

	gettimeofday(&start, NULL);

	min_dist = sqrt(min_cluster_distance(root, clusters, attributes));

	gettimeofday(&end, NULL);

	print_elapsed_time(start, end, "Time to calculate minimum cluster distance");

	printf("Minimum cluster distance is %f\n", min_dist);

	printf("Dunn index = %.40f\n", min_dist / max_diam);
}	// #pragma offload

	gettimeofday(&start, NULL);

        write_results(result, data_points, argv[2]);

        gettimeofday(&end, NULL);

        print_elapsed_time(start, end, "Time to write output data");

	return(0);
}