cmake_minimum_required(VERSION 3.20)
project(FilterProject LANGUAGES C CXX)
add_executable(program main.cpp filter.F90)
set_source_files_properties(filter.F90 PROPERTIES LANGUAGE C)
target_compile_options(program PRIVATE
    -O3                     # 启用最高级别优化
    -march=native            # 针对本地硬件生成优化代码
    -ftree-vectorize         # 启用自动向量化
    -funroll-loops
)

#include <vector>
#include <utility>
#include <cstdlib>
#include <algorithm>
#include <omp.h>
typedef int element_t;
size_t lcs(element_t* arr_1, element_t* arr_2, size_t len_1, size_t len_2) {
    std::vector<size_t> previous(len_2 + 1, 0);  // previous row
    std::vector<size_t> current(len_2 + 1, 0);   // current row
    // 对内存访问顺序进行优化，保证数据按列访问
    for (size_t i = 1; i <= len_1; ++i) {
        #pragma omp parallel for num_threads(16)
        for (size_t j = 1; j <= len_2; ++j) {
            if (arr_1[i - 1] == arr_2[j - 1]) {
                current[j] = previous[j - 1] + 1;
            } else {
                current[j] = std::max(current[j - 1], previous[j]);
            }
        }
        previous.swap(current); // 交换 current 和 previous
    }
    return previous[len_2];  // 结果在最后的 previous 行
}

#include <utility>
#include <cstdlib>
#include <algorithm>
#include <omp.h>
#include <arm_neon.h>
typedef int element_t;
size_t lcs(element_t* arr_1, element_t* arr_2, size_t len_1, size_t len_2) {
    size_t* previous = (size_t*)calloc(len_2 + 1, sizeof(size_t));
    size_t* current = (size_t*)calloc(len_2 + 1, sizeof(size_t));

    for (size_t i = 1; i <= len_1; ++i) {
        // Parallelize the inner loop for each row, ensure correct data handling
        #pragma omp parallel for
        for (size_t j = 1; j <= len_2; j += 4) {  // Process 4 elements at a time with NEON
            // Load 4 elements from arr_2[j-1..j+3]
            int32x4_t arr_2_values = vld1q_s32(&arr_2[j - 1]);

            // Process 4 elements of arr_1[i-1]
            size_t prev_val = previous[j - 1]; // Initialize with previous[j-1]
            int32x4_t current_vals = vld1q_s32(&arr_1[i - 1]);

            // Compare 4 elements at once
            uint32x4_t match_mask = vceqq_s32(arr_2_values, current_vals);

            // Find maximum in current and previous rows, update current
            for (size_t k = 0; k < 4 && j + k <= len_2; ++k) {
                if (arr_1[i - 1] == arr_2[j + k - 1]) {
                    current[j + k] = previous[j + k - 1] + 1;
                } else {
                    current[j + k] = std::max(current[j + k - 1], previous[j + k]);
                }
            }    
        }
       
        // Ensure all threads are synchronized before swapping
        #pragma omp barrier

        // Swap current and previous arrays after all threads finish updating
        std::swap(current, previous);
    }
    size_t result = previous[len_2];
    free(previous);
    free(current);

    return result;
}

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <vector_types.h>
#include <cuda_runtime.h>

typedef double d_t;
struct d3_t { d_t x, y, z; };

__device__ d3_t operator-(const d3_t& a, const d3_t& b) {
    return {a.x - b.x, a.y - b.y, a.z - b.z};
}

__device__ d_t norm(const d3_t& v) {
    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
}

__constant__ d3_t c_src;

__global__ void computeDist1(const d3_t* __restrict__ mir, d_t* dist1, int64_t mirn) {
    int64_t j = blockIdx.x * blockDim.x + threadIdx.x;
    if (j >= mirn) return;
    dist1[j] = norm(mir[j] - c_src);
}

__global__ void computeData(const d3_t* mir, const d_t* dist1,
                            const d3_t* sen, int64_t mirn, int64_t senn,
                            d_t* data, int tile_size) {
    extern __shared__ char shared_mem[];
    d3_t* shared_mir = (d3_t*)shared_mem;
    d_t* shared_dist1 = (d_t*)(shared_mem + tile_size * sizeof(d3_t));

    int64_t sen_idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (sen_idx >= senn) return;

    const d3_t sensor = sen[sen_idx];
    const d_t freq_factor = 6.283185307179586 * 2000;
    d_t a = 0.0, b = 0.0;

    // 分块处理所有mir元素
    for (int64_t tile_start = 0; tile_start < mirn; tile_start += tile_size) {
        int64_t tile_end = min(tile_start + tile_size, mirn);
        int64_t num_in_tile = tile_end - tile_start;

        // 协作加载当前tile到共享内存
        for (int64_t j = threadIdx.x; j < num_in_tile; j += blockDim.x) {
            int64_t global_idx = tile_start + j;
            if (global_idx < mirn) {
                shared_mir[j] = mir[global_idx];
                shared_dist1[j] = dist1[global_idx];
            }
        }
        __syncthreads();

        // 处理当前tile中的元素
        for (int64_t j = 0; j < num_in_tile; ++j) {
            d3_t mirror = shared_mir[j];
            d_t dist1_val = shared_dist1[j];
            d_t dist2 = norm(mirror - sensor);
            d_t phase = freq_factor * (dist1_val + dist2);

            d_t cos_val, sin_val;
            sincos(phase, &sin_val, &cos_val);
            a += cos_val;
            b += sin_val;
        }
        __syncthreads(); // 确保下一tile加载前完成计算
    }

    data[sen_idx] = hypot(a, b);
}


int main() {
    FILE* fi = fopen("in.data", "rb");
    d3_t src;
    int64_t mirn, senn;
    d3_t *mir, *sen;

    fread(&src, sizeof(d3_t), 1, fi);
    fread(&mirn, sizeof(int64_t), 1, fi);
    mir = (d3_t*)malloc(mirn * sizeof(d3_t));
    fread(mir, sizeof(d3_t), mirn, fi);

    fread(&senn, sizeof(int64_t), 1, fi);
    sen = (d3_t*)malloc(senn * sizeof(d3_t));
    fread(sen, sizeof(d3_t), senn, fi);
    fclose(fi);

    // 设备内存分配
    d3_t *d_mir, *d_sen;
    d_t *d_data, *d_dist1;
    cudaMalloc(&d_mir, mirn * sizeof(d3_t));
    cudaMalloc(&d_sen, senn * sizeof(d3_t));
    cudaMalloc(&d_data, senn * sizeof(d_t));
    cudaMalloc(&d_dist1, mirn * sizeof(d_t)); // 预计算dist1的数组

    // 将src拷贝到常量内存
    cudaMemcpyToSymbol(c_src, &src, sizeof(d3_t));

    // 创建多个CUDA流
    const int num_streams = 4;
    cudaStream_t streams[num_streams];
    for (int i = 0; i < num_streams; ++i) {
        cudaStreamCreate(&streams[i]);
    }

    // 同步传输mir数据并预计算dist1
    cudaMemcpy(d_mir, mir, mirn * sizeof(d3_t), cudaMemcpyHostToDevice);
    dim3 block_pre(512);
    dim3 grid_pre((mirn + block_pre.x - 1) / block_pre.x);
    computeDist1<<<grid_pre, block_pre>>>(d_mir, d_dist1, mirn);
    cudaDeviceSynchronize(); // 确保预计算完成

    // 检查共享内存需求
    int max_shared_mem_per_block;
    cudaDeviceGetAttribute(&max_shared_mem_per_block, cudaDevAttrMaxSharedMemoryPerBlock, 0);
    size_t per_element_size = sizeof(d3_t) + sizeof(d_t); // 每个mir元素需要的共享内存
    size_t tile_size = max_shared_mem_per_block / per_element_size;
    if (tile_size == 0) tile_size = 1; // 确保至少处理一个元素

    // 分块处理时使用二维拷贝优化
    size_t chunk = (senn + 3) / 4;
    for (int i = 0; i < num_streams; ++i) {
        cudaStreamCreate(&streams[i]);
        size_t offset = i * chunk;
        size_t valid = min(chunk, senn - offset);
        if (valid == 0) continue;

        // 异步拷贝传感器数据
        cudaMemcpyAsync(d_sen + offset, sen + offset,
                        valid * sizeof(d3_t),
                        cudaMemcpyHostToDevice,
                        streams[i]);

        dim3 block(512);
        dim3 grid((valid + block.x - 1) / block.x);
       
        // 计算所需共享内存大小
        size_t shared_mem_size = tile_size * (sizeof(d3_t) + sizeof(d_t));
       
        computeData<<<grid, block, shared_mem_size, streams[i]>>>(
            d_mir, d_dist1, d_sen + offset,
            mirn, valid, d_data + offset, tile_size
        );
    }

    // 同步所有流并销毁
    for (int i = 0; i < num_streams; ++i) {
        cudaStreamSynchronize(streams[i]);
        cudaStreamDestroy(streams[i]);
    }

    // 拷贝结果回主机
    d_t* host_data = (d_t*)malloc(senn * sizeof(d_t));
    cudaMemcpy(host_data, d_data, senn * sizeof(d_t), cudaMemcpyDeviceToHost);

    // 输出结果
    FILE* fo = fopen("out.data", "wb");
    fwrite(host_data, sizeof(d_t), senn, fo);
    fclose(fo);

    // 释放资源
    free(mir);
    free(sen);
    free(host_data);
    cudaFree(d_mir);
    cudaFree(d_sen);
    cudaFree(d_data);

    return 0;
}

