/* -*- Mode: C; c-basic-offset:4 ; indent-tabs-mode:nil ; -*- */

/*

 *  Changes to this example

 *  (C) 2001 by Argonne National Laboratory.

 *      See COPYRIGHT in top-level directory.

 */

/*

 * This example is taken from MPI-The complete reference, Vol 1,

 * pages 222-224.

 *

 * Lines after the "--CUT HERE--" were added to make this into a complete

 * test program.

 */

/* Specify the maximum number of errors to report. */

#define MAX_ERRORS 10

#include "mpi.h"

#include "mpitest.h"

#include <stdio.h>

#include <stdlib.h>

#define MAX_SIZE 64

MPI_Datatype transpose_type(int M, int m, int n, MPI_Datatype type);

MPI_Datatype submatrix_type(int N, int m, int n, MPI_Datatype type);

void Transpose(float *localA, float *localB, int M, int N, MPI_Comm comm);

void Transpose(float *localA, float *localB, int M, int N, MPI_Comm comm)

/* transpose MxN matrix A that is block distributed (1-D) on

   processes of comm onto block distributed matrix B  */

{

    int i, j, extent, myrank, p, n[2], m[2];

    int lasti, lastj;

    int *sendcounts, *recvcounts;

    int *sdispls, *rdispls;

    MPI_Datatype xtype[2][2], stype[2][2], *sendtypes, *recvtypes;

    MTestPrintfMsg(2, "M = %d, N = %d\n", M, N);

    /* compute parameters */

    MPI_Comm_size(comm, &p);

    MPI_Comm_rank(comm, &myrank);

    extent = sizeof(float);

    /* allocate arrays */

    sendcounts = (int *) malloc(p * sizeof(int));

    recvcounts = (int *) malloc(p * sizeof(int));

    sdispls = (int *) malloc(p * sizeof(int));

    rdispls = (int *) malloc(p * sizeof(int));

    sendtypes = (MPI_Datatype *) malloc(p * sizeof(MPI_Datatype));

    recvtypes = (MPI_Datatype *) malloc(p * sizeof(MPI_Datatype));

    /* compute block sizes */

    m[0] = M / p;

    m[1] = M - (p - 1) * (M / p);

    n[0] = N / p;

    n[1] = N - (p - 1) * (N / p);

    /* compute types */

    for (i = 0; i <= 1; i++)

        for (j = 0; j <= 1; j++) {

            xtype[i][j] = transpose_type(N, m[i], n[j], MPI_FLOAT);

            stype[i][j] = submatrix_type(M, m[i], n[j], MPI_FLOAT);

        }

    /* prepare collective operation arguments */

    lasti = myrank == p - 1;

    for (j = 0; j < p; j++) {

        lastj = j == p - 1;

        sendcounts[j] = 1;

        sdispls[j] = j * n[0] * extent;

        sendtypes[j] = xtype[lasti][lastj];

        recvcounts[j] = 1;

        rdispls[j] = j * m[0] * extent;

        recvtypes[j] = stype[lastj][lasti];

    }

    /* communicate */

    MTestPrintfMsg(2, "Begin Alltoallw...\n");

    /* -- Note that the book incorrectly uses &localA and &localB

     * as arguments to MPI_Alltoallw */

    MPI_Alltoallw(localA, sendcounts, sdispls, sendtypes,

                  localB, recvcounts, rdispls, recvtypes, comm);

    MTestPrintfMsg(2, "Done with Alltoallw\n");

    /* Free buffers */

    free(sendcounts);

    free(recvcounts);

    free(sdispls);

    free(rdispls);

    free(sendtypes);

    free(recvtypes);

    /* Free datatypes */

    for (i = 0; i <= 1; i++)

        for (j = 0; j <= 1; j++) {

            MPI_Type_free(&xtype[i][j]);

            MPI_Type_free(&stype[i][j]);

        }

}

/* Define an n x m submatrix in a n x M local matrix (this is the

   destination in the transpose matrix */

MPI_Datatype submatrix_type(int M, int m, int n, MPI_Datatype type)

/* computes a datatype for an mxn submatrix within an MxN matrix

   with entries of type type */

{

    /* MPI_Datatype subrow; */

    MPI_Datatype submatrix;

    /* The book, MPI: The Complete Reference, has the wrong type constructor

     * here.  Since the stride in the vector type is relative to the input

     * type, the stride in the book's code is n times as long as is intended.

     * Since n may not exactly divide N, it is better to simply use the

     * blocklength argument in Type_vector */

    /*

     * MPI_Type_contiguous(n, type, &subrow);

     * MPI_Type_vector(m, 1, N, subrow, &submatrix);

     */

    MPI_Type_vector(n, m, M, type, &submatrix);

    MPI_Type_commit(&submatrix);

    /* Add a consistency test: the size of submatrix should be

     * n * m * sizeof(type) and the extent should be ((n-1)*M+m) * sizeof(type) */

    {

        int tsize;

        MPI_Aint textent, lb;

        MPI_Type_size(type, &tsize);

        MPI_Type_get_extent(submatrix, &lb, &textent);

        if (textent != tsize * (M * (n - 1) + m)) {

            fprintf(stderr, "Submatrix extent is %ld, expected %ld (%d,%d,%d)\n",

                    (long) textent, (long) (tsize * (M * (n - 1) + m)), M, n, m);

        }

    }

    return (submatrix);

}

