/*
   Copyright (C) 2013, Northwestern University

   Please email questions to Wei-keng Liao <wkliao@eecs.northwestern.edu>

   This example program illustrates the interoperability of PnetCDF and netCDF.
   It uses netCDF-4 APIs to read 2D arrays in parallel from a shared file in the
   classical netCDF format. The global 2D array is partitioned among processes
   in a block-block fashion along both X and Y dimensions. To enable PnetCDF to
   carry out parallel I/O underneath, add the NC_PNETCDF option to the file open
   mode when opening the file. Compile and run commands are given below, along
   with the results of running command "ncdump" to the input file.

   1. At the time of this program was developed, the following libraries were
      used.
            HDF5    version 1.8.10
            netCDF  version 4.4.0 and
            PnetCDF version 1.4.0

   2. To build netCDF, the following configure options were used.
        ./configure --prefix=/usr/local \
                    --disable-shared \
                    --enable-netcdf-4 \
                    --enable-pnetcdf \
                    FC=mpif90 CXX=mpicxx CC=mpicc

   3. To compile and link this example program:
      mpicc -O2 -o nc4_pnc_get.c ./nc4_pnc_get \
           -I/path/PnetCDF/include -I/path/netCDF/include -I/path/HDF5/include \
           -L/path/PnetCDF/lib -L/path/netCDF/lib -L/path/HDF5/lib \
           -lnetcdf -lhdf5_hl -lhdf5 -lpnetcdf -lz -lcurl -ldl -lm

   4. Input file:
      Run this program using the output file generated from nc4_pnc_put.c.
      The output file from nc4_pnc_put.c is in CDF-1 format. The screenshot
      of running command "ncdump" to the output file is given below.
      % ncdump testfile.nc
      netcdf testfile {
      dimensions:
                Y = 10 ;
                X = 10 ;
      variables:
                int var(Y, X) ;
      data:

       var =
        10000, 10001, 10002, 10003, 10004, 10100, 10101, 10102, 10103, 10104,
        10005, 10006, 10007, 10008, 10009, 10105, 10106, 10107, 10108, 10109,
        10010, 10011, 10012, 10013, 10014, 10110, 10111, 10112, 10113, 10114,
        10015, 10016, 10017, 10018, 10019, 10115, 10116, 10117, 10118, 10119,
        10020, 10021, 10022, 10023, 10024, 10120, 10121, 10122, 10123, 10124,
        10200, 10201, 10202, 10203, 10204, 10300, 10301, 10302, 10303, 10304,
        10205, 10206, 10207, 10208, 10209, 10305, 10306, 10307, 10308, 10309,
        10210, 10211, 10212, 10213, 10214, 10310, 10311, 10312, 10313, 10314,
        10215, 10216, 10217, 10218, 10219, 10315, 10316, 10317, 10318, 10319,
        10220, 10221, 10222, 10223, 10224, 10320, 10321, 10322, 10323, 10324 ;
      }

   5. run command (of running on 4 processes):
      mpiexec -machinefile hostfile -n 4 nc4_pnc_get testfile.nc

      rank  0 get start= 0  0 count= 5  5
      rank  1 get start= 0  5 count= 5  5
      rank  2 get start= 5  0 count= 5  5
      rank  3 get start= 5  5 count= 5  5
      P0 [0, 0...4]: 10000 10001 10002 10003 10004 
         [1, 0...4]: 10005 10006 10007 10008 10009 
         [2, 0...4]: 10010 10011 10012 10013 10014 
         [3, 0...4]: 10015 10016 10017 10018 10019 
         [4, 0...4]: 10020 10021 10022 10023 10024 
      P1 [0, 0...4]: 10100 10101 10102 10103 10104 
         [1, 0...4]: 10105 10106 10107 10108 10109 
         [2, 0...4]: 10110 10111 10112 10113 10114 
         [3, 0...4]: 10115 10116 10117 10118 10119 
         [4, 0...4]: 10120 10121 10122 10123 10124 
      P2 [0, 0...4]: 10200 10201 10202 10203 10204 
         [1, 0...4]: 10205 10206 10207 10208 10209 
         [2, 0...4]: 10210 10211 10212 10213 10214 
         [3, 0...4]: 10215 10216 10217 10218 10219 
         [4, 0...4]: 10220 10221 10222 10223 10224 
      P3 [0, 0...4]: 10300 10301 10302 10303 10304 
         [1, 0...4]: 10305 10306 10307 10308 10309 
         [2, 0...4]: 10310 10311 10312 10313 10314 
         [3, 0...4]: 10315 10316 10317 10318 10319 
         [4, 0...4]: 10320 10321 10322 10323 10324 
*/

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

