apr8-parteplo.c

#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#include <math.h>
#include "l17.h"

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

    int size;
    rc = MPI_Comm_size(MPI_COMM_WORLD, &size);
    L17_CHECK_ERROR(MPI_Comm_size, rc);

    int rank;
    rc = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    L17_CHECK_ERROR(MPI_Comm_rank, rc);

    double *memu, *u, *memunew, *unew, delta, maxdelta, h, tau;
    double eps = 1.0e-6;
    int N, i;
    FILE *ff;

    if (argc != 2) {
        printf("Usage: exefile npoints\n");
        MPI_Abort(MPI_COMM_WORLD, 1);
    }

    N = atoi(argv[1]);
    if (N == 0) {
        printf("Set N to 1000\n");
        N = 1000;
    }
    else {
        printf("Set N to %d\n", N);
    }

    int jbeg, jend, sizej;
    jbeg = rank * ((N - 1) / size) + 1;
    jend = (rank == size - 1) ? (N - 1) : ((rank + 1) * ((N - 1) / size));
    sizej = jend - jbeg + 1;

    if ((memu = (double*) malloc((sizej + 2) * sizeof(double))) == NULL) {
        printf("Error. Couldn't allocate memory.\n");
        MPI_Abort(MPI_COMM_WORLD, 2);
    }

    if ((memunew = (double*) malloc((sizej + 2) * sizeof(double))) == NULL) {
        free(memu);
        printf("Error. Couldn't allocate memory.\n");
        MPI_Abort(MPI_COMM_WORLD, 3);
    }

    u = memu - jbeg + 1;
    unew = memunew - jbeg + 1;

    for (i = jbeg; i <= jend; ++i)
        u[i] = 0;

    if (rank == 0)
        unew[0] = u[0] = 1.0;

    if (rank == size - 1)
        unew[N] = u[N] = 0.0;

    h = 1.0 / (double) N;
    tau = 0.5 * h * h;

    while (1) {
        for (i = jbeg; i <= jend; ++i)
            unew[i] = u[i] + (tau / (h * h)) * (u[i - 1] - 2.0 * u[i] + u[i + 1]);

        maxdelta = 0.0;
        for (i = jbeg; i <= jend; ++i) {
            delta = fabs(unew[i] - u[i]);
            maxdelta = (delta > maxdelta) ? delta : maxdelta;
        }

        rc = MPI_Reduce((rank == 0) ? MPI_IN_PLACE : &maxdelta,
                        &maxdelta, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
        L17_CHECK_ERROR(MPI_Reduce, rc);

        rc = MPI_Bcast(&maxdelta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        L17_CHECK_ERROR(MPI_Reduce, rc);

        if (maxdelta < eps)
            break;

        for (i = jbeg; i <= jend; ++i)
            u[i] = unew[i];

        if (rank > 0) {
            rc = MPI_Send(u + jbeg, 1, MPI_DOUBLE, rank - 1, 1, MPI_COMM_WORLD);
            L17_CHECK_ERROR(MPI_Reduce, rc);
        }

        if (rank < size - 1) {
            rc = MPI_Recv(u + jend + 1, 1, MPI_DOUBLE, rank + 1, MPI_ANY_TAG, MPI_COMM_WORLD, NULL);
            L17_CHECK_ERROR(MPI_Reduce, rc);
        }

        if (rank < size - 1) {
            rc = MPI_Send(u + jend, 1, MPI_DOUBLE, rank + 1, 1, MPI_COMM_WORLD);
            L17_CHECK_ERROR(MPI_Reduce, rc);
        }

        if (rank > 0) {
            rc = MPI_Recv(u + jbeg - 1, 1, MPI_DOUBLE, rank - 1, MPI_ANY_TAG, MPI_COMM_WORLD, NULL);
            L17_CHECK_ERROR(MPI_Reduce, rc);
        }
    }

    free(memu);

    int* recvcounts;
    if (rank == 0) {
        if ((recvcounts = (int*) malloc(size * sizeof(int))) == NULL) {
            printf("Error. Couldn't allocate memory.\n");
            free(memunew);
            MPI_Abort(MPI_COMM_WORLD, 4);
        }
    }

    rc = MPI_Gather(&sizej, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, MPI_COMM_WORLD);
    L17_CHECK_ERROR(MPI_Gather, rc);

    int* displs;
    if (rank == 0) {
        if ((displs = (int*) malloc(size * sizeof(int))) == NULL) {
            printf("Error. Couldn't allocate memory.\n");
            free(memunew);
            free(recvcounts);
            MPI_Abort(MPI_COMM_WORLD, 5);
        }

        displs[0] = 1;
        for (i = 1; i < size; ++i)
            displs[i] = displs[i - 1] + recvcounts[i - 1];
    }

    double *recvbuf;
    if (rank == 0) {
        if ((recvbuf = (double*) malloc((N + 1) * sizeof(double))) == NULL) {
            printf("Error. Couldn't allocate memory.\n");
            free(memunew);
            free(recvcounts);
            free(displs);
            MPI_Abort(MPI_COMM_WORLD, 6);
        }
    }

    rc = MPI_Gatherv(unew + jbeg, sizej, MPI_DOUBLE, recvbuf, recvcounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    L17_CHECK_ERROR(MPI_Gatherv, rc);

    free(memunew);
    if (rank == 0) {
        recvbuf[0] = 1.0;
        recvbuf[N] = 0.0;

        free(recvcounts);
        free(displs);

        ff = fopen("apr8-parteplo.dat", "w+");
        if (ff == NULL) {
            printf("Couldn't open file.\n");
            free(recvbuf);
            MPI_Abort(MPI_COMM_WORLD, 7);
        }

        for (i = 0; i <= N; ++i)
            fprintf(ff, "%f ", recvbuf[i]);

        free(recvbuf);
        fclose(ff);
    }

    MPI_Finalize();

    return 0;
}