g++ -O3 -march=native vanity.cpp -o vanity -lsecp256k1 -lcrypto -fopenmp -lpthread
./vanity

#include <stdio.h>
#include<omp.h>
#include <arm_sve.h>
#include <stddef.h>
#include "hpl-ai.h"

#define A(i, j) *HPLAI_INDEX2D(A, (i), (j), lda)

void sgetrf_nopiv(int m, int n, float *A, int lda) {
    int j;
    int nb = 32;
    int jb = nb;

    // Use unblock code.
    if (nb > m || nb > n) {
        sgetrf2_nopiv(m, n, A, lda);
        return;
    }

    int min_mn = m < n ? m : n;

    for (j = 0; j < min_mn; j += nb) {
        if (min_mn - j < nb) {
            jb = min_mn - j;
        }

        // Factor panel
        sgetrf2_nopiv(m - j, jb, &A[j * lda + j], lda);

        if (j + jb < n) {
            // strsm optimization with SVE
            svbool_t pg = svptrue_b32(); // Create a mask for all elements
            svfloat32_t vec_A_jj = svld1_f32(pg, &A[j * lda + j]);
            svfloat32_t vec_A_jj_next = svld1_f32(pg, &A[(j + jb) * lda + j]);

            // Create a vector with all elements set to 1.0f
            svfloat32_t one = svdup_f32(1.0f);
            
            // Perform matrix multiply add (svmla)
            svfloat32_t result = svmla_f32_x(pg, vec_A_jj, vec_A_jj_next, one);
            svst1_f32(pg, &A[(j + jb) * lda + j], result);

            // GEMM with SVE
            if (j + jb < m) {
                svfloat32_t vec_A_next_row = svld1_f32(pg, &A[(j + jb) * lda + j]);
                svfloat32_t vec_A_next_col = svld1_f32(pg, &A[j * lda + j + jb]);
                svfloat32_t updated_result = svmul_f32_x(pg, vec_A_next_row, vec_A_next_col);
                svst1_f32(pg, &A[(j + jb) * lda + j + jb], updated_result);
            }
        }
    }
}

import Pkg
Pkg.activate(temp=true)
Pkg.add("包")

@inbounds function topk(data::AbstractVector{T}, k) where T
    # your code here
end

using Base.Threads

# 原子类型定义，用于线程安全的操作
mutable struct Atomic{T}; @atomic x::T; end

@inbounds function topk(data::AbstractVector{T}, k) where T
    # 将数值和下标配对
    indexed_data = [(i, data[i]) for i in 1:length(data)]

    # 数据划分为8份
    num_threads = nthreads()
    chunk_size = div(length(data), num_threads)
    chunks = [indexed_data[(i-1)*chunk_size+1:min(i*chunk_size, length(data))] for i in 1:num_threads]

    # 定义一个并行函数，处理每个线程的排序和选择最大K个
    function parallel_sort(chunk)
        return sort(chunk, by = x -> (-x[2], x[1]))[1:k]
    end

    # 创建一个数组来存储每个线程的局部结果
    results = Vector{Tuple{Int, T}}(undef, k * num_threads)

    # 使用多个线程并行化排序操作
    @threads for i in 1:num_threads
        local_sorted = parallel_sort(chunks[i])
        # 将结果写入共享数组
        for j in 1:k
            results[(i-1)*k + j] = local_sorted[j]
        end
    end

    # 合并所有线程的结果并进行排序
    final_sorted = sort(results, by = x -> (-x[2], x[1]))

    # 返回前K个下标
    return [final_sorted[i][1] for i in 1:k]
end

import Pkg
Pkg.activate(temp=true)
Pkg.add("DataStructures")
using Base.Threads
using DataStructures
 
# 并行化选择TopK的下标的核心计算函数
function parallel_topk_indices(chunk::Vector{Tuple{Int, T}}, k::Int) where T
    # 使用最大堆来选择TopK的下标
    max_heap = BinaryMaxHeap{Tuple{T, Int}}()
   
   for (idx, value) in chunk
        # 插入堆中，堆按值的大小自动排序（这里是负值因为BinaryMaxHeap是最大堆）
        push!(max_heap, (-value, idx))  # 负值以模拟最大堆行为
        # 如果堆的大小超过k，移除最小的元素
        if length(max_heap) > k
            pop!(max_heap)
        end
    end
 
    # 从堆中提取元素，并获取下标
    result = Vector{Int}(undef, length(max_heap))
    @simd for i in 1:length(max_heap)
        (_, idx) = pop!(max_heap)
        result[i] = idx
    end
    return result
end
 
# 主函数：计算前K大的数的下标
@inbounds function topk(data::Vector{T}, k::Int) where T
    # 数据划分为线程数份
    num_threads = nthreads()
    chunk_size = div(length(data), num_threads)
 
    # 预分配存储结果的数组
    results = Vector{Int}(undef, k * num_threads)
    all_indices = Vector{Int}(undef, k * num_threads)  # 确保大小足够大
   
    # 使用多个线程并行化排序操作
    @threads for i in 1:num_threads
        # 计算每个线程的分片范围
        start_idx = (i - 1) * chunk_size + 1
        end_idx = min(i * chunk_size, length(data))
 
        # 提取当前分片的数据
        local_chunk = [(idx, data[idx]) for idx in start_idx:end_idx]
 
        # 对当前分片进行TopK下标选择
        local_topk_indices = parallel_topk_indices(local_chunk, k)
 
        # 将结果写入共享数组（按线程对应位置）
        offset = (i - 1) * k
        for j in 1:k
            results[offset + j] = local_topk_indices[j]
        end
    end
 
    # 合并所有线程的结果
    # 只合并有效的K个元素，避免超出数组边界
    all_indices[1:k*num_threads] .= results[1:k*num_threads]
 
    # 使用这些下标从原始数据中获取对应的值
    values = data[all_indices]
 
    # 对值进行排序，并选择前K个元素对应的下标
    sort_perm = sortperm(values, rev=true)
    final_topk_indices = all_indices[sort_perm[1:k]]
 
    return final_topk_indices
end

import Pkg
Pkg.activate(temp=true)
Pkg.add("DataStructures")
using Base.Threads
using DataStructures

# 并行化选择TopK下标的核心计算函数
function parallel_topk_indices(data::Vector{T}, start_idx::Int, end_idx::Int, k::Int) where T
    # 处理无效的k值或空区间
    k = max(0, min(k, end_idx - start_idx + 1))
    k == 0 && return Int[]
   
    # 使用最小堆维护当前最大的k个元素
    min_heap = BinaryMinHeap{Tuple{T, Int}}()
   
    @inbounds for idx in start_idx:end_idx
        value = data[idx]
        if length(min_heap) < k
            push!(min_heap, (value, idx))
        else
            if value > first(min_heap)[1]
                pop!(min_heap)
                push!(min_heap, (value, idx))
            end
        end
    end
   
    # 将堆中元素按降序输出
    n = length(min_heap)
    result = Vector{Int}(undef, n)
    for i in n:-1:1
        result[i] = pop!(min_heap)[2]
    end
    return result