/* Extract an m x n submatrix within an m x N matrix and transpose it.

   Assume storage by rows; the defined datatype accesses by columns */

MPI_Datatype transpose_type(int N, int m, int n, MPI_Datatype type)

/* computes a datatype for the transpose of an mxn matrix

   with entries of type type */

{

    MPI_Datatype subrow, subrow1, submatrix;

    MPI_Aint lb, extent;

    MPI_Type_vector(m, 1, N, type, &subrow);

    MPI_Type_get_extent(type, &lb, &extent);

    MPI_Type_create_resized(subrow, 0, extent, &subrow1);

    MPI_Type_contiguous(n, subrow1, &submatrix);

    MPI_Type_commit(&submatrix);

    MPI_Type_free(&subrow);

    MPI_Type_free(&subrow1);

    /* Add a consistency test: the size of submatrix should be

     * n * m * sizeof(type) and the extent should be ((m-1)*N+n) * sizeof(type) */

    {

        int tsize;

        MPI_Aint textent, llb;

        MPI_Type_size(type, &tsize);

        MPI_Type_get_true_extent(submatrix, &llb, &textent);

        if (textent != tsize * (N * (m - 1) + n)) {

            fprintf(stderr, "Transpose Submatrix extent is %ld, expected %ld (%d,%d,%d)\n",

                    (long) textent, (long) (tsize * (N * (m - 1) + n)), N, n, m);

        }

    }

    return (submatrix);

}

/* -- CUT HERE -- */

int main(int argc, char *argv[])

{

    int gM, gN, lm, lmlast, ln, lnlast, i, j, errs = 0;

    int size, rank;

    float *localA, *localB;

    MPI_Comm comm;

    MTest_Init(&argc, &argv);

    comm = MPI_COMM_WORLD;

    MPI_Comm_size(comm, &size);

    MPI_Comm_rank(comm, &rank;;

    gM = 20;

    gN = 30;

    /* Each block is lm x ln in size, except for the last process,

     * which has lmlast x lnlast */

    lm = gM / size;

    lmlast = gM - (size - 1) * lm;

    ln = gN / size;

    lnlast = gN - (size - 1) * ln;

    /* Create the local matrices.

     * Initialize the input matrix so that the entries are

     * consequtive integers, by row, starting at 0.

     */

    if (rank == size - 1) {

        localA = (float *) malloc(gN * lmlast * sizeof(float));

        localB = (float *) malloc(gM * lnlast * sizeof(float));

        for (i = 0; i < lmlast; i++) {

            for (j = 0; j < gN; j++) {

                localA[i * gN + j] = (float) (i * gN + j + rank * gN * lm);

            }

        }

    } else {

        localA = (float *) malloc(gN * lm * sizeof(float));

        localB = (float *) malloc(gM * ln * sizeof(float));

        for (i = 0; i < lm; i++) {

            for (j = 0; j < gN; j++) {

                localA[i * gN + j] = (float) (i * gN + j + rank * gN * lm);

            }

        }

    }

    MTestPrintfMsg(2, "Allocated local arrays\n");

    /* Transpose */

    Transpose(localA, localB, gM, gN, comm);

    /* check the transposed matrix

     * In the global matrix, the transpose has consequtive integers,

     * organized by columns.

     */

    if (rank == size - 1) {

        for (i = 0; i < lnlast; i++) {

            for (j = 0; j < gM; j++) {

                int expected = i + gN * j + rank * ln;

                if ((int) localB[i * gM + j] != expected) {

                    if (errs < MAX_ERRORS)

                        printf("Found %d but expected %d\n", (int) localB[i * gM + j], expected);

                    errs++;

                }

            }

        }

    } else {

        for (i = 0; i < ln; i++) {

            for (j = 0; j < gM; j++) {

                int expected = i + gN * j + rank * ln;

                if ((int) localB[i * gM + j] != expected) {

                    if (errs < MAX_ERRORS)

                        printf("Found %d but expected %d\n", (int) localB[i * gM + j], expected);

                    errs++;

                }

            }

        }

    }

    /* Free storage */

    free(localA);

    free(localB);

    MTest_Finalize(errs);

    return MTestReturnValue(errs);

}