diff --git a/.gitignore b/.gitignore index a8635fedcfff415fe3e6a59a549cef8fe0c17a1b..a5298177b7336dd4da5b8c29c12383a2b3ed3c37 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,8 @@ bb scalability_results check.py var.py +test** +.agents/ +.codex/ +__pycache__/ +prompt.txt diff --git a/read_borders.py b/read_borders.py new file mode 100644 index 0000000000000000000000000000000000000000..f833b9d53527fd8395a72b1c8d54fee40a020686 --- /dev/null +++ b/read_borders.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +"""Read DADP borders files. + +The file written by write_borders_binary() is: + + repeated n_clusters times: + idx_t count + count * sparse_border_t + +where sparse_border_t is: + + idx_t i, j, idx; + float_t density, error; +""" + +from __future__ import annotations + +import argparse +import csv +import json +import struct +import sys +from pathlib import Path + + +def formats(idx_bytes: int, float_bytes: int) -> tuple[str, str]: + if idx_bytes == 8: + idx_fmt = "Q" + elif idx_bytes == 4: + idx_fmt = "I" + else: + raise ValueError("idx_bytes must be 4 or 8") + + if float_bytes == 8: + float_fmt = "d" + elif float_bytes == 4: + float_fmt = "f" + else: + raise ValueError("float_bytes must be 4 or 8") + + count_fmt = "<" + idx_fmt + + # Match the common C ABI layout for: + # idx_t i, j, idx; float_t density, error; + # The only padding case here is uint32_t indices followed by doubles. + pad = "4x" if idx_bytes == 4 and float_bytes == 8 else "" + border_fmt = "<" + idx_fmt * 3 + pad + float_fmt * 2 + return count_fmt, border_fmt + + +def read_exact(f, n: int, what: str) -> bytes: + data = f.read(n) + if len(data) != n: + raise EOFError(f"unexpected EOF while reading {what}: got {len(data)} of {n} bytes") + return data + + +def read_borders(path: Path, n_clusters: int | None, idx_bytes: int, float_bytes: int): + count_fmt, border_fmt = formats(idx_bytes, float_bytes) + count_size = struct.calcsize(count_fmt) + border_size = struct.calcsize(border_fmt) + + clusters = [] + with path.open("rb") as f: + c = 0 + while n_clusters is None or c < n_clusters: + count_buf = f.read(count_size) + if not count_buf: + break + if len(count_buf) != count_size: + raise EOFError(f"truncated count for cluster {c}") + + (count,) = struct.unpack(count_fmt, count_buf) + borders = [] + for _ in range(count): + i, j, idx, density, error = struct.unpack( + border_fmt, + read_exact(f, border_size, f"border record for cluster {c}"), + ) + borders.append( + { + "cluster_block": c, + "i": i, + "j": j, + "idx": idx, + "density": density, + "error": error, + } + ) + + clusters.append(borders) + c += 1 + + leftover = f.read(1) + if n_clusters is not None and leftover: + raise ValueError("file has extra bytes after the requested number of clusters") + + if n_clusters is not None and len(clusters) != n_clusters: + raise EOFError(f"expected {n_clusters} clusters, found {len(clusters)}") + + return clusters + + +def resolve_inputs(input_path: Path, n_clusters: int | None) -> tuple[Path, int | None]: + if input_path.is_dir(): + borders_path = input_path / "borders" + metadata_path = input_path / "metadata.json" + if n_clusters is None and metadata_path.exists(): + with metadata_path.open() as f: + n_clusters = int(json.load(f)["n_clusters"]) + else: + borders_path = input_path + + return borders_path, n_clusters + + +def write_csv(clusters, out): + writer = csv.DictWriter( + out, + fieldnames=["cluster_block", "i", "j", "idx", "density", "error"], + ) + writer.writeheader() + for borders in clusters: + writer.writerows(borders) + + +def main() -> int: + parser = argparse.ArgumentParser(description="Read DADP borders") + parser.add_argument("input", type=Path, help="output directory or path to borders") + parser.add_argument("--n-clusters", type=int, help="number of cluster blocks to read") + parser.add_argument("--idx-bytes", type=int, choices=(4, 8), default=8) + parser.add_argument("--float-bytes", type=int, choices=(4, 8), default=8) + parser.add_argument("--csv", type=Path, help="write all borders as CSV") + parser.add_argument("--json", type=Path, help="write nested cluster list as JSON") + parser.add_argument("--limit", type=int, default=20, help="rows to print when no output file is given") + args = parser.parse_args() + + borders_path, n_clusters = resolve_inputs(args.input, args.n_clusters) + clusters = read_borders(borders_path, n_clusters, args.idx_bytes, args.float_bytes) + + if args.csv: + with args.csv.open("w", newline="") as f: + write_csv(clusters, f) + + if args.json: + with args.json.open("w") as f: + json.dump(clusters, f, indent=2) + + if not args.csv and not args.json: + total = sum(len(borders) for borders in clusters) + print(f"clusters: {len(clusters)}") + print(f"borders: {total}") + printed = 0 + for borders in clusters: + for border in borders: + if printed >= args.limit: + return 0 + print(border) + printed += 1 + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/common/common.c b/src/common/common.c index 1106d6492377c3e2d1c2dc3c2d9d921b60b0792c..09d1655e1e12bb03a1fe42db68dbcbbc3017ae81 100644 --- a/src/common/common.c +++ b/src/common/common.c @@ -2,6 +2,7 @@ #include "mpi.h" #include #include +#include "../adp/adp.h" #define ARRAY_INCREMENT 100 @@ -26,6 +27,7 @@ void get_context(global_context_t* ctx) ctx -> dims = 0; ctx -> k = 300; ctx -> z = 3; + ctx -> parallel_dump = MY_FALSE; } void get_dataset_diagnostics(global_context_t* ctx, float_t* data) @@ -726,7 +728,7 @@ void test_file_path(const char* fname) void test_distributed_file_path(global_context_t* ctx, const char* fname) { - char out_path_w_proc_name[PATH_LEN]; + char out_path_w_proc_name[PATH_LEN]; snprintf(out_path_w_proc_name, PATH_LEN, "%s.%d", fname, ctx -> mpi_rank); FILE* file = fopen(out_path_w_proc_name,"w"); if(!file) @@ -741,3 +743,378 @@ void test_distributed_file_path(global_context_t* ctx, const char* fname) fclose(file); } + +static void* gather_ordered_buffer(global_context_t* ctx, const void* buffer, + size_t el_size, uint64_t n, + uint64_t* total_out) +{ + void* out = NULL; + uint64_t total = 0; + uint64_t local_n = n; + uint64_t* counts = NULL; + uint64_t* displs = NULL; + uint64_t* already_recv = NULL; + + MPI_Barrier(ctx->mpi_communicator); + MPI_Reduce(&local_n, &total, 1, MPI_UINT64_T, MPI_SUM, 0, ctx->mpi_communicator); + + if(I_AM_MASTER) + { + out = malloc(el_size * total); + counts = (uint64_t*)malloc(ctx->world_size * sizeof(uint64_t)); + displs = (uint64_t*)malloc(ctx->world_size * sizeof(uint64_t)); + already_recv = (uint64_t*)calloc(ctx->world_size, sizeof(uint64_t)); + CHECK_ALLOCATION_NO_CTX(out); + CHECK_ALLOCATION_NO_CTX(counts); + CHECK_ALLOCATION_NO_CTX(displs); + CHECK_ALLOCATION_NO_CTX(already_recv); + } + + MPI_Gather(&local_n, 1, MPI_UINT64_T, counts, 1, MPI_UINT64_T, 0, ctx->mpi_communicator); + + if(I_AM_MASTER) + { + displs[0] = 0; + for(int i = 0; i < ctx->world_size; ++i) + { + counts[i] *= el_size; + } + for(int i = 1; i < ctx->world_size; ++i) + { + displs[i] = displs[i - 1] + counts[i - 1]; + } + + if(counts[0] > 0) + { + memcpy(out, buffer, counts[0]); + } + + for(int r = 1; r < ctx->world_size; ++r) + { + while(already_recv[r] < counts[r]) + { + MPI_Status status; + int count_recv; + + MPI_Probe(r, r, ctx->mpi_communicator, &status); + MPI_Get_count(&status, MPI_BYTE, &count_recv); + MPI_Recv((char*)out + displs[r] + already_recv[r], + count_recv, MPI_BYTE, r, r, + ctx->mpi_communicator, MPI_STATUS_IGNORE); + already_recv[r] += (uint64_t)count_recv; + } + } + } + else + { + uint64_t already_sent = 0; + uint64_t bytes_to_send = n * el_size; + uint64_t default_msg_len = 100000000; + + while(already_sent < bytes_to_send) + { + int count_send = (int)MIN(default_msg_len, bytes_to_send - already_sent); + MPI_Send((char*)buffer + already_sent, count_send, MPI_BYTE, 0, + ctx->mpi_rank, ctx->mpi_communicator); + already_sent += (uint64_t)count_send; + } + } + + if(I_AM_MASTER) + { + free(counts); + free(displs); + free(already_recv); + } + + if(total_out) *total_out = total; + MPI_Barrier(ctx->mpi_communicator); + return out; +} + +static void remap_sparse_border(sparse_border_t* b, clusters_t* clusters, + idx_t* global_og_idxs, uint64_t n_points) +{ + b->i = clusters->centers.data[b->i]; + b->j = clusters->centers.data[b->j]; + if(b->idx < n_points) + { + b->idx = global_og_idxs[b->idx]; + } +} + +typedef struct assignment_pair_t +{ + idx_t idx; + int assignment; +} assignment_pair_t; + +static idx_t original_rank_start(global_context_t* ctx, int rank) +{ + idx_t base = ctx->n_points / ctx->world_size; + idx_t rem = ctx->n_points % ctx->world_size; + return rank * base + MIN((idx_t)rank, rem); +} + +static idx_t original_rank_count(global_context_t* ctx, int rank) +{ + idx_t base = ctx->n_points / ctx->world_size; + idx_t rem = ctx->n_points % ctx->world_size; + return base + ((idx_t)rank < rem); +} + +static int original_owner(global_context_t* ctx, idx_t original_idx) +{ + idx_t base = ctx->n_points / ctx->world_size; + idx_t rem = ctx->n_points % ctx->world_size; + idx_t larger_count = base + 1; + idx_t larger_total = rem * larger_count; + + if(original_idx < larger_total) + { + return original_idx / larger_count; + } + return rem + (original_idx - larger_total) / base; +} + +static void write_reordered_assignments_binary(global_context_t* ctx, + int* local_assignments, + const char* filepath) +{ + int* send_counts = (int*)calloc(ctx->world_size, sizeof(int)); + int* recv_counts = (int*)calloc(ctx->world_size, sizeof(int)); + int* send_displs = (int*)calloc(ctx->world_size, sizeof(int)); + int* recv_displs = (int*)calloc(ctx->world_size, sizeof(int)); + CHECK_ALLOCATION_NO_CTX(send_counts); + CHECK_ALLOCATION_NO_CTX(recv_counts); + CHECK_ALLOCATION_NO_CTX(send_displs); + CHECK_ALLOCATION_NO_CTX(recv_displs); + + for(idx_t i = 0; i < ctx->local_n_points; ++i) + { + int owner = original_owner(ctx, ctx->og_idxs[i]); + send_counts[owner]++; + } + + MPI_Alltoall(send_counts, 1, MPI_INT, recv_counts, 1, MPI_INT, + ctx->mpi_communicator); + + for(int i = 1; i < ctx->world_size; ++i) + { + send_displs[i] = send_displs[i - 1] + send_counts[i - 1]; + recv_displs[i] = recv_displs[i - 1] + recv_counts[i - 1]; + } + + int total_send = 0; + int total_recv = 0; + for(int i = 0; i < ctx->world_size; ++i) + { + total_send += send_counts[i]; + total_recv += recv_counts[i]; + } + + assignment_pair_t* send_pairs = (assignment_pair_t*)malloc(total_send * sizeof(assignment_pair_t)); + assignment_pair_t* recv_pairs = (assignment_pair_t*)malloc(total_recv * sizeof(assignment_pair_t)); + int* offsets = (int*)calloc(ctx->world_size, sizeof(int)); + CHECK_ALLOCATION_NO_CTX(send_pairs); + CHECK_ALLOCATION_NO_CTX(recv_pairs); + CHECK_ALLOCATION_NO_CTX(offsets); + + for(idx_t i = 0; i < ctx->local_n_points; ++i) + { + idx_t original_idx = ctx->og_idxs[i]; + int owner = original_owner(ctx, original_idx); + int pos = send_displs[owner] + offsets[owner]++; + send_pairs[pos].idx = original_idx; + send_pairs[pos].assignment = local_assignments[i]; + } + + int* send_counts_bytes = (int*)malloc(ctx->world_size * sizeof(int)); + int* recv_counts_bytes = (int*)malloc(ctx->world_size * sizeof(int)); + int* send_displs_bytes = (int*)malloc(ctx->world_size * sizeof(int)); + int* recv_displs_bytes = (int*)malloc(ctx->world_size * sizeof(int)); + CHECK_ALLOCATION_NO_CTX(send_counts_bytes); + CHECK_ALLOCATION_NO_CTX(recv_counts_bytes); + CHECK_ALLOCATION_NO_CTX(send_displs_bytes); + CHECK_ALLOCATION_NO_CTX(recv_displs_bytes); + + for(int i = 0; i < ctx->world_size; ++i) + { + send_counts_bytes[i] = send_counts[i] * sizeof(assignment_pair_t); + recv_counts_bytes[i] = recv_counts[i] * sizeof(assignment_pair_t); + send_displs_bytes[i] = send_displs[i] * sizeof(assignment_pair_t); + recv_displs_bytes[i] = recv_displs[i] * sizeof(assignment_pair_t); + } + + MPI_Alltoallv(send_pairs, send_counts_bytes, send_displs_bytes, MPI_BYTE, + recv_pairs, recv_counts_bytes, recv_displs_bytes, MPI_BYTE, + ctx->mpi_communicator); + + idx_t original_start = original_rank_start(ctx, ctx->mpi_rank); + idx_t original_count = original_rank_count(ctx, ctx->mpi_rank); + int* reordered = (int*)malloc(original_count * sizeof(int)); + CHECK_ALLOCATION_NO_CTX(reordered); + + for(idx_t i = 0; i < original_count; ++i) + { + reordered[i] = -1; + } + + for(int i = 0; i < total_recv; ++i) + { + idx_t pos = recv_pairs[i].idx - original_start; + if(pos < original_count) + { + reordered[pos] = recv_pairs[i].assignment; + } + } + + big_ordered_buffer_to_file(ctx, reordered, sizeof(int), original_count, filepath); + + free(send_counts); + free(recv_counts); + free(send_displs); + free(recv_displs); + free(send_pairs); + free(recv_pairs); + free(offsets); + free(send_counts_bytes); + free(recv_counts_bytes); + free(send_displs_bytes); + free(recv_displs_bytes); + free(reordered); +} + +void write_metadata_json(global_context_t* ctx, clusters_t* clusters, const char* output_dir) +{ + if(!(I_AM_MASTER)) return; + + char filepath[PATH_LEN]; + snprintf(filepath, PATH_LEN, "%s/metadata.json", output_dir); + + FILE* f = fopen(filepath, "w"); + if(!f) + { + printf("Cannot open file %s\n", filepath); + return; + } + + time_t now = time(NULL); + const char* tree_type_str = (ctx->local_tree_type == KD) ? "kd" : "vp"; + + fprintf(f, "{\n"); + fprintf(f, " \"n_points\": %lu,\n", ctx->n_points); + fprintf(f, " \"dims\": %u,\n", ctx->dims); + fprintf(f, " \"k\": %lu,\n", ctx->k); + fprintf(f, " \"z\": %f,\n", ctx->z); + fprintf(f, " \"n_clusters\": %lu,\n", clusters->centers.count); + fprintf(f, " \"timestamp\": %ld,\n", (long)now); + fprintf(f, " \"input_file\": \"%s\",\n", ctx->input_data_file); + fprintf(f, " \"n_mpi_procs\": %d,\n", ctx->world_size); + fprintf(f, " \"tree_type\": \"%s\"\n", tree_type_str); + fprintf(f, "}\n"); + + fclose(f); +} + +void write_centers_binary(global_context_t* ctx, clusters_t* clusters, const char* output_dir) +{ + if(!(I_AM_MASTER)) return; + + char filepath[PATH_LEN]; + snprintf(filepath, PATH_LEN, "%s/centers", output_dir); + + FILE* f = fopen(filepath, "wb"); + if(!f) + { + printf("Cannot open file %s\n", filepath); + return; + } + + fwrite(clusters->centers.data, sizeof(idx_t), clusters->centers.count, f); + fclose(f); +} + +void write_borders_binary(global_context_t* ctx, clusters_t* clusters, const char* output_dir) +{ + uint64_t total_og_idxs = 0; + idx_t* global_og_idxs = (idx_t*)gather_ordered_buffer(ctx, ctx->og_idxs, + sizeof(idx_t), + ctx->local_n_points, + &total_og_idxs); + + if(!(I_AM_MASTER)) return; + + if(total_og_idxs != ctx->n_points) + { + printf("Expected %lu original indices, gathered %lu\n", + (uint64_t)ctx->n_points, total_og_idxs); + } + + char filepath[PATH_LEN]; + snprintf(filepath, PATH_LEN, "%s/borders", output_dir); + + FILE* f = fopen(filepath, "wb"); + if(!f) + { + printf("Cannot open file %s\n", filepath); + free(global_og_idxs); + return; + } + + for(idx_t c = 0; c < clusters->centers.count; ++c) + { + idx_t count = clusters->sparse_borders[c].count; + fwrite(&count, sizeof(idx_t), 1, f); + + if(count > 0) + { + sparse_border_t* remapped = (sparse_border_t*)malloc(count * sizeof(sparse_border_t)); + for(idx_t i = 0; i < count; ++i) + { + remapped[i] = clusters->sparse_borders[c].data[i]; + remap_sparse_border(&remapped[i], clusters, global_og_idxs, total_og_idxs); + } + fwrite(remapped, sizeof(sparse_border_t), count, f); + free(remapped); + } + } + + fclose(f); + free(global_og_idxs); +} + +void write_assignments_binary(global_context_t* ctx, int* local_assignments, + const char* output_dir, int parallel_dump) +{ + char filepath[PATH_LEN]; + + if(parallel_dump) + { + snprintf(filepath, PATH_LEN, "%s/assignments", output_dir); + distributed_buffer_to_file(ctx, local_assignments, sizeof(int), + ctx->local_n_points, filepath); + } + else + { + snprintf(filepath, PATH_LEN, "%s/assignments", output_dir); + write_reordered_assignments_binary(ctx, local_assignments, filepath); + } +} + +void write_cluster_results(global_context_t* ctx, clusters_t* clusters, + int* local_assignments, const char* output_dir, + int parallel_dump) +{ + write_metadata_json(ctx, clusters, output_dir); + write_centers_binary(ctx, clusters, output_dir); + write_borders_binary(ctx, clusters, output_dir); + write_assignments_binary(ctx, local_assignments, output_dir, parallel_dump); + + if(parallel_dump) + { + char filepath[PATH_LEN]; + snprintf(filepath, PATH_LEN, "%s/og_idxs", output_dir); + distributed_buffer_to_file(ctx, ctx -> og_idxs, sizeof(idx_t), ctx -> local_n_points, filepath); + } +} diff --git a/src/common/common.h b/src/common/common.h index c4a882758346dffed188ae5c8067636a00829780..aa5dc0c54d3f9a19e0e0328080ca4722915ca2c0 100644 --- a/src/common/common.h +++ b/src/common/common.h @@ -194,6 +194,7 @@ typedef struct global_context_t char input_data_file[DEFAULT_STR_LEN]; char output_assignment_file[DEFAULT_STR_LEN]; char output_data_file[DEFAULT_STR_LEN]; + int parallel_dump; } global_context_t; typedef struct pointset_t @@ -243,3 +244,8 @@ void get_dataset_diagnostics(global_context_t* ctx, float_t* data); void test_distributed_file_path(global_context_t* ctx, const char* fname); void distributed_buffer_to_file(global_context_t* ctx, void* buffer, size_t el_size, uint64_t n, const char* fname); + +typedef struct clusters_t clusters_t; +void write_cluster_results(global_context_t* ctx, clusters_t* clusters, + int* local_assignments, const char* output_dir, + int parallel_dump); diff --git a/src/main/main.c b/src/main/main.c index deb2d209305cd2d50d9bece3b6c4c2e7e08c450e..48f54eaeb82cf6636ddae15aff7e9feaec759387 100644 --- a/src/main/main.c +++ b/src/main/main.c @@ -24,6 +24,7 @@ struct option long_options[] = {"kngbh", optional_argument, 0, 'k'}, {"z", optional_argument, 0, 'z'}, {"knn-strategy", optional_argument, 0, 's'}, + {"parallel-out", no_argument, 0, 'p'}, {"help", optional_argument, 0, 'h'}, {0, 0, 0, 0} }; @@ -39,7 +40,8 @@ const char* help = "Distributed Advanced Density Peak\n"\ " of the datapoints. Produces files like `data.[rank]`, `og_idx.[rank]` and `assignment.[rank]`\n"\ " -k --kngbh (optional) number of nearest neighbors to compute\n"\ " -s --knn-strategy (optional) strategy for local knn (allowed values 'kd' for kdtree, 'vp' for ball-tree (vantage point)\n"\ - " -z --z (optional) number of nearest neighbors to compute\n"; + " -z --z (optional) z parameter for density correction\n" + " -p --parallel-out (optional) parallel dump of assignments\n"; void parse_args(global_context_t* ctx, int argc, char** argv) { @@ -48,7 +50,7 @@ void parse_args(global_context_t* ctx, int argc, char** argv) int input_file_set = 0; int input_type_set = 0; - while((opt = getopt_long(argc, argv, "i:t:d:o:k:z:hs:", long_options, NULL)) != -1) + while((opt = getopt_long(argc, argv, "i:t:d:o:k:z:hs:p", long_options, NULL)) != -1) { switch(opt) { @@ -105,6 +107,10 @@ void parse_args(global_context_t* ctx, int argc, char** argv) } break; + case 'p': + ctx -> parallel_dump = MY_TRUE; + break; + default: mpi_printf(ctx, "%s\n", help); MPI_Finalize(); @@ -243,6 +249,26 @@ void mpiio_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) int halo = MY_TRUE; float_t tol = 0.002; + if(I_AM_MASTER) + { + char out_file[200]; + + snprintf(out_file, 200, "%s/metadata.json", ctx->output_data_file); + test_file_path(out_file); + + snprintf(out_file, 200, "%s/centers", ctx->output_data_file); + test_file_path(out_file); + + snprintf(out_file, 200, "%s/borders", ctx->output_data_file); + test_file_path(out_file); + + if(!(ctx->parallel_dump)) + { + snprintf(out_file, 200, "%s/assignments", ctx->output_data_file); + test_file_path(out_file); + } + } + if(ctx -> world_size <= 32) { if(I_AM_MASTER) @@ -252,9 +278,6 @@ void mpiio_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); test_file_path(out_file); - snprintf(out_file, 200, "%s/assignment", ctx->output_data_file); - test_file_path(out_file); - snprintf(out_file, 200, "%s/data", ctx->output_data_file); test_file_path(out_file); } @@ -265,13 +288,16 @@ void mpiio_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); test_distributed_file_path(ctx, out_file); - snprintf(out_file, 200, "%s/assignment", ctx->output_data_file); + snprintf(out_file, 200, "%s/data", ctx->output_data_file); test_distributed_file_path(ctx, out_file); + } - snprintf(out_file, 200, "%s/data", ctx->output_data_file); + if(ctx->parallel_dump) + { + char out_file[200]; + snprintf(out_file, 200, "%s/assignments", ctx->output_data_file); test_distributed_file_path(ctx, out_file); } - TIME_START; ctx->local_data = mpiio_read_data_file(ctx, ctx->input_data_file, ctx->dims, @@ -408,7 +434,7 @@ void mpiio_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) TIME_START; int* cl = (int*)MY_MALLOC(2 * ctx -> local_n_points * sizeof(int)); - float_t* data_to_write = (float_t*)MY_MALLOC(ctx -> local_n_points * ctx -> k * sizeof(float_t)); + float_t* data_to_write = (float_t*)MY_MALLOC(ctx -> local_n_points * ctx -> dims * sizeof(float_t)); switch(ctx -> local_tree_type) { @@ -442,34 +468,28 @@ void mpiio_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) break; } - if(ctx -> world_size <= 32) - { - char out_file[200]; - - snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); - big_ordered_buffer_to_file(ctx, ctx->og_idxs, sizeof(idx_t), ctx -> local_n_points, out_file); - - snprintf(out_file, 200, "%s/assignment", ctx->output_data_file); - big_ordered_buffer_to_file(ctx, cl, sizeof(int), ctx -> local_n_points, out_file); - - snprintf(out_file, 200, "%s/data", ctx->output_data_file); - big_ordered_buffer_to_file(ctx, data_to_write, sizeof(double), ctx -> local_n_points * ctx -> dims, out_file); - } - else - { - char out_file[200]; - - snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); - distributed_buffer_to_file(ctx, ctx->og_idxs, sizeof(idx_t), ctx -> local_n_points, out_file); - - snprintf(out_file, 200, "%s/assignment", ctx->output_data_file); - distributed_buffer_to_file(ctx, cl, sizeof(int), ctx -> local_n_points, out_file); - - snprintf(out_file, 200, "%s/data", ctx->output_data_file); - distributed_buffer_to_file(ctx, data_to_write, sizeof(double), ctx -> local_n_points * ctx -> dims, out_file); - } - - + // if(ctx -> world_size <= 32) + // { + // char out_file[200]; + // + // snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); + // big_ordered_buffer_to_file(ctx, ctx->og_idxs, sizeof(idx_t), ctx -> local_n_points, out_file); + // + // snprintf(out_file, 200, "%s/data", ctx->output_data_file); + // big_ordered_buffer_to_file(ctx, data_to_write, sizeof(float_t), ctx -> local_n_points * ctx -> dims, out_file); + // } + // else + // { + // char out_file[200]; + // + // snprintf(out_file, 200, "%s/og_idx", ctx->output_data_file); + // distributed_buffer_to_file(ctx, ctx->og_idxs, sizeof(idx_t), ctx -> local_n_points, out_file); + // + // snprintf(out_file, 200, "%s/data", ctx->output_data_file); + // distributed_buffer_to_file(ctx, data_to_write, sizeof(float_t), ctx -> local_n_points * ctx -> dims, out_file); + // } + + write_cluster_results(ctx, &clusters, cl, ctx->output_data_file, ctx->parallel_dump); elapsed_time = TIME_STOP; LOG_WRITE("Write results to file", elapsed_time); diff --git a/src/tree/tree.c b/src/tree/tree.c index 3f8b3e4ad95097795fb71b503f2a075b7ec8a5fc..8d0579fd67d70e609e61f79fcd618d85ac4ca6d0 100644 --- a/src/tree/tree.c +++ b/src/tree/tree.c @@ -290,40 +290,15 @@ void compute_bounding_box_data(global_context_t *ctx, float_t *data, float_t* lb /* compute minimum and maximum for each dimensions, store them in local bb */ /* each processor on its own */ - // this moves memory maybe there is a better way - #pragma omp parallel + // Use OpenMP 4.5+ array reduction with runtime size + #pragma omp parallel for reduction(min:lb[:ctx->dims]) reduction(max:ub[:ctx->dims]) + for (size_t i = 0; i < n_points; ++i) { - float_t* pvt_lb = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t)); - float_t* pvt_ub = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t)); - for (size_t d = 0; d < ctx->dims; ++d) + for (size_t d = 0; d < ctx->dims; ++d) { - pvt_lb[d] = FLT_MAX; - pvt_ub[d] = -FLT_MAX; + lb[d] = MIN(data[i*ctx->dims + d], lb[d]); + ub[d] = MAX(data[i*ctx->dims + d], ub[d]); } - - #pragma omp for - for (size_t i = 0; i < n_points; ++i) - { - for (size_t d = 0; d < ctx->dims; ++d) - { - pvt_lb[d] = MIN(data[i*ctx->dims + d], pvt_lb[d]); - pvt_ub[d] = MAX(data[i*ctx->dims + d], pvt_ub[d]); - } - } - - #pragma omp critical (bounding_box_reduction) - { - for (size_t d = 0; d < ctx->dims; ++d) - { - lb[d] = MIN(pvt_lb[d], lb[d]); - ub[d] = MAX(pvt_ub[d], ub[d]); - } - } - - - free(pvt_lb); - free(pvt_ub); - } /*get the bounding box */ @@ -342,54 +317,28 @@ void compute_bounding_box_data(global_context_t *ctx, float_t *data, float_t* lb } void compute_bounding_box_pointset(global_context_t *ctx, pointset_t *ps) { - #define local_data ps->datapoints - #define lb ps->lb_box - #define ub ps->ub_box - - for (size_t d = 0; d < ps->dims; ++d) - { - ps->lb_box[d] = FLT_MAX; - ps->ub_box[d] = -FLT_MAX; + // Use local pointers instead of macros to enable efficient OpenMP reduction + point_t* local_data = ps->datapoints; + float_t* lb = ps->lb_box; + float_t* ub = ps->ub_box; + + // Initialize bounding boxes + for (size_t d = 0; d < ps->dims; ++d) { + lb[d] = FLT_MAX; + ub[d] = -FLT_MAX; } - /* compute minimum and maximum for each dimensions, store them in local bb */ /* each processor on its own */ - // this moves memory maybe there is a better way - #pragma omp parallel + #pragma omp parallel for reduction(min:lb[:ps->dims]) reduction(max:ub[:ps->dims]) + for (size_t i = 0; i < ps->n_points; ++i) { - float_t* pvt_lb = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t)); - float_t* pvt_ub = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t)); - for (size_t d = 0; d < ps->dims; ++d) - { - pvt_lb[d] = FLT_MAX; - pvt_ub[d] = -FLT_MAX; - } - - #pragma omp for - for (size_t i = 0; i < ps->n_points; ++i) - { - for (size_t d = 0; d < ps->dims; ++d) - { - pvt_lb[d] = MIN(local_data[i].data[d], pvt_lb[d]); - pvt_ub[d] = MAX(local_data[i].data[d], pvt_ub[d]); - } - } - - #pragma omp critical (bounding_box_reduction) + for (size_t d = 0; d < ps->dims; ++d) { - for (size_t d = 0; d < ps->dims; ++d) - { - lb[d] = MIN(pvt_lb[d], lb[d]); - ub[d] = MAX(pvt_ub[d], ub[d]); - } + lb[d] = MIN(local_data[i].data[d], lb[d]); + ub[d] = MAX(local_data[i].data[d], ub[d]); } - - - free(pvt_lb); - free(pvt_ub); - } /*get the bounding box */