end

# 主函数：计算前K大的数的下标
@inbounds function topk(data::Vector{T}, k::Int) where T
    len = length(data)
    k = min(k, len)
    len == 0 && return Int[]
   
    num_threads = nthreads()
    chunk_size = div(len, num_threads)
   
    # 存储各线程结果的数组
    thread_results = Vector{Vector{Int}}(undef, num_threads)
   
    @threads for i in 1:num_threads
        start_idx = (i-1)*chunk_size + 1
        end_idx = i == num_threads ? len : i * chunk_size
       
        # 直接处理数据分片，避免创建临时数组
        thread_results[i] = parallel_topk_indices(data, start_idx, end_idx, k)
    end
   
    # 合并所有结果
    all_indices = vcat(thread_results...)
   
    # 当总候选数不超过k时直接返回
    length(all_indices) <= k && return all_indices
   
    # 高效部分排序获取前k大元素的下标
    values = data[all_indices]
    sort_order = partialsortperm(values, 1:k, rev=true)
    return all_indices[sort_order]
end

 for (int i = 0; i < n; i++) {    
   for (int j = 0; j < n; j++) {
        fo << grid[i * n + j] << " ";    
       }    
      fo << std::endl;
}

. /etc/profile >/dev/null #加载环境krun init >/dev/nullif [ "$JOB_COMPLETION_INDEX" -eq 0 ]; then
    krun host -n -H > hostfile
    mpirun --hostfile hostfile -np 64 -ppn 16 /data/forest /problem/input.dat /data/output.datelse
    sleep infinityfi

. /etc/profile >/dev/null #加载环境krun init >/dev/null
krun host -n -H > hostfile
mpirun --hostfile hostfile -np 4 -ppn 1 /data/forest /problem/input.dat /data/output.dat

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <unordered_map>

struct Event {
    int step;
    int type;
    int params[4];
};

int main(int argc, char** argv) {
    MPI_Init(&argc, &argv);

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    int n, m, t;
    std::unordered_map<int, Event> events;

    int* global_grid = nullptr;
    int* event_data = nullptr;

    if (rank == 0) {
        FILE* fin = fopen(argv[1], "r");
        fscanf(fin, "%d %d %d", &n, &m, &t);

        global_grid = new int[n * n];
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                fscanf(fin, "%d", &global_grid[i * n + j]);
            }
        }

        events.reserve(m);
        for (int k = 0; k < m; ++k) {
            Event e;
            fscanf(fin, "%d %d", &e.step, &e.type);
            if (e.type == 1) {
                fscanf(fin, "%d %d", &e.params[0], &e.params[1]);
                e.params[2] = e.params[3] = 0;
            } else {
                fscanf(fin, "%d %d %d %d", &e.params[0], &e.params[1], &e.params[2], &e.params[3]);
            }
            events[e.step] = e;
        }
        fclose(fin);

        event_data = new int[m * 6];
        int idx = 0;
        for (const auto& pair : events) {
            const Event& e = pair.second;
            event_data[idx++] = e.step;
            event_data[idx++] = e.type;
            event_data[idx++] = e.params[0];
            event_data[idx++] = e.params[1];
            event_data[idx++] = e.params[2];
            event_data[idx++] = e.params[3];
        }
    }

    MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&m, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&t, 1, MPI_INT, 0, MPI_COMM_WORLD);

    if (rank != 0) {
        event_data = new int[m * 6];
    }
    MPI_Bcast(event_data, m * 6, MPI_INT, 0, MPI_COMM_WORLD);

    if (rank != 0) {
        for (int i = 0; i < m; ++i) {
            Event e;
            e.step = event_data[i * 6];
            e.type = event_data[i * 6 + 1];
            e.params[0] = event_data[i * 6 + 2];
            e.params[1] = event_data[i * 6 + 3];
            e.params[2] = event_data[i * 6 + 4];
            e.params[3] = event_data[i * 6 + 5];
            events[e.step] = e;
        }
        delete[] event_data;
    } else {
        delete[] event_data;
    }

    int remainder = n % size;
    int local_rows = n / size + (rank < remainder ? 1 : 0);
    int start_row = (rank < remainder) ? (rank * (n / size + 1)) : (remainder * (n / size + 1) + (rank - remainder) * (n / size));
    int end_row = start_row + local_rows - 1;

    int* local_data = new int[local_rows * n];
    if (rank == 0) {
        int* ptr = global_grid;
        for (int dest = 0; dest < size; ++dest) {
            int dest_remainder = n % size;
            int dest_local_rows = n / size + (dest < dest_remainder ? 1 : 0);
            int dest_start_row = (dest < dest_remainder) ? (dest * (n / size + 1)) : (dest_remainder * (n / size + 1) + (dest - dest_remainder) * (n / size));

            if (dest == 0) {
                memcpy(local_data, ptr, dest_local_rows * n * sizeof(int));
            } else {
                MPI_Send(ptr, dest_local_rows * n, MPI_INT, dest, 0, MPI_COMM_WORLD);
            }
            ptr += dest_local_rows * n;
        }
        delete[] global_grid;
    } else {
        MPI_Recv(local_data, local_rows * n, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }

    int** current = new int*[local_rows + 2];
    int** next = new int*[local_rows + 2];
    for (int i = 0; i < local_rows + 2; ++i) {
        current[i] = new int[n];
        next[i] = new int[n];
        if (i >= 1 && i <= local_rows) {
            memcpy(current[i], local_data + (i - 1) * n, n * sizeof(int));
        } else {
            memset(current[i], 0, n * sizeof(int));
        }
    }
    delete[] local_data;

    for (int step = 1; step <= t; ++step) {
        auto it = events.find(step);
        if (it != events.end()) {
            const Event& e = it->second;
            if (e.type == 1) {
                int x = e.params[0];
                int y = e.params[1];
                if (x >= start_row && x <= end_row) {
                    int local_x = x - start_row + 1;
                    if (current[local_x][y] == 1) {
                        current[local_x][y] = 2;
                    }
                }
            } else {
                int sx = e.params[0];
                int sy = e.params[1];
                int ex = e.params[2];
                int ey = e.params[3];
                int overlap_sx = std::max(sx, start_row);
                int overlap_ex = std::min(ex, end_row);
                int overlap_sy = std::max(sy, 0);
                int overlap_ey = std::min(ey, n - 1);
                if (overlap_sx <= overlap_ex && overlap_sy <= overlap_ey) {
                    for (int x = overlap_sx; x <= overlap_ex; ++x) {
                        int local_x = x - start_row + 1;
                        for (int y = overlap_sy; y <= overlap_ey; ++y) {
                            if (current[local_x][y] == 3) {
                                current[local_x][y] = 1;
                            }
                        }
                    }
                }
            }
        }

        if (rank > 0) {
            MPI_Sendrecv(current[1], n, MPI_INT, rank - 1, 0,
                        current[0], n, MPI_INT, rank - 1, 0,
                        MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        }
        if (rank < size - 1) {
            MPI_Sendrecv(current[local_rows], n, MPI_INT, rank + 1, 0,
                        current[local_rows + 1], n, MPI_INT, rank + 1, 0,
                        MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        }

        for (int i = 1; i <= local_rows; ++i) {
            for (int j = 0; j < n; ++j) {
                int old_val = current[i][j];
                if (old_val == 2) {
                    next[i][j] = 3;
                } else if (old_val == 1) {
                    int up = (i - 1 >= 0) ? current[i - 1][j] : 0;
                    int down = (i + 1 < local_rows + 2) ? current[i + 1][j] : 0;
                    int left = (j > 0) ? current[i][j - 1] : 0;
                    int right = (j < n - 1) ? current[i][j + 1] : 0;
                    next[i][j] = (up == 2 || down == 2 || left == 2 || right == 2) ? 2 : 1;
                } else {
                    next[i][j] = old_val;
                }
            }
        }

        memcpy(next[0], current[0], n * sizeof(int));
        memcpy(next[local_rows + 1], current[local_rows + 1], n * sizeof(int));
        std::swap(current, next);
    }

    int* sendbuf = new int[local_rows * n];
    for (int i = 0; i < local_rows; ++i) {
        memcpy(sendbuf + i * n, current[i + 1], n * sizeof(int));
    }

    if (rank == 0) {
        int** final_grid = new int*[n];
        for (int i = 0; i < n; ++i) {
            final_grid[i] = new int[n];
        }
        for (int i = 0; i < local_rows; ++i) {
            memcpy(final_grid[start_row + i], sendbuf + i * n, n * sizeof(int));
        }
        delete[] sendbuf;

        for (int src = 1; src < size; ++src) {
            int src_remainder = n % size;
            int src_local_rows = n / size + (src < src_remainder ? 1 : 0);
            int* recvbuf = new int[src_local_rows * n];
            MPI_Recv(recvbuf, src_local_rows * n, MPI_INT, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

            int src_start_row = (src < src_remainder) ? (src * (n / size + 1)) : (src_remainder * (n / size + 1) + (src - src_remainder) * (n / size));
            for (int i = 0; i < src_local_rows; ++i) {
                memcpy(final_grid[src_start_row + i], recvbuf + i * n, n * sizeof(int));
            }
            delete[] recvbuf;
        }

        FILE* fout = fopen(argv[2], "w");
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                if (j != 0) fprintf(fout, " ");
                fprintf(fout, "%d", final_grid[i][j]);
            }
            fprintf(fout, "\n");
        }
        fclose(fout);

        for (int i = 0; i < n; ++i) {
            delete[] final_grid[i];
        }
        delete[] final_grid;
    } else {
        MPI_Send(sendbuf, local_rows * n, MPI_INT, 0, 0, MPI_COMM_WORLD);
        delete[] sendbuf;
    }

    for (int i = 0; i < local_rows + 2; ++i) {
        delete[] current[i];
        delete[] next[i];
    }
    delete[] current;
    delete[] next;

    MPI_Finalize();
    return 0;
}

./flute data.in.bin data.out.bin

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <algorithm>
#include <cassert>
#include <omp.h>
using std::max;
using std::min;

typedef unsigned DTYPE;

typedef struct
{
  DTYPE x, y;   // starting point of the branch
  int n;   // index of neighbor
} Branch;

typedef struct
{
  int deg;   // degree
  DTYPE length;   // total wirelength
  Branch *branch;   // array of tree branches
} Tree;

// Major functions
void readLUT();
DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc);
//Macro: DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
Tree flute(int d, DTYPE x[], DTYPE y[], int acc);
//Macro: Tree flutes(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
DTYPE wirelength(Tree t);
void printtree(Tree t);