#define FATAL_ERR {if(err!=NC_NOERR) {printf("Error at line=%d: %s Aborting ...\n", __LINE__, nc_strerror(err)); goto fn_exit;}}
#define ERR {if(err!=NC_NOERR)printf("Error at line=%d: %s\n", __LINE__, nc_strerror(err));}

int main(int argc, char* argv[])
{
    int i, j, rank, nprocs, err, ncid, omode, varid, mpi_namelen, verbose=1;
    int dimid[2], psizes[2], *buf, format;
    size_t start[2], count[2], gsizes[2], lsizes[2];
    MPI_Info info=MPI_INFO_NULL;
    char mpi_name[MPI_MAX_PROCESSOR_NAME], *stdout_buf, *ptr;

    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Get_processor_name(mpi_name,&mpi_namelen);

    if (argc != 2) {
        if (!rank) printf("Usage: %s filename\n",argv[0]);
        goto fn_exit;
    }
    if (verbose) printf("rank %2d runs on host %s\n",rank,mpi_name);

    /* open file and indicate the parallel I/O method:
     * Depending on the file format, netCDF picks automatically PnetCDF or
     * HDF5 to perform parallel I/O. No additional flag is needed.
     * Using NC_MPIIO is no longer required, as it has been deprecated since
     * NetCDF 4.6.2
     */
    omode = NC_NOWRITE | NC_MPIIO;
    err = nc_open_par(argv[1], omode, MPI_COMM_WORLD, info, &ncid); FATAL_ERR

    /* free info object */
    if (info != MPI_INFO_NULL) MPI_Info_free(&info);

    /* create a block-block 2D data partitioning pattern */
    psizes[0] = psizes[1] = 0;
    MPI_Dims_create(nprocs, 2, psizes);

    /* inquire dimension IDs and lengths */
    err = nc_inq_dimid(ncid, "Y", &dimid[0]); ERR
    err = nc_inq_dimid(ncid, "X", &dimid[1]); ERR

    err = nc_inq_dimlen(ncid, dimid[0], &gsizes[0]); ERR
    err = nc_inq_dimlen(ncid, dimid[1], &gsizes[1]); ERR
    if (verbose)
        printf("global dimension Y is of size %zd, X is of size %zd\n", gsizes[0], gsizes[1]);

    /* obtain variable ID */
    err = nc_inq_varid(ncid, "var", &varid); ERR

    /* set the access method to use MPI collective I/O */
    err = nc_inq_format(ncid, &format); ERR
    if (format == NC_FORMAT_CLASSIC || format == NC_FORMAT_64BIT_OFFSET ||
        format == NC_FORMAT_CDF5) {
        err = nc_var_par_access(ncid, NC_GLOBAL, NC_COLLECTIVE); ERR
    } else {
        err = nc_var_par_access(ncid, varid, NC_COLLECTIVE); ERR
    }

    /* allocate space for the local 2D subarray */
    lsizes[0] = gsizes[0] / psizes[0];
    lsizes[1] = gsizes[1] / psizes[1];
    buf = (int*) calloc(lsizes[0] * lsizes[1], sizeof(int));
    stdout_buf = (char*) malloc(lsizes[0] * lsizes[1] * 16);

    /* each process reads its subarray from the file */
    start[0] = lsizes[0] * (rank / psizes[1]);
    start[1] = lsizes[1] * (rank % psizes[1]);
    count[0] = lsizes[0];
    count[1] = lsizes[1];

    /* parallel read: each process reads a subarray from the file */
    err = nc_get_vara_int(ncid, varid, start, count, buf); ERR
    if (verbose)
        printf("rank %2d get start=%2zd %2zd count=%2zd %2zd\n",
               rank,start[0],start[1],count[0],count[1]);

    /* close file */
    err = nc_close(ncid); ERR

    /* check and valid the contents of read values */
    ptr = stdout_buf;
    for (j=0; j<lsizes[0]; j++) {
        if (j == 0) sprintf(ptr, "P%d [%d, 0...%zd]: ", rank, j, lsizes[1]-1);
        else        sprintf(ptr,  "   [%d, 0...%zd]: ",       j, lsizes[1]-1);
        ptr += strlen(ptr);
        for (i=0; i<lsizes[1]; i++) {
            sprintf(ptr, "%2d ",buf[j*lsizes[1]+i]);
            ptr += strlen(ptr);
        }
        sprintf(ptr,"\n");
        ptr++;
    }
    printf("%s",stdout_buf);
    free(stdout_buf);
    free(buf);

fn_exit:
    MPI_Finalize();
    return 0;
}