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


typedef struct buffers_s
{
    double *buf1;
    double *buf2;
    int size;
} Buffers;



void BuffersMPI(Buffers *buf, MPI_Datatype *bufmpi)
{
    int block_lengths[2];                        // # of elt. in each block
    MPI_Aint displacements[2];                   // displac.
    MPI_Datatype typelist[2];                    // list of types
    MPI_Aint start_address, address;            // use for calculating displac.
    MPI_Datatype myType;

    block_lengths[0] = buf->size;
    block_lengths[1] = buf->size;

    typelist[0] = MPI_DOUBLE;
    typelist[1] = MPI_DOUBLE;

    displacements[0] = 0;
    MPI_Address(buf->buf1, &start_address);
    MPI_Address(buf->buf2, &address);
    displacements[1] = address - start_address;

    MPI_Type_struct(2,block_lengths, displacements,typelist,bufmpi);
    MPI_Type_commit(bufmpi);
}




void buffersmpisum(Buffers *in, Buffers *out, int *len, MPI_Datatype *typeptr)
{
    int i;

    for (i=0; i < *len; i++)
    {
        out->buf1[i] += in->buf1[i];
        out->buf2[i] += in->buf2[i];
    }
}




int main(int argc, char **argv)
{
    Buffers buffers_1, buffers_2;
    double *pack, *unpack;
    double size;
    double start, end;
    int pos;
    int i;
    int rank;
    int cnt;
    int ierr;
    MPI_Datatype Buffers_mpi;
    MPI_Op buffersumop;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    cnt = 100000;

    size = 10;

    buffers_1.buf1   = calloc(size, sizeof *buffers_1.buf1);
    buffers_1.buf2   = calloc(size, sizeof *buffers_1.buf2);
    buffers_2.buf1   = calloc(size, sizeof *buffers_2.buf1);
    buffers_2.buf2   = calloc(size, sizeof *buffers_2.buf2);
    pack   = calloc(3*size, sizeof *pack);
    unpack = calloc(3*size, sizeof *unpack);

    for (i=0; i<size; i++)
    {
        buffers_1.buf1[i] = 2;
        buffers_1.buf2[i] = 4;
    }

    MPI_Barrier(MPI_COMM_WORLD);

// ============================================================================
//                      TWO CONSECUTIVE REDUCTIONS
// ============================================================================


    start = MPI_Wtime();

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

        MPI_Reduce(buffers_1.buf1,               // send buffer
                   buffers_2.buf1,               // reduced data goes here
                   size,                          // size of the buffers
                   MPI_DOUBLE,                    // type of the data
                   MPI_SUM,                       // operation
                   0,                             // proc 0 only will recieve data
                   MPI_COMM_WORLD);               // all procs do the reduction


        MPI_Reduce(buffers_1.buf2,               // send buffer
                   buffers_2.buf2,               // reduced data goes here
                   size,                          // size of the buffers
                   MPI_DOUBLE,                    // type of the data
                   MPI_SUM,                       // operation
                   0,                             // proc 0 only will recieve data
                   MPI_COMM_WORLD);              // all procs do the reduction

    }
    end = MPI_Wtime();

    if(rank==0)
    {
        printf("%f seconds\n", end-start);
    }

    // empty bufffers_2
    memset(buffers_2.buf1, 0, size * sizeof *buffers_2.buf1);
    memset(buffers_2.buf2, 0, size * sizeof *buffers_2.buf2);



// ============================================================================
//                              PACK/UNPACK
// ============================================================================



    MPI_Barrier(MPI_COMM_WORLD);
    start = MPI_Wtime();
    for(i=0; i < cnt; i++)
    {
        pos = 0;                           // start to pack at pack+0

        ierr = MPI_Pack(buffers_1.buf1,              // data to be packed
                 size,                                // size of buf1
                 MPI_DOUBLE,                          // type of these elements
                 pack,                                // start adress of the pack
                 2*size*sizeof *pack,                 // size of the pack buffer
                 &pos,                                // pack+pos is the start position
                 MPI_COMM_WORLD);                     // every proc will use 'pack'



    // pack the second buffer

        MPI_Pack(buffers_1.buf2,                     // data to be packed
                 size,                                // size of buf1
                 MPI_DOUBLE,                          // type of these elements
                 pack,                                // start adress of the pack
                 2*size*sizeof *pack,                 // size of the pack buffer
                 &pos,                                // pack+pos is the start position
                 MPI_COMM_WORLD);                     // every proc will use 'pack'

        // now reduce the pack

        MPI_Reduce(pack,               // send buffer
                   unpack,             // reduced data goes here
                   2*size,             // size of the buffers
                   MPI_DOUBLE,         // type of the data
                   MPI_SUM,            // operation
                   0,                  // proc 0 only will recieve data
                   MPI_COMM_WORLD);    // all procs do the reduction

    }

    end = MPI_Wtime();
    if(rank==0)
    {
        printf("%f seconds\n", end-start);
    }



    if (rank == 0)
    {
        for (i=0; i<2*size; i++)
        {
            printf("unpack[%d] = %f\n",i,unpack[i]);
        }
    }


    // empty buffers_2
    memset(buffers_2.buf1, 0, size * sizeof *buffers_2.buf1);
    memset(buffers_2.buf2, 0, size * sizeof *buffers_2.buf2);


// ============================================================================
//                       MPI_Type_struct
// ============================================================================


    MPI_Barrier(MPI_COMM_WORLD);
    start = MPI_Wtime();

    BuffersMPI(&buffers_1, &Buffers_mpi);
    ierr = MPI_Op_create((MPI_User_function *)buffersmpisum,1,&buffersumop);

    if (ierr == MPI_SUCCESS) printf("op create success\n");

    for (i=0; i<cnt; i++)
    {
        MPI_Reduce(&buffers_1,               // send buffer
                   &buffers_2,               // reduced data goes here
                    size,                      // size of the buffers
                    Buffers_mpi,                  // type of the data
                    buffersumop,                     // operation
                    0,                           // proc 0 only will recieve data
                    MPI_COMM_WORLD);             // all procs do the reduction

    }

    end = MPI_Wtime();
    if(rank==0)
    {
        printf("%f seconds\n", end-start);
    }



    if (rank == 0)
    {
        for (i=0; i<  size; i++)
        {
            printf("buffers_2.buf1[%d] = %f \t buffers_2.buf2[%d] = %f\n",
                    i,buffers_2.buf1[i], i, buffers_2.buf2[i]);
        }
    }



    free(buffers_1.buf1);
    free(buffers_1.buf2);
    free(buffers_2.buf1);
    free(buffers_2.buf2);
    free(pack);
    free(unpack);


    MPI_Finalize();

    return 0;
}