#define POWVFILE "POWV9.dat"    // LUT for POWV (Wirelength Vector)
#define PORTFILE "PORT9.dat"    // LUT for PORT (Routing Tree)
#define D 9        // LUT is used for d <= D, D <= 9
#define FLUTEROUTING 1   // 1 to construct routing, 0 to estimate WL only
#define REMOVE_DUPLICATE_PIN 0  // Remove dup. pin for flute_wl() & flute()
#define ACCURACY 8  // Default accuracy is 8
//#define MAXD 2008840  // max. degree of a net that can be handled
#define MAXD 1000  // max. degree of a net that can be handled

// Other useful functions
DTYPE flutes_wl_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
DTYPE flutes_wl_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
Tree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
Tree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
Tree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);

#if REMOVE_DUPLICATE_PIN==1
#else
#define flutes_wl(d, xs, ys, s, acc) flutes_wl_ALLD(d, xs, ys, s, acc)
#define flutes(d, xs, ys, s, acc) flutes_ALLD(d, xs, ys, s, acc)
#endif

#define flutes_wl_ALLD(d, xs, ys, s, acc) flutes_wl_LMD(d, xs, ys, s, acc)
#define flutes_ALLD(d, xs, ys, s, acc) flutes_LMD(d, xs, ys, s, acc)

#define flutes_wl_LMD(d, xs, ys, s, acc)                                \
  (d<=D ? flutes_wl_LD(d, xs, ys, s) : flutes_wl_MD(d, xs, ys, s, acc))
#define flutes_LMD(d, xs, ys, s, acc)                             \
  (d<=D ? flutes_LD(d, xs, ys, s) : flutes_MD(d, xs, ys, s, acc))

#define ADIFF(x,y) ((x)>(y)?(x-y):(y-x))  // Absolute difference

#if D<=7
#elif D==8
#elif D==9
#define MGROUP 362880/4  // Max. # of groups, 9! = 362880
#define MPOWV 79  // Max. # of POWVs per group
#endif
int numgrp[10]={0,0,0,0,6,30,180,1260,10080,90720};

typedef struct point {
  DTYPE x, y;
  int o;     //ZZ: temporarily store discretized x value indices
} POINT;

typedef POINT* POINTptr;

struct csoln
{
  unsigned char parent;
  unsigned char seg[12];  // Add: 0..i, Sub: j..11; seg[i+1]=seg[j-1]=0
  unsigned char row[D-2], col[D-2];
  unsigned char neighbor[2*D-2];
};

struct csoln *LUT[D+1][MGROUP];  // storing 4 .. D
int numsoln[D+1][MGROUP];

void readLUT();
Tree flute(int d, DTYPE x[], DTYPE y[], int acc);
Tree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
Tree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
Tree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
Tree dmergetree(Tree t1, Tree t2);
Tree hmergetree(Tree t1, Tree t2, int s[]);
Tree vmergetree(Tree t1, Tree t2);
void printtree(Tree t);

void readLUT()
{
  static bool LUTread = false;
  if(LUTread) return;
  LUTread = true;

  FILE *fpwv, *fprt;
  struct csoln *p;
  int d, i, j, k, kk, ns, nn, ne;
  unsigned char line[99], *linep, c;
  unsigned char charnum[256];
  unsigned char divd[256], modd[256], div16[256], mod16[256], dsq;

  for (i=0; i<=255; i++) {
    if ('0'<=i && i<='9')
      charnum[i] = i - '0';
    else if (i>='A')
      charnum[i] = i - 'A' + 10;
    else // if (i=='$' || i=='\n' || ... )
      charnum[i] = 0;

    div16[i] = i/16;
    mod16[i] = i%16;
  }

  fpwv=fopen(POWVFILE, "r");
  if (fpwv == NULL) {
    printf("Error in opening %s\n", POWVFILE);
    printf("Please make sure %s and %s\n"
           "are in the current working directory.\n", POWVFILE, PORTFILE);
    exit(1);
  }

#if FLUTEROUTING==1
  fprt=fopen(PORTFILE, "r");
  if (fprt == NULL) {
    printf("Error in opening %s\n", PORTFILE);
    printf("Please make sure %s and %s\n(found in the Flute directory of UMpack)\n"
           "are in the current working directory.\n", POWVFILE, PORTFILE);
    exit(1);
  }
#endif

  for (d=4; d<=D; d++) {
    for (i=0; i<=255; i++) {
      divd[i] = i/d;
      modd[i] = i%d;
    }
    dsq=d*d;

    fscanf(fpwv, "d=%d\n", &d);
#if FLUTEROUTING==1
    fscanf(fprt, "d=%d\n", &d);
#endif
    for (k=0; k<numgrp[d]; k++) {
      ns = (int) charnum[fgetc(fpwv)];

      if (ns==0) {  // same as some previous group
        fscanf(fpwv, "%d\n", &kk);
        numsoln[d][k] = numsoln[d][kk];
        LUT[d][k] = LUT[d][kk];
      }
      else {
        fgetc(fpwv);  // '\n'
        numsoln[d][k] = ns;
        p = (struct csoln*) malloc(ns*sizeof(struct csoln));
        LUT[d][k] = p;
        for (i=1; i<=ns; i++) {
          linep = (unsigned char*)fgets((char*)line, 99, fpwv);
          p->parent = charnum[*(linep++)];
          j = 0;
          while ((p->seg[j++] = charnum[*(linep++)]) != 0) ;
          j = 10;
          while ((p->seg[j--] = charnum[*(linep++)]) != 0) ;
#if FLUTEROUTING==1
          nn = 2*d-2;
          ne = 2*d-3;
          fread(line, 1, d-2, fprt); linep=line;
          for (j=d; j<nn; j++) {
            c = *(linep++);
            if (c >= dsq) {
              c -= dsq;
              p->neighbor[divd[c]] = j;
              ne--;
            }
            p->row[j-d] = divd[c];
            p->col[j-d] = modd[c];
            p->neighbor[j] = j;  // initialized
          }
          fread(line, 1, ne+1, fprt); linep=line;  // last char = \n
          for (j=0; j<ne; j++) {
            c = *(linep++);
            p->neighbor[div16[c]] = mod16[c];
          }
#endif
          p++;
        }
      }
    }
  }
}

