高性能计算

卷帘矩阵划分和顺序划分

卷帘矩阵

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#define R_SIZE 1024

void initialize_matrices(int *a, int *b, int size) {
    for (int i = 0; i < size; i++) {
        a[i] = rand() % 10000;
        b[i] = rand() % 10000;
    }
}

int main(int argc, char *argv[]) {
    int rank, size;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    int block_size = R_SIZE / size;
    if (R_SIZE % size != 0) {
        if (rank == 0) {
            printf("R_SIZE must be divisible by size\n");
        }
        MPI_Finalize();
        return 1;
    }

    int *a = (int *)malloc(R_SIZE * R_SIZE * sizeof(int));
    int *b = (int *)malloc(R_SIZE * R_SIZE * sizeof(int));
    int *c = (int *)malloc(R_SIZE * R_SIZE * sizeof(int));
    int *local_a = (int *)malloc(block_size * R_SIZE * sizeof(int));
    int *local_c = (int *)malloc(block_size * R_SIZE * sizeof(int));

    if (rank == 0) {
        srand(time(NULL));
        initialize_matrices(a, b, R_SIZE * R_SIZE);
    }

    MPI_Bcast(b, R_SIZE * R_SIZE, MPI_INT, 0, MPI_COMM_WORLD);

    // 卷帘方式散列矩阵 a
    for (int i = 0; i < block_size; i++) {
        for (int j = 0; j < R_SIZE; j++) {
            local_a[i * R_SIZE + j] = a[(rank * block_size + i) * R_SIZE + j];
        }
    }

    for (int i = 0; i < block_size * R_SIZE; i++) {
        local_c[i] = 0;
    }

    clock_t start, end;
    if (rank == 0) {
        start = clock();
    }

    for (int i = 0; i < block_size; i++) {
        for (int j = 0; j < R_SIZE; j++) {
            for (int k = 0; k < R_SIZE; k++) {
                local_c[i * R_SIZE + j] += local_a[i * R_SIZE + k] * b[k * R_SIZE + j];
            }
        }
    }

    // 卷帘方式汇聚结果矩阵 c
    for (int i = 0; i < block_size; i++) {
        for (int j = 0; j < R_SIZE; j++) {
            c[(rank * block_size + i) * R_SIZE + j] = local_c[i * R_SIZE + j];
        }
    }

    if (rank == 0) {
        end = clock();
        printf("Time taken: %f\n", (double)(end - start) / CLOCKS_PER_SEC);
        // print_matrix(c, R_SIZE);
    }

    if (rank == 0) {
        free(a);
        free(b);
        free(c);
    }
    free(local_a);
    free(local_c);

    MPI_Finalize();
    return 0;
}

顺序划分

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>

#define M 1024  // 行数
#define N 512   // 列数
#define P 256   // 列数

void initialize_matrices(int *a, int *b, int rows_a, int cols_a, int cols_b) {
    for (int i = 0; i < rows_a * cols_a; i++) {
        a[i] = rand() % 10000;
    }
    for (int i = 0; i < cols_a * cols_b; i++) {
        b[i] = rand() % 10000;
    }
}

int main(int argc, char *argv[]) {
    int rank, size;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    int block_size = M / size;
    if (M % size != 0) {
        if (rank == 0) {
            printf("M must be divisible by size\n");
        }
        MPI_Finalize();
        return 1;
    }

    int *a = (int *)malloc(M * N * sizeof(int));
    int *b = (int *)malloc(N * P * sizeof(int));
    int *c = (int *)malloc(M * P * sizeof(int));
    int *local_a = (int *)malloc(block_size * N * sizeof(int));
    int *local_c = (int *)malloc(block_size * P * sizeof(int));

    if (rank == 0) {
        srand(time(NULL));
        initialize_matrices(a, b, M, N, P);
    }

    MPI_Bcast(b, N * P, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Scatter(a, block_size * N, MPI_INT, local_a, block_size * N, MPI_INT, 0, MPI_COMM_WORLD);

    for (int i = 0; i < block_size * P; i++) {
        local_c[i] = 0;
    }

    // 矩阵乘法
    for (int i = 0; i < block_size; i++) {
        for (int j = 0; j < P; j++) {
            for (int k = 0; k < N; k++) {
                local_c[i * P + j] += local_a[i * N + k] * b[k * P + j];
            }
        }
    }

    // 汇总结果
    MPI_Gather(local_c, block_size * P, MPI_INT, c, block_size * P, MPI_INT, 0, MPI_COMM_WORLD);

    if (rank == 0) {
        // 打印结果矩阵（仅适用于小矩阵）
        // for (int i = 0; i < M; i++) {
        //     for (int j = 0; j < P; j++) {
        //         printf("%d ", c[i * P + j]);
        //     }
        //     printf("\n");
        // }

        // 释放内存
        free(a);
        free(b);
        free(c);
    }

    free(local_a);
    free(local_c);

    MPI_Finalize();
    return 0;
}