/****************************************************************************
* FILE: datatypes.c
* DESCRIPTION:
*  This code compares timings of MPI non-contiguous derived datatypes
*  against sending the same data without using derived datatypes.  The 
*  comparisons all send a vector of 1000 elements of type MPI_DOUBLE, with a 
*  stride of 24 between each element. The four tests use:
*    1. MPI_Type_vector 
*    2. MPI_Type_struct 
*    3. Hand packing and unpacking of the strided vector to/from a contiguous 
*       vector 
*    4. Individual send/recv for each strided element of the vector
*   
*  Cases 1 and 2 use MPI non-contiguous derived datatypes, cases 3 and 4 do
*  not.  Case 3 demonstrates how to improve on the MPI non-contiguous derived 
*  datatypes, while case 4 demonstrates a worst case method of sending
*  strided data elements.
*
* SOURCE: Adapted from example code provided with the "Performance Tuning
*   MPI" tutorial materials by William Gropp and Ewing Lusk, located at
*   www.mcs.anl.gov/mpi/tutorial 
* LAST REVISED: 01/09/99 Blaise Barney
******************************************************************************/

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

#define NUMBER_OF_TESTS 10

int main( argc, argv )
int argc;
char **argv;
{
    MPI_Datatype vec1, vec_n;
    int          blocklens[2];
    MPI_Aint     indices[2];
    MPI_Datatype old_types[2];

    double       *buf, *lbuf, t1, t2, tmin;
    int          rank, n, stride, i, j, k, nloop;
    MPI_Status   status;

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

    n	   = 1000;
    stride = 24;
    nloop  = 100000/n;

    buf = (double *) malloc( n * stride * sizeof(double) );
    if (!buf) {
	fprintf( stderr, "Could not allocate send/recv buffer of size %d\n",
		 n * stride );
	MPI_Abort( MPI_COMM_WORLD, 1 );
    }
    lbuf = (double *) malloc( n * sizeof(double) );
    if (!lbuf) {
	fprintf( stderr, "Could not allocate send/recv lbuffer of size %d\n",
		 n );
	MPI_Abort( MPI_COMM_WORLD, 1 );
    }

    if (rank == 0) 
	printf( "Kind\tn\tstride\ttime (sec)\tRate (MB/sec)\n" );

    /*************************** Test 1 ************************************/
    /* Use a fixed vector type */
    MPI_Type_vector( n, 1, stride, MPI_DOUBLE, &vec1 );
    MPI_Type_commit( &vec1 );

    tmin = 1000;
    for (k=0; k<NUMBER_OF_TESTS; k++) {
	if (rank == 0) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    t1 = MPI_Wtime();
	    for (j=0; j<nloop; j++) {
		MPI_Send( buf, 1, vec1, 1, k, MPI_COMM_WORLD );
		MPI_Recv( buf, 1, vec1, 1, k, MPI_COMM_WORLD, &status );
	    }
	    t2 = (MPI_Wtime() - t1) / nloop;
	    if (t2 < tmin) tmin = t2;
	}
	else if (rank == 1) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    for (j=0; j<nloop; j++) {
		MPI_Recv( buf, 1, vec1, 0, k, MPI_COMM_WORLD, &status );
		MPI_Send( buf, 1, vec1, 0, k, MPI_COMM_WORLD );
	    }
	}
    }
    /* Convert to half the round-trip time */
    tmin = tmin / 2.0;
    if (rank == 0) {
	printf( "Vector\t%d\t%d\t%f\t%f\n", 
		n, stride, tmin, n * sizeof(double) * 1.0e-6 / tmin );
    }
    MPI_Type_free( &vec1 );

    /*************************** Test 2 ************************************/
    /* Use a variable vector type */
    blocklens[0] = 1;
    blocklens[1] = 1;
    indices[0]   = 0;
    indices[1]   = stride * sizeof(double);
    old_types[0] = MPI_DOUBLE;
    old_types[1] = MPI_UB;
    MPI_Type_struct( 2, blocklens, indices, old_types, &vec_n );
    MPI_Type_commit( &vec_n );

    tmin = 1000;
    for (k=0; k<NUMBER_OF_TESTS; k++) {
	if (rank == 0) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    t1 = MPI_Wtime();
	    for (j=0; j<nloop; j++) {
		MPI_Send( buf, n, vec_n, 1, k, MPI_COMM_WORLD );
		MPI_Recv( buf, n, vec_n, 1, k, MPI_COMM_WORLD, &status );
	    }
	    t2 = (MPI_Wtime() - t1) / nloop;
	    if (t2 < tmin) tmin = t2;
	}
	else if (rank == 1) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    for (j=0; j<nloop; j++) {
		MPI_Recv( buf, n, vec_n, 0, k, MPI_COMM_WORLD, &status );
		MPI_Send( buf, n, vec_n, 0, k, MPI_COMM_WORLD );
	    }
	}
    }
    /* Convert to half the round-trip time */
    tmin = tmin / 2.0;
    if (rank == 0) {
	printf( "Struct\t%d\t%d\t%f\t%f\n", 
		n, stride, tmin, n * sizeof(double) * 1.0e-6 / tmin );
    }

    MPI_Type_free( &vec_n );

    /*************************** Test 3 ************************************/
    /* Use user-packing with known stride */
    tmin = 1000;
    for (k=0; k<NUMBER_OF_TESTS; k++) {
	if (rank == 0) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    t1 = MPI_Wtime();
	    for (j=0; j<nloop; j++) {
		for (i=0; i<n; i++) 
		    lbuf[i] = buf[i*stride];
		MPI_Send( lbuf, n, MPI_DOUBLE, 1, k, MPI_COMM_WORLD );
		MPI_Recv( lbuf, n, MPI_DOUBLE, 1, k, MPI_COMM_WORLD, &status );
		for (i=0; i<n; i++) 
		    buf[i*stride] = lbuf[i];
	    }
	    t2 = (MPI_Wtime() - t1) / nloop;
	    if (t2 < tmin) tmin = t2;
	}
	else if (rank == 1) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    for (j=0; j<nloop; j++) {
		MPI_Recv( lbuf, n, MPI_DOUBLE, 0, k, MPI_COMM_WORLD, &status );
		for (i=0; i<n; i++) 
		    buf[i*stride] = lbuf[i];
		for (i=0; i<n; i++) 
		    lbuf[i] = buf[i*stride];
		MPI_Send( lbuf, n, MPI_DOUBLE, 0, k, MPI_COMM_WORLD );
	    }
	}
    }
    /* Convert to half the round-trip time */
    tmin = tmin / 2.0;
    if (rank == 0) {
	printf( "User\t%d\t%d\t%f\t%f\n", 
		n, stride, tmin, n * sizeof(double) * 1.0e-6 / tmin );
    }

    /*************************** Test 4 ************************************/
    /* Worst case: use individual send/recv for each element sent */
    tmin = 1000;
    for (k=0; k<NUMBER_OF_TESTS; k++) {
	if (rank == 0) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    t1 = MPI_Wtime();
	    for (j=0; j<nloop; j++) {
              for (i=0; i<n; i++) {
                MPI_Send(&buf[i*stride], 1, MPI_DOUBLE, 1, k, MPI_COMM_WORLD);
                MPI_Recv(&buf[i*stride], 1, MPI_DOUBLE, 1, k, MPI_COMM_WORLD, 
                         &status);
              }
	    }
	    t2 = (MPI_Wtime() - t1) / nloop;
	    if (t2 < tmin) tmin = t2;
	}
	else if (rank == 1) {
	    /* Make sure both processes are ready */
            MPI_Barrier(MPI_COMM_WORLD);

	    for (j=0; j<nloop; j++) {
              for (i=0; i<n; i++) {
		MPI_Recv(&buf[i*stride], 1, MPI_DOUBLE, 0, k, MPI_COMM_WORLD, 
                         &status);
		MPI_Send(&buf[i*stride], 1, MPI_DOUBLE, 0, k, MPI_COMM_WORLD);
                }
	    }
	}
    }
    /* Convert to half the round-trip time */
    tmin = tmin / 2.0;
    if (rank == 0) {
	printf( "Worst case\t%d\t%d\t%f\t%f\n", 
		n, stride, tmin, n * sizeof(double) * 1.0e-6 / tmin );
    }

    MPI_Finalize( );
    return 0;
}