static int compare_x(const void *a, const void *b) {
    const POINT *pa = *(const POINTptr*)a;
    const POINT *pb = *(const POINTptr*)b;
    if (pa->x < pb->x) return -1;
    if (pa->x > pb->x) return 1;
    if (pa->y < pb->y) return -1;
    if (pa->y > pb->y) return 1;
    return (pa < pb) ? -1 : 1; // 确保稳定性
}

// 比较函数，用于按y排序ptp数组
static int compare_y(const void *a, const void *b) {
    const POINT *pa = *(const POINTptr*)a;
    const POINT *pb = *(const POINTptr*)b;
    if (pa->y < pb->y) return -1;
    if (pa->y > pb->y) return 1;
    if (pa->x < pb->x) return -1;
    if (pa->x > pb->x) return 1;
    return (pa < pb) ? -1 : 1;
}

Tree flute(int d, DTYPE x[], DTYPE y[], int acc) {
    unsigned allocateSize = MAXD;
    if (d > MAXD) {
        allocateSize = d + 1;
    }

    // 使用对齐的内存分配以提高SIMD效率
    DTYPE *xs, *ys;
    int *s;
    POINT *pt;
    POINTptr *ptp;

    posix_memalign((void **)&xs, 64, sizeof(DTYPE) * allocateSize);
    posix_memalign((void **)&ys, 64, sizeof(DTYPE) * allocateSize);
    posix_memalign((void **)&s, 64, sizeof(int) * allocateSize);
    pt = (POINT *)malloc(sizeof(POINT) * allocateSize);
    ptp = (POINTptr *)malloc(sizeof(POINTptr) * allocateSize);

    for (int i = 0; i < d; i++) {
        pt[i].x = x[i];
        pt[i].y = y[i];
        ptp[i] = &pt[i];
    }

    // 使用快速排序按x排序
    qsort(ptp, d, sizeof(POINTptr), compare_x);

    // 使用SIMD填充xs数组
#pragma omp simd aligned(xs : 64)
    for (int i = 0; i < d; i++) {
        xs[i] = ptp[i]->x;
        ptp[i]->o = i;
    }

    // 使用快速排序按y排序
    qsort(ptp, d, sizeof(POINTptr), compare_y);

    // 使用SIMD填充ys和s数组
#pragma omp simd aligned(ys, s : 64)
    for (int i = 0; i < d; i++) {
        ys[i] = ptp[i]->y;
        s[i] = ptp[i]->o;
    }

    Tree t = flutes(d, xs, ys, s, acc);

    // 释放内存
    free(xs);
    free(ys);
    free(s);
    free(pt);
    free(ptp);

    return t;
}

// xs[] and ys[] are coords in x and y in sorted order
// s[] is a list of nodes in increasing y direction
//   if nodes are indexed in the order of increasing x coord
//   i.e., s[i] = s_i in defined in paper
// The points are (xs[s[i]], ys[i]) for i=0..d-1
//             or (xs[i], ys[si[i]]) for i=0..d-1
// ZZ: RDP denotes Remove Duplicate Pin, simple, like std::unique
Tree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
{
  int i, j, ss;

  for (i=0; i<d-1; i++) {
    if (xs[s[i]]==xs[s[i+1]] && ys[i]==ys[i+1]) {
      if (s[i] < s[i+1])
        ss = s[i+1];
      else {
        ss = s[i];
        s[i] = s[i+1];
      }
      for (j=i+2; j<d; j++) {
        ys[j-1] = ys[j];
        s[j-1] = s[j];
      }
      for (j=ss+1; j<d; j++)
        xs[j-1] = xs[j];
      for (j=0; j<=d-2; j++)
        if (s[j] > ss) s[j]--;
      i--;
      d--;
    }
  }
  return flutes_ALLD(d, xs, ys, s, acc);
}

// For low-degree, i.e., 2 <= d <= D
Tree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[])
{
  int k, pi, i, j;
  struct csoln *rlist, *bestrlist;
  DTYPE dd[2*D-2];  // 0..D-2 for v, D-1..2*D-3 for h
  DTYPE minl, sum;
  DTYPE l[MPOWV+1];
  int hflip;
  Tree t;

  t.deg = d;
  t.branch = (Branch *) malloc((2*d-2)*sizeof(Branch));
  if (d == 2) {
    // printf("s[0]=%d, s[1]=%d, xs[0,1]=[%u,%u], ys[0,1]=[%u,%u]\n",
    //        s[0], s[1], xs[0], xs[1], ys[0], ys[1]);
    minl = xs[1]-xs[0]+ys[1]-ys[0];
    t.branch[0].x = xs[s[0]];
    t.branch[0].y = ys[0];
    t.branch[0].n = 1;
    t.branch[1].x = xs[s[1]];
    t.branch[1].y = ys[1];
    t.branch[1].n = 1;
  }
  else if (d == 3) {
    minl = xs[2]-xs[0]+ys[2]-ys[0];
    t.branch[0].x = xs[s[0]];
    t.branch[0].y = ys[0];
    t.branch[0].n = 3;
    t.branch[1].x = xs[s[1]];
    t.branch[1].y = ys[1];
    t.branch[1].n = 3;
    t.branch[2].x = xs[s[2]];
    t.branch[2].y = ys[2];
    t.branch[2].n = 3;
    t.branch[3].x = xs[1];
    t.branch[3].y = ys[1];
    t.branch[3].n = 3;
  }
  else {
    k = 0;
    if (s[0] < s[2]) k++;
    if (s[1] < s[2]) k++;

    for (i=3; i<=d-1; i++) {  // p0=0 always, skip i=1 for symmetry
      pi = s[i];
      for (j=d-1; j>i; j--)
        if (s[j] < s[i])
          pi--;
      k = pi + (i+1)*k;
    }

    if (k < numgrp[d]) { // no horizontal flip
      hflip = 0;
      for (i=1; i<=d-3; i++) {
        dd[i]=ys[i+1]-ys[i];
        dd[d-1+i]=xs[i+1]-xs[i];
      }
    }
    else {
      hflip = 1;
      k=2*numgrp[d]-1-k;
      for (i=1; i<=d-3; i++) {
        dd[i]=ys[i+1]-ys[i];
        dd[d-1+i]=xs[d-1-i]-xs[d-2-i];
      }
    }

    minl = l[0] = xs[d-1]-xs[0]+ys[d-1]-ys[0];
    rlist = LUT[d][k];
    for (i=0; rlist->seg[i]>0; i++)
      minl += dd[rlist->seg[i]];
    bestrlist = rlist;
    l[1] = minl;
    j = 2;
    while (j <= numsoln[d][k]) {
      rlist++;
      sum = l[rlist->parent];
      for (i=0; rlist->seg[i]>0; i++)
        sum += dd[rlist->seg[i]];
      for (i=10; rlist->seg[i]>0; i--)
        sum -= dd[rlist->seg[i]];
      if (sum < minl) {
        minl = sum;
        bestrlist = rlist;
      }
      l[j++] = sum;
    }

    t.branch[0].x = xs[s[0]];
    t.branch[0].y = ys[0];
    t.branch[1].x = xs[s[1]];
    t.branch[1].y = ys[1];
    for (i=2; i<d-2; i++) {
      t.branch[i].x = xs[s[i]];
      t.branch[i].y = ys[i];
      t.branch[i].n = bestrlist->neighbor[i];
    }
    t.branch[d-2].x = xs[s[d-2]];
    t.branch[d-2].y = ys[d-2];
    t.branch[d-1].x = xs[s[d-1]];
    t.branch[d-1].y = ys[d-1];
    if (hflip) {
      if (s[1] < s[0]) {
        t.branch[0].n = bestrlist->neighbor[1];
        t.branch[1].n = bestrlist->neighbor[0];
      }
      else {
        t.branch[0].n = bestrlist->neighbor[0];
        t.branch[1].n = bestrlist->neighbor[1];
      }
      if (s[d-1] < s[d-2]) {
        t.branch[d-2].n = bestrlist->neighbor[d-1];
        t.branch[d-1].n = bestrlist->neighbor[d-2];
      }
      else {
        t.branch[d-2].n = bestrlist->neighbor[d-2];
        t.branch[d-1].n = bestrlist->neighbor[d-1];
      }
      for (i=d; i<2*d-2; i++) {
        t.branch[i].x = xs[d-1-bestrlist->col[i-d]];
        t.branch[i].y = ys[bestrlist->row[i-d]];
        t.branch[i].n = bestrlist->neighbor[i];
      }
    }
    else {  // !hflip
      if (s[0] < s[1]) {
        t.branch[0].n = bestrlist->neighbor[1];
        t.branch[1].n = bestrlist->neighbor[0];
      }
      else {
        t.branch[0].n = bestrlist->neighbor[0];
        t.branch[1].n = bestrlist->neighbor[1];
      }
      if (s[d-2] < s[d-1]) {
        t.branch[d-2].n = bestrlist->neighbor[d-1];
        t.branch[d-1].n = bestrlist->neighbor[d-2];
      }
      else {
        t.branch[d-2].n = bestrlist->neighbor[d-2];
        t.branch[d-1].n = bestrlist->neighbor[d-1];
      }
      for (i=d; i<2*d-2; i++) {
        t.branch[i].x = xs[bestrlist->col[i-d]];
        t.branch[i].y = ys[bestrlist->row[i-d]];
        t.branch[i].n = bestrlist->neighbor[i];
      }
    }
  }
  t.length = minl;

  return t;
}

