c++ - 计算文件中的 react 数

Question

我正在尝试使用 MPI c/c++ 和 OpenMP c/c++ 库来计算给定数据集中有多少个矩形。数据集是一个包含 1 和 0 的文本文件，它非常大，大约有 100000x100000 行 0 和 1，但为了简单起见，我将提供该集中的数据示例。

我最初想将数据放入整数 2d 数组，这里是我的数据集的样子

0 0 0 0 1 1 1 00 0 0 0 1 1 1 00 1 1 1 0 0 0 00 1 1 1 0 0 0 00 1 1 1 0 1 1 00 0 0 0 0 1 1 00 0 0 1 1 1 1 10 0 0 1 1 1 1 1

我应该得到 4 个矩形的坐标。

1. 右上角的那个(2 x 3 矩形)
2. 左边那个(3x3 矩形)
3. 2个重叠的矩形

需要的每个矩形的坐标是:左上角&右下角

考虑到使用 mpi/openmp 库来提高效率，我该如何解决这个问题，我是否需要在线程上划分矩阵？

Answer 1

我认为以下内容可以完成这项工作。请注意，如果我没记错的话，这个(新的，请参阅历史记录)算法应该缩放 O(N)。

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

//--some helper functions and structs--
typedef struct Line_s {
    long x0;
    long nx;
    long y0;
} Line;

typedef struct LineVect_s {
  Line *lines;
  long size;
  long capacity;
} LineVect;

LineVect* newLineVect(void) {
    LineVect *vect = malloc(sizeof(LineVect));
    vect->lines = malloc(sizeof(Line));
    vect->size = 0;
    vect->capacity = 1;
    return vect;
}

void deleteLineVect(LineVect* self) {
    free(self->lines);
    free(self);
}

void LineVect_pushBack(LineVect* self, Line line) {
    self->size++;
    if (self->size > self->capacity) {
        self->capacity *= 2;
        self->lines = realloc(self->lines, self->capacity * sizeof(Line));
    }
    self->lines[self->size - 1] = line;
}

typedef struct Block_s {
    long x0, y0;
    long nx, ny;
    char val;
} Block;

typedef struct BlockVect_s {
   Block *blocks;
   long size;
   long capacity;
} BlockVect;

BlockVect* newBlockVect(void) {
    BlockVect *vect = malloc(sizeof(BlockVect));
    vect->blocks = malloc(sizeof(Block));
    vect->size = 0;
    vect->capacity = 1;
    return vect;
}

void deleteBlockVect(BlockVect* self) {
    free(self->blocks);
    free(self);
}

void BlockVect_pushBack(BlockVect* self, Block block) {
    self->size++;
    if (self->size > self->capacity) {
        self->capacity *= 2;
        self->blocks = realloc(self->blocks, self->capacity * sizeof(Block));
    }
    self->blocks[self->size - 1] = block;
}


LineVect** find_lines(char* data, long nx, long ny) {
    LineVect **slices = malloc(ny * sizeof(LineVect*));
    //essentially each y-slice is independent here
#pragma omp parallel for
    for (long y = 0; y < ny; y++) {
        slices[y] = newLineVect();
        char prev_val = 0;
        long xstart = 0;
        for (long x = 0; x < nx; x++) {
            char val = data[nx*y + x];
            if (val == 1 && prev_val == 0) {
                xstart = x;
            } else if (val == 0 && prev_val == 1) {
                Line line;
                line.y0 = -1;
                line.x0 = xstart;
                line.nx = x - xstart;
//              printf("Found line at y=%d from x=%d to %d\n", y, xstart, x-1); 
                LineVect_pushBack(slices[y], line);
            }
            
            if (val == 1 && x == nx-1) {
                Line line;
                line.y0 = -1;
                line.x0 = xstart;
                line.nx = x - xstart + 1;
//              printf("Found line at y=%d from x=%d to %d\n", y, xstart, x);
                LineVect_pushBack(slices[y], line);
            }
            prev_val = val;
        }
    }
    return slices;
}

BlockVect* find_blocks(LineVect** slices, long ny) {
    BlockVect* blocks = newBlockVect();
    //for each line
    for (long y = 0; y < ny; y++) {
        LineVect* slice = slices[y];
        long l2 = 0;
        for (long l = 0; l < slice->size; l++) {
            Line *line = slice->lines + l;
            if (line->y0 == -1) {
                line->y0 = y;
            }
            char match = 0;
            //check next slice if there is any
            if (y != ny-1) {
                //check all the remaining lines in the next slice
                if (l2 < slices[y+1]->size) {
                    Line *line2 = slices[y+1]->lines + l2;
                    //note that we exploit the sorted nature of the lines (x0 must increase with l2)
                    for (; l2 < slices[y+1]->size && line2->x0 < line->x0; l2++) {
                        line2 = slices[y+1]->lines + l2;
                    }
                    //matching line2 found -> mark it as same block (y0 cascades)
                    if (line->x0 == line2->x0 && line->nx == line2->nx) {
                        match = 1;
                        line2->y0 = line->y0;
                    }
                }
            }
            //current line didn't find a matching line hence store it (and all the previously collected lines) as a block
            if (match == 0) {
                Block block;
                block.x0 = line->x0;
                block.nx = line->nx;
                block.y0 = line->y0;
                block.ny = y - line->y0 + 1;
                BlockVect_pushBack(blocks, block);
            }
        }
    }
    return blocks;
}

#define Nx 30000
#define Ny 30000

char * write_blocks(const BlockVect* blocks) {
    char * data = calloc(Ny, Nx*sizeof(char));
    for (long i = 0; i < blocks->size; i++) {
        Block *block = blocks->blocks + i;
        for (long y = block->y0; y < block->y0 + block->ny; y++) {
            for (long x = block->x0; x < block->x0 + block->nx; x++) {
                data[Nx*y + x] = 1;
            }
        }
    }
    return data;
}

int main() {
    struct timeval t0;
    gettimeofday(&t0, NULL);
    
//  char data[Ny][Nx] = {
//      {0, 0, 0, 0, 1, 1, 1, 0},
//      {0, 0, 0, 0, 1, 1, 1, 0},
//      {0, 1, 1, 1, 0, 0, 0, 1},
//      {0, 1, 1, 1, 0, 0, 0, 1},
//      {0, 1, 1, 1, 0, 1, 1, 0},
//      {0, 0, 0, 0, 0, 1, 1, 0},
//      {0, 1, 0, 1, 1, 1, 1, 1},
//      {0, 0, 0, 1, 1, 1, 1, 1}
//  };
    
    srand(t0.tv_usec);
    char * input = calloc(Ny, Nx*sizeof(char));
    printf("Filling...\n");
    for (long y = 0; y < Ny; y++) {
        for (long x = 0; x < Nx; x++) {
            input[Nx*y + x] = rand() % 10 != 0; //data[y][x];
        }
    }
    printf("Computing...\n");
    struct timeval t;
    struct timeval t_new;
    gettimeofday(&t, NULL);
    
    LineVect** slices = find_lines(input, Nx, Ny);
    
    gettimeofday(&t_new, NULL);
    printf("Finding lines took %lu milliseconds\n", (t_new.tv_sec - t.tv_sec)*1000 + (t_new.tv_usec -t.tv_usec) / 1000);
    
    gettimeofday(&t, NULL);
    
    BlockVect* blocks = find_blocks(slices, Ny);
    
    gettimeofday(&t_new, NULL);
    printf("Finding blocks took %lu milliseconds\n", (t_new.tv_sec - t.tv_sec)*1000 + (t_new.tv_usec -t.tv_usec) / 1000);
    
    long n_lines = 0;
    for (long y = 0; y < Ny; y++) {
        n_lines += slices[y]->size;
        deleteLineVect(slices[y]);
    }
    free(slices);
    
    printf("Done computation of %ld lines and %ld blocks\n", n_lines, blocks->size);
    
    char * output = write_blocks(blocks);
    deleteBlockVect(blocks);
    
    char pass = 1;
    for (long y = 0; y < Ny; y++) {
        for (long x = 0; x < Nx; x++) {
            if (input[Nx*y + x] != output[Nx*y + x]) {
                printf("ERROR at x=%ld and y=%ld!\n", x, y);
                pass = 0;
            }
        }
    }
    printf("Plausibility check %s\n", pass ? "PASSED" : "FAILED");
    
    //free all the rest
    free(input);
    free(output);
    gettimeofday(&t_new, NULL);
    printf("Whole run took %.3lf seconds\n", (t_new.tv_sec - t0.tv_sec) + (t_new.tv_usec -t0.tv_usec) * 1e-6);
}

---

缩放测试:

---

在我的笔记本电脑上(Intel(R) Core(TM) i5-4310U @ 2.0 GHz，2 个物理内核):

OMP_NUM_THREADS=1

Filling...
Computing...
Finding lines took 2973 milliseconds
Finding blocks took 2035 milliseconds
Done computation of 81012713 lines and 80743987 blocks
Plausibility check PASSED
Whole run took 16.981 seconds

OMP_NUM_THREADS=2

Filling...
Computing...
Finding lines took 1725 milliseconds
Finding blocks took 2019 milliseconds
Done computation of 81027281 lines and 80757086 blocks
Plausibility check PASSED
Whole run took 15.392 seconds

---

在我的电脑上(Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz，4 个物理内核):

OMP_NUM_THREADS=1

Filling...
Computing...
Finding lines took 3656 milliseconds
Finding blocks took 1631 milliseconds
Done computation of 81024019 lines and 80754472 blocks
Plausibility check PASSED
Whole run took 13.611 seconds

OMP_NUM_THREADS=4

Filling...
Computing...
Finding lines took 1007 milliseconds
Finding blocks took 1637 milliseconds
Done computation of 81036637 lines and 80766687 blocks
Plausibility check PASSED
Whole run took 11.055 seconds

OMP_NUM_THREADS=8

Filling...
Computing...
Finding lines took 812 milliseconds
Finding blocks took 1655 milliseconds
Done computation of 81025147 lines and 80754509 blocks
Plausibility check PASSED
Whole run took 11.137 seconds

这么看来，平行部分的缩放(寻线)还是相当不错的。运行的大部分时间是通过随机数(填充)实际生成样本数据。相比之下，rects 的查找速度相当快。