DPHPC Homework 1

DPHPC-HW1.c

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

#define SIZEX 2000
#define SIZEY 2000

#define IDX(a, x, y) a[(x) + (y) * (size_x + 2)]

enum tag {
    RIGHT,
    DOWN,
    LEFT,
    UP,
};

int main(int argc, char **argv) {

    // Initialize MPI
    MPI_Init(&argc, &argv);

    // Get the number of processes and rank
    int size, rank;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    // Calculate layout of blocking
    // ASSUMPTION: the number of processes is square
    int dim, pos_x, pos_y, size_x, size_y;
    dim = (int) sqrt(size);
    pos_x = rank % dim;
    pos_y = rank / dim;
    size_x = SIZEX / dim;
    size_y = SIZEY / dim;

    // Allocate memory
    double *A = (double *) calloc((size_x + 2) * (size_y + 2), sizeof(double));
    double *B = (double *) calloc((size_x + 2) * (size_y + 2), sizeof(double));
    assert(A);
    assert(B);

    // Set initial value
    int init_x = (SIZEX / 2) - (pos_x * size_x);
    int init_y = (SIZEY / 2) - (pos_y * size_y);
    if (init_x >= 0 && init_x < size_x && init_y >= 0 && init_y < size_y) {
        IDX(A, init_x + 1, init_y + 1) = 10000.0;
    }

    // Create MPI column datatype
    MPI_Datatype column_type;
    MPI_Type_vector(size_y, 1, size_x + 2, MPI_DOUBLE, &column_type);
    MPI_Type_commit(&column_type);

    // Iterate over steps
    for (int steps = 0; steps < 500; steps++) {

        // Send last column to the right
        if (pos_x < dim - 1) {
            MPI_Send(&IDX(A, size_x, 1), 1, column_type, rank + 1, RIGHT,
                     MPI_COMM_WORLD);
        }

        // Receive last column from the left
        if (pos_x > 0) {
            MPI_Recv(&IDX(A, 0, 1), 1, column_type, rank - 1, RIGHT, MPI_COMM_WORLD, NULL);
        }

        // Send last row to below
        if (pos_y < dim - 1) {
            MPI_Send(&IDX(A, 1, size_y), size_x, MPI_DOUBLE, rank + dim, DOWN, MPI_COMM_WORLD);
        }

        // Receive last row from above
        if (pos_y > 0) {
            MPI_Recv(&IDX(A, 1, 0), size_x, MPI_DOUBLE, rank - dim, DOWN, MPI_COMM_WORLD, NULL);
        }

        // Send first column to the left
        if (pos_x > 0) {
            MPI_Send(&IDX(A, 1, 1), 1, column_type, rank - 1, LEFT, MPI_COMM_WORLD);
        }

        // Receive first column from the right
        if (pos_x < dim - 1) {
            MPI_Recv(&IDX(A, size_x + 1, 1), 1, column_type, rank + 1, LEFT, MPI_COMM_WORLD, NULL);
        }

        // Send first row to above
        if (pos_y > 0) {
            MPI_Send(&IDX(A, 1, 1), size_x, MPI_DOUBLE, rank - dim, UP, MPI_COMM_WORLD);
        }

        // Receive first row from below
        if (pos_y < dim - 1) {
            MPI_Recv(&IDX(A, 1, size_y + 1), size_x, MPI_DOUBLE, rank + dim, UP, MPI_COMM_WORLD, NULL);
        }

        // Execute next step
        for (int i = 1; i <= size_x; i++) {
            for (int j = 1; j <= size_y; j++) {
                IDX(B, i, j) =
                        (IDX(A, i, j) + IDX(A, i + 1, j) + IDX(A, i - 1, j) +
                         IDX(A, i, j + 1) + IDX(A, i, j - 1)) / 5.0;
            }
        }

        // Switch matrices
        double *temp = B;
        B = A;
        A = temp;

    }

    // Sum all values and reduce
    double sum = 0.0, global_sum;
    for (int i = 1; i <= size_x; i++) {
        for (int j = 1; j <= size_y; j++) {
            sum += IDX(A, i, j);
        }
    }
    MPI_Reduce(&sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

    // Print the result at the root
    if (!rank) {
        printf("sum = %0.20lf\n", global_sum);
    }

    // Finalize MPI
    MPI_Finalize();

    return 0;
}

Makefile

.PHONY: run all

stencil: DPHPC-HW1.c
	mpicc -o $@ $^

run: stencil
	mpirun -np 4 $<

all: stencil