// For medium-degree, i.e., D+1 <= d <= D2
Tree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
{
  unsigned allocateSize = MAXD;
  if (d > MAXD)
    allocateSize = d+1;
  DTYPE *x1      = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);
  DTYPE *x2      = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);
  DTYPE *y1      = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);
  DTYPE *y2      = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);
  int   *si      = (int*)   malloc(sizeof(int)*allocateSize);
  int   *s1      = (int*)   malloc(sizeof(int)*allocateSize);
  int   *s2      = (int*)   malloc(sizeof(int)*allocateSize);
  float *score   = (float*) malloc(sizeof(float)*2*allocateSize);
  float *penalty = (float*) malloc(sizeof(float)*allocateSize);
  DTYPE *distx   = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);
  DTYPE *disty   = (DTYPE*) malloc(sizeof(DTYPE)*allocateSize);

//  DTYPE x1[MAXD], x2[MAXD], y1[MAXD], y2[MAXD];
//  int si[MAXD], s1[MAXD], s2[MAXD];
//  float score[2*MAXD], penalty[MAXD];
//  DTYPE distx[MAXD], disty[MAXD];
  float pnlty, dx, dy;
  DTYPE ll, minl, coord1, coord2;
  int i, r, p, maxbp, bestbp, bp, nbp, ub, lb, n1, n2, nn1, nn2, ms, newacc;
  Tree t, t1, t2, bestt1, bestt2;
  int mins, maxs, minsi, maxsi;
  DTYPE xydiff;

  #pragma omp simd
  for(i = 0; i < (int)allocateSize; ++i) {
    distx[i] = disty[i] = 0;
  }

  if (s[0] < s[d-1]) {
    ms = max(s[0], s[1]);
    for (i=2; i<=ms; i++)
      ms = max(ms, s[i]);
    if (ms <= d-3) {
      for (i=0; i<=ms; i++) {
        x1[i] = xs[i];
        y1[i] = ys[i];
        s1[i] = s[i];
      }
      x1[ms+1] = xs[ms];
      y1[ms+1] = ys[ms];
      s1[ms+1] = ms+1;

      s2[0] = 0;
      for (i=1; i<=d-1-ms; i++)
        s2[i] = s[i+ms]-ms;

      t1 = flutes_LMD(ms+2, x1, y1, s1, acc);
      t2 = flutes_LMD(d-ms, xs+ms, ys+ms, s2, acc);
      t = dmergetree(t1, t2);
      free(t1.branch);
      free(t2.branch);

      free(x1);
      free(x2);
      free(y1);
      free(y2);
      free(si);
      free(s1);
      free(s2);
      free(score);
      free(penalty);
      free(distx);
      free(disty);
      return t;
    }
  }
  else {  // (s[0] > s[d-1])
    ms = min(s[0], s[1]);
    for (i=2; i<=d-1-ms; i++)
      ms = min(ms, s[i]);
    if (ms >= 2) {
      x1[0] = xs[ms];
      y1[0] = ys[0];
      s1[0] = s[0]-ms+1;
      for (i=1; i<=d-1-ms; i++) {
        x1[i] = xs[i+ms-1];
        y1[i] = ys[i];
        s1[i] = s[i]-ms+1;
      }
      x1[d-ms] = xs[d-1];
      y1[d-ms] = ys[d-1-ms];
      s1[d-ms] = 0;

      s2[0] = ms;
      for (i=1; i<=ms; i++)
        s2[i] = s[i+d-1-ms];

      t1 = flutes_LMD(d+1-ms, x1, y1, s1, acc);
      t2 = flutes_LMD(ms+1, xs, ys+d-1-ms, s2, acc);
      t = dmergetree(t1, t2);
      free(t1.branch);
      free(t2.branch);

      free(x1);
      free(x2);
      free(y1);
      free(y2);
      free(si);
      free(s1);
      free(s2);
      free(score);
      free(penalty);
      free(distx);
      free(disty);
      return t;
    }
  }

  // Find inverse si[] of s[]
  for (r=0; r<d; r++)
    si[s[r]] = r;

  // Determine breaking directions and positions dp[]
  lb=max((d-acc)/5, 2);
  ub=d-1-lb;

  // Compute scores
#define AA 0.6  // 2.0*BB
#define BB 0.3
  //#define CCC max(0.425-0.005*d-0.015*acc, 0.1)
  //#define CCC max(0.43-0.005*d-0.01*acc, 0.1)
#define CCC max(0.41-0.005*d, 0.1)
#define DD 4.8
  float DDD = DD/(d-1);

  // Compute penalty[]
  dx = CCC*(xs[d-2]-xs[1])/(d-3);
  dy = CCC*(ys[d-2]-ys[1])/(d-3);
  for (r = d/2, pnlty = 0; r>=2; r--, pnlty += dx)
    penalty[r] = pnlty,  penalty[d-1-r] = pnlty;
  penalty[1] = pnlty,  penalty[d-2] = pnlty;
  penalty[0] = pnlty,  penalty[d-1] = pnlty;
  for (r = d/2-1, pnlty = dy; r>=2; r--, pnlty += dy)
    penalty[s[r]] += pnlty,  penalty[s[d-1-r]] += pnlty;
  penalty[s[1]] += pnlty,  penalty[s[d-2]] += pnlty;
  penalty[s[0]] += pnlty,  penalty[s[d-1]] += pnlty;
  //#define CC 0.16
  //#define v(r) ((r==0||r==1||r==d-2||r==d-1) ? d-3 : abs(d-1-r-r))
  //    for (r=0; r<d; r++)
  //        penalty[r] = v(r)*dx + v(si[r])*dy;

  // Compute distx[], disty[]
  xydiff = (xs[d-1] - xs[0]) - (ys[d-1] - ys[0]);
  if (s[0] < s[1])
    mins = s[0], maxs = s[1];
  else mins = s[1], maxs = s[0];
  if (si[0] < si[1])
    minsi = si[0], maxsi = si[1];
  else minsi = si[1], maxsi = si[0];
  for (r=2; r<=ub; r++) {
    if (s[r] < mins)
      mins = s[r];
    else if (s[r] > maxs)
      maxs = s[r];
    distx[r] = xs[maxs] - xs[mins];
    if (si[r] < minsi)
      minsi = si[r];
    else if (si[r] > maxsi)
      maxsi = si[r];
    disty[r] = ys[maxsi] - ys[minsi] + xydiff;
  }

  if (s[d-2] < s[d-1])
    mins = s[d-2], maxs = s[d-1];
  else mins = s[d-1], maxs = s[d-2];
  if (si[d-2] < si[d-1])
    minsi = si[d-2], maxsi = si[d-1];
  else minsi = si[d-1], maxsi = si[d-2];
  #pragma omp simd
  for (r=d-3; r>=lb; r--) {
    if (s[r] < mins)
      mins = s[r];
    else if (s[r] > maxs)
      maxs = s[r];
    distx[r] += xs[maxs] - xs[mins];
    if (si[r] < minsi)
      minsi = si[r];
    else if (si[r] > maxsi)
      maxsi = si[r];
    disty[r] += ys[maxsi] - ys[minsi];
  }

  nbp=0;
  #pragma omp simd
  for (r=lb; r<=ub; r++) {
    if (si[r]<=1)
      score[nbp] = (xs[r+1] - xs[r-1]) - penalty[r]
        - AA*(ys[2]-ys[1]) - DDD*disty[r];
    else if (si[r]>=d-2)
      score[nbp] = (xs[r+1] - xs[r-1]) - penalty[r]
        - AA*(ys[d-2]-ys[d-3]) - DDD*disty[r];
    else score[nbp] = (xs[r+1] - xs[r-1]) - penalty[r]
           - BB*(ys[si[r]+1]-ys[si[r]-1]) - DDD*disty[r];
    nbp++;

    if (s[r]<=1)
      score[nbp] = (ys[r+1] - ys[r-1]) - penalty[s[r]]
        - AA*(xs[2]-xs[1]) - DDD*distx[r];
    else if (s[r]>=d-2)
      score[nbp] = (ys[r+1] - ys[r-1]) - penalty[s[r]]
        - AA*(xs[d-2]-xs[d-3]) - DDD*distx[r];
    else score[nbp] = (ys[r+1] - ys[r-1]) - penalty[s[r]]
           - BB*(xs[s[r]+1]-xs[s[r]-1]) - DDD*distx[r];
    nbp++;
  }

  if (acc <= 3)
    newacc = 1;
  else {
    newacc = acc/2;
    if (acc >= nbp) acc = nbp-1;
  }

  minl = (DTYPE) INT_MAX;
  bestt1.branch = bestt2.branch = NULL;
  for (i=0; i<acc; i++) {
    maxbp = 0;
    for (bp=1; bp<nbp; bp++)
      if (score[maxbp] < score[bp]) maxbp = bp;
    score[maxbp] = -9e9f;

#define BreakPt(bp) ((bp)/2+lb)
#define BreakInX(bp) ((bp)%2==0)
    p = BreakPt(maxbp);
    // Breaking in p
    if (BreakInX(maxbp)) {  // break in x
      n1 = n2 = 0;
      for (r=0; r<d; r++) {
        if (s[r] < p) {
          s1[n1] = s[r];
          y1[n1] = ys[r];
          n1++;
        }
        else if (s[r] > p) {
          s2[n2] = s[r]-p;
          y2[n2] = ys[r];
          n2++;
        }
        else { // if (s[r] == p)  i.e.,  r = si[p]
          s1[n1] = p;  s2[n2] = 0;
          y1[n1] = y2[n2] = ys[r];
          nn1 = n1;  nn2 = n2;
          n1++;  n2++;
        }
      }

      t1 = flutes_LMD(p+1, xs, y1, s1, newacc);
      t2 = flutes_LMD(d-p, xs+p, y2, s2, newacc);
      ll = t1.length + t2.length;
      coord1 = t1.branch[t1.branch[nn1].n].y;
      coord2 = t2.branch[t2.branch[nn2].n].y;
      if (t2.branch[nn2].y > max(coord1, coord2))
        ll -= t2.branch[nn2].y - max(coord1, coord2);
      else if (t2.branch[nn2].y < min(coord1, coord2))
        ll -= min(coord1, coord2) - t2.branch[nn2].y;
    }
    else {  // if (!BreakInX(maxbp))
      n1 = n2 = 0;
      for (r=0; r<d; r++) {
        if (si[r] < p) {
          s1[si[r]] = n1;
          x1[n1] = xs[r];
          n1++;
        }
        else if (si[r] > p) {
          s2[si[r]-p] = n2;
          x2[n2] = xs[r];
          n2++;
        }
        else { // if (si[r] == p)  i.e.,  r = s[p]
          s1[p] = n1;  s2[0] = n2;
          x1[n1] = x2[n2] = xs[r];
          n1++;  n2++;
        }
      }

      t1 = flutes_LMD(p+1, x1, ys, s1, newacc);
      t2 = flutes_LMD(d-p, x2, ys+p, s2, newacc);
      ll = t1.length + t2.length;
      coord1 = t1.branch[t1.branch[p].n].x;
      coord2 = t2.branch[t2.branch[0].n].x;
      if (t2.branch[0].x > max(coord1, coord2))
        ll -= t2.branch[0].x - max(coord1, coord2);
      else if (t2.branch[0].x < min(coord1, coord2))
        ll -= min(coord1, coord2) - t2.branch[0].x;
    }
    if (minl > ll) {
      minl = ll;
      free(bestt1.branch);
      free(bestt2.branch);
      bestt1 = t1;
      bestt2 = t2;
      bestbp = maxbp;
    }
    else {
      free(t1.branch);
      free(t2.branch);
    }
  }

  if (BreakInX(bestbp))
    t = hmergetree(bestt1, bestt2, s);
  else t = vmergetree(bestt1, bestt2);
  free(bestt1.branch);
  free(bestt2.branch);

  free(x1);
  free(x2);
  free(y1);
  free(y2);
  free(si);
  free(s1);
  free(s2);
  free(score);
  free(penalty);
  free(distx);
  free(disty);
  return t;
}

Tree dmergetree(Tree t1, Tree t2)
{
  int i, d, prev, curr, next, offset1, offset2;
  Tree t;

  t.deg = d = t1.deg + t2.deg - 2;
  t.length = t1.length + t2.length;
  t.branch = (Branch *) malloc((2*d-2)*sizeof(Branch));
  offset1 = t2.deg-2;
  offset2 = 2*t1.deg-4;

  for (i=0; i<=t1.deg-2; i++) {
    t.branch[i].x = t1.branch[i].x;
    t.branch[i].y = t1.branch[i].y;
    t.branch[i].n = t1.branch[i].n + offset1;
  }
  for (i=t1.deg-1; i<=d-1; i++) {
    t.branch[i].x = t2.branch[i-t1.deg+2].x;
    t.branch[i].y = t2.branch[i-t1.deg+2].y;
    t.branch[i].n = t2.branch[i-t1.deg+2].n + offset2;
  }
  for (i=d; i<=d+t1.deg-3; i++) {
    t.branch[i].x = t1.branch[i-offset1].x;
    t.branch[i].y = t1.branch[i-offset1].y;
    t.branch[i].n = t1.branch[i-offset1].n + offset1;
  }
  for (i=d+t1.deg-2; i<=2*d-3; i++) {
    t.branch[i].x = t2.branch[i-offset2].x;
    t.branch[i].y = t2.branch[i-offset2].y;
    t.branch[i].n = t2.branch[i-offset2].n + offset2;
  }

  prev = t2.branch[0].n + offset2;
  curr = t1.branch[t1.deg-1].n + offset1;
  next = t.branch[curr].n;
  while (curr != next) {
    t.branch[curr].n = prev;
    prev = curr;
    curr = next;
    next = t.branch[curr].n;
  }
  t.branch[curr].n = prev;

  return t;
}

Tree hmergetree(Tree t1, Tree t2, int s[])
{
  int i, prev, curr, next, extra, offset1, offset2;
  int p, ii, n1, n2, nn1, nn2;
  DTYPE coord1, coord2;
  Tree t;

  t.deg = t1.deg + t2.deg - 1;
  t.length = t1.length + t2.length;
  t.branch = (Branch *) malloc((2*t.deg-2)*sizeof(Branch));
  offset1 = t2.deg-1;
  offset2 = 2*t1.deg-3;

  p = t1.deg - 1;
  n1 = n2 = 0;
  for (i=0; i<t.deg; i++) {
    if (s[i] < p) {
      t.branch[i].x = t1.branch[n1].x;
      t.branch[i].y = t1.branch[n1].y;
      t.branch[i].n = t1.branch[n1].n + offset1;
      n1++;
    }
    else if (s[i] > p) {
      t.branch[i].x = t2.branch[n2].x;
      t.branch[i].y = t2.branch[n2].y;
      t.branch[i].n = t2.branch[n2].n + offset2;
      n2++;
    }
    else {
      t.branch[i].x = t2.branch[n2].x;
      t.branch[i].y = t2.branch[n2].y;
      t.branch[i].n = t2.branch[n2].n + offset2;
      nn1 = n1;  nn2 = n2;  ii = i;
      n1++;  n2++;
    }
  }
  for (i=t.deg; i<=t.deg+t1.deg-3; i++) {
    t.branch[i].x = t1.branch[i-offset1].x;
    t.branch[i].y = t1.branch[i-offset1].y;
    t.branch[i].n = t1.branch[i-offset1].n + offset1;
  }
  for (i=t.deg+t1.deg-2; i<=2*t.deg-4; i++) {
    t.branch[i].x = t2.branch[i-offset2].x;
    t.branch[i].y = t2.branch[i-offset2].y;
    t.branch[i].n = t2.branch[i-offset2].n + offset2;
  }
  extra = 2*t.deg-3;
  coord1 = t1.branch[t1.branch[nn1].n].y;
  coord2 = t2.branch[t2.branch[nn2].n].y;
  if (t2.branch[nn2].y > max(coord1, coord2)) {
    t.branch[extra].y = max(coord1, coord2);
    t.length -= t2.branch[nn2].y - t.branch[extra].y;
  }
  else if (t2.branch[nn2].y < min(coord1, coord2)) {
    t.branch[extra].y = min(coord1, coord2);
    t.length -= t.branch[extra].y - t2.branch[nn2].y;
  }
  else t.branch[extra].y = t2.branch[nn2].y;
  t.branch[extra].x = t2.branch[nn2].x;
  t.branch[extra].n = t.branch[ii].n;
  t.branch[ii].n = extra;

  prev = extra;
  curr = t1.branch[nn1].n + offset1;
  next = t.branch[curr].n;
  while (curr != next) {
    t.branch[curr].n = prev;
    prev = curr;
    curr = next;
    next = t.branch[curr].n;
  }
  t.branch[curr].n = prev;

  return t;
}

Tree vmergetree(Tree t1, Tree t2)
{
  int i, prev, curr, next, extra, offset1, offset2;
  DTYPE coord1, coord2;
  Tree t;

  t.deg = t1.deg + t2.deg - 1;
  t.length = t1.length + t2.length;
  t.branch = (Branch *) malloc((2*t.deg-2)*sizeof(Branch));
  offset1 = t2.deg-1;
  offset2 = 2*t1.deg-3;

  for (i=0; i<=t1.deg-2; i++) {
    t.branch[i].x = t1.branch[i].x;
    t.branch[i].y = t1.branch[i].y;
    t.branch[i].n = t1.branch[i].n + offset1;
  }
  for (i=t1.deg-1; i<=t.deg-1; i++) {
    t.branch[i].x = t2.branch[i-t1.deg+1].x;
    t.branch[i].y = t2.branch[i-t1.deg+1].y;
    t.branch[i].n = t2.branch[i-t1.deg+1].n + offset2;
  }
  for (i=t.deg; i<=t.deg+t1.deg-3; i++) {
    t.branch[i].x = t1.branch[i-offset1].x;
    t.branch[i].y = t1.branch[i-offset1].y;
    t.branch[i].n = t1.branch[i-offset1].n + offset1;
  }
  for (i=t.deg+t1.deg-2; i<=2*t.deg-4; i++) {
    t.branch[i].x = t2.branch[i-offset2].x;
    t.branch[i].y = t2.branch[i-offset2].y;
    t.branch[i].n = t2.branch[i-offset2].n + offset2;
  }
  extra = 2*t.deg-3;
  coord1 = t1.branch[t1.branch[t1.deg-1].n].x;
  coord2 = t2.branch[t2.branch[0].n].x;
  if (t2.branch[0].x > max(coord1, coord2)) {
    t.branch[extra].x = max(coord1, coord2);
    t.length -= t2.branch[0].x - t.branch[extra].x;
  }
  else if (t2.branch[0].x < min(coord1, coord2)) {
    t.branch[extra].x = min(coord1, coord2);
    t.length -= t.branch[extra].x - t2.branch[0].x;
  }
  else t.branch[extra].x = t2.branch[0].x;
  t.branch[extra].y = t2.branch[0].y;
  t.branch[extra].n = t.branch[t1.deg-1].n;
  t.branch[t1.deg-1].n = extra;

  prev = extra;
  curr = t1.branch[t1.deg-1].n + offset1;
  next = t.branch[curr].n;
  while (curr != next) {
    t.branch[curr].n = prev;
    prev = curr;
    curr = next;
    next = t.branch[curr].n;
  }
  t.branch[curr].n = prev;

  return t;
}

int main(int argc, const char **argv) {
  if(argc != 3) {
    printf("usage: %s <input.bin> <output.bin>\n", argv[0]);
    return 255;
  }
  readLUT();

  FILE *input_fp = fopen(argv[1], "rb");
  assert(input_fp);
  fseek(input_fp, 0, SEEK_END);
  unsigned long input_bin_size = ftell(input_fp);
  fseek(input_fp, 0, SEEK_SET);
  int *buf_input = (int *)malloc(input_bin_size);
  assert(fread(buf_input, 1, input_bin_size, input_fp) == input_bin_size);
  fclose(input_fp);

  int num_nets = buf_input[0];
  int *pin_st = buf_input + 1;
  int tot_num_pins = pin_st[num_nets];
  DTYPE *data_x = (DTYPE *)(pin_st + num_nets + 1);
  DTYPE *data_y = (DTYPE *)(data_x + tot_num_pins);

  DTYPE *result_x = (DTYPE *)calloc(tot_num_pins * 2, sizeof(DTYPE));
  DTYPE *result_y = (DTYPE *)calloc(tot_num_pins * 2, sizeof(DTYPE));
  int *result_n = (int *)calloc(tot_num_pins * 2, sizeof(int));

  #pragma omp parallel for
  for(int net_i = 0; net_i < num_nets; ++net_i) {
    Tree t = flute(pin_st[net_i + 1] - pin_st[net_i],
                   data_x + pin_st[net_i],
                   data_y + pin_st[net_i],
                   ACCURACY);
    // printf("%d %d\n", t.deg, t.length);
    // for(int i = 0; i < 2 * t.deg - 2; ++i) {
    //   printf("%d %d %d\n", t.branch[i].x, t.branch[i].y, t.branch[i].n);
    // }
    int offset = pin_st[net_i] * 2;
    for(int i = 0; i < 2 * t.deg - 2; ++i) {
      result_x[offset + i] = t.branch[i].x;
      result_y[offset + i] = t.branch[i].y;
      result_n[offset + i] = t.branch[i].n;
    }
    free(t.branch);
  }

  FILE *output_fp = fopen(argv[2], "wb");
  assert(output_fp);
  fwrite(result_x, sizeof(DTYPE), tot_num_pins, output_fp);
  fwrite(result_y, sizeof(DTYPE), tot_num_pins, output_fp);
  fwrite(result_n, sizeof(int), tot_num_pins, output_fp);
  fclose(output_fp);

  free(result_x);
  free(result_y);
  free(result_n);
  free(buf_input);
  return 0;
}