/*
 * The program presented here can be compiled for serial or parallel
 * execution using h5cc or h5pcc, respectively. The program receives
 * arguments from the command line as follows:
 * 
 * Serial execution:
 * 
 * a.out <t|h>
 * 
 * where t and h indicate contiguous and chunked storage, correspondingly.
 * 
 * Parallel execution:
 * 
 * mpirun .np num-processes a.out <t|h> <i|c>
 * 
 * where i and c specify independent and collective access, respectively. For
 * a proper parallel execution, the number of processes must be a power
 * of two.
 * 
 * Scalable parallel filesystems like GPFS can yield parallel speedups
 * close to the number of used processes. This can be verified by
 * comparing serial and parallel performance using contiguous storage and
 * independent access. Filesystems like Lustre and CXFS can still achieve
 * high performance by using relatively large stripe sizes (>1 MB) and
 * optimizations like collective access.
 */

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "hdf5.h"
#define NX    8192
#define NY    8192       /* dataset dimensions */
#define HY    256        /* width of vertical hyperslab */
int main (int argc, char **argv) {
    hid_t   file_id, dset_id, filespace, memspace, fapl, dxpl, dcpl;
    hsize_t dimsf[2], count[2], offset[2], chunk_dims[2]={256, 256};
    int     *data;
    unsigned long i, group;
    int         mpi_size, mpi_rank;
    fapl = H5Pcreate(H5P_FILE_ACCESS);
    dcpl = H5Pcreate(H5P_DATASET_CREATE);
    dxpl = H5Pcreate(H5P_DATASET_XFER);
    if (argc != 3) {
	fprintf(stderr, "Usage: %s <t|h> <i|c>\n"
	    "where\tt and h are for contiguous and chunked storage\n"
	    "\ti and c are for independent and collective access\n",
	    argv[0]
	    );
	exit(1);
    }
    if (!strcmp(argv[1], "h"))
        H5Pset_chunk ( dcpl, 2, chunk_dims);
    else
        dcpl = H5P_DEFAULT;
#ifdef H5_HAVE_PARALLEL
    MPI_Comm comm  = MPI_COMM_WORLD;
    MPI_Info info  = MPI_INFO_NULL;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(comm, &mpi_size);
    MPI_Comm_rank(comm, &mpi_rank);
    H5Pset_fapl_mpio(fapl, comm, info);
    if (!strcmp(argv[2], "c"))
        H5Pset_dxpl_mpio(dxpl, H5FD_MPIO_COLLECTIVE);
    else
        dxpl = H5P_DEFAULT;
#else
    mpi_size = 1;
    mpi_rank = 0;
    fapl = dxpl = H5P_DEFAULT;
#endif
    file_id = H5Fcreate("file1.h5", H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
    dimsf[0] = NX;
    dimsf[1] = NY;
    filespace = H5Screate_simple(2, dimsf, NULL);
    dset_id = H5Dcreate(file_id, "dataset1", H5T_NATIVE_INT, filespace,
                        H5P_DEFAULT, dcpl, H5P_DEFAULT);
    group  = dimsf[1]/(mpi_size*HY);
    count[0] = dimsf[0];
    count[1] = HY;
    memspace = H5Screate_simple(2, count, NULL);
    data = (int *)malloc(count[0]*count[1]*sizeof(int));
    for (i=0; i < count[0]*count[1]; i++){
        data[i] = mpi_rank + 10;
    }
    offset[0] = 0;
    for (i=0; i<group; i++) {
        offset[1] = (mpi_size*HY)*i + mpi_rank * HY;
        H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
        H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace, dxpl, data);
    }
    free(data);
    H5Dclose(dset_id);
    H5Sclose(filespace);
    H5Sclose(memspace);
    H5Pclose(dxpl);
    H5Pclose(dcpl);
    H5Pclose(fapl);
    H5Fclose(file_id);
#ifdef H5_HAVE_PARALLEL
    MPI_Finalize();
#endif
    return 0;
}