llc_90/code-async-noseaice/readtile_mpiio.c

//
// MPI IO for MITgcm
//

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

#include "SIZE.h"


// lat-lon-cap decomposition has 13 square facets
// facetElements1D is typically 1080, or 2160, or 4320


//////////////////////////////////////////////////////////////
// These values filled in during "initSizesAndTypes()"
MPI_Datatype  fieldElementalTypeMPI;
size_t  sizeofFieldElementalType;

long int tileSizeX;
long int tileSizeY;
long int xGhosts;
long int yGhosts;

long int facetElements1D;
long int facetBytes2D;
long int facetTilesInX;
long int facetTilesInY;
long int tilesPerFacet;

long int fieldZlevelSizeInBytes;
long int tileZlevelSizeInBytes;
long int ghostedTileZlevelSizeInBytes;

// The first 7 facets all use the same style of layout; here we call
// them "section1".  The last 6 facets share a layout style, but this
// style is different than section1.  Here, we call them "section2".
MPI_Datatype  section1_ioShape2D, section2_ioShape2D;
MPI_Datatype  tileShape2D, ghostedTileShape2D;
MPI_Info  ioHints;
//////////////////////////////////////////////////////////////


int
getSizeOfMPIType(MPI_Datatype mpi_type)
{
    switch (mpi_type) {

      case MPI_INT:  case MPI_FLOAT:  case MPI_REAL4:
        return 4;
      break;

      case MPI_LONG_INT:  case MPI_DOUBLE:  case MPI_REAL8:
        return 8;
      break;

      default:
        assert(("unexpected mpi elemental type", 0));
      break;

    }
    return -1;
}


void
createMPItypes(void)
{
    // Create a type with the "shape" of a section1, 2D tile
    MPI_Type_vector(tileSizeY, tileSizeX, facetElements1D,
                    fieldElementalTypeMPI, &section1_ioShape2D);
    MPI_Type_commit(&section1_ioShape2D);

    // Create a type with the "shape" of a section2, 2D tile
    MPI_Type_vector(tileSizeY, tileSizeX, 3*facetElements1D,
                    fieldElementalTypeMPI, &section2_ioShape2D);
    MPI_Type_commit(&section2_ioShape2D);


    // Create a type that describes a 2D tile in memory
    MPI_Type_vector(tileSizeY, tileSizeX, tileSizeX,
                    fieldElementalTypeMPI, &tileShape2D);
    MPI_Type_commit(&tileShape2D);

    // Create a type that describes a 2D tile in memory with ghost-cells.
    int fullDims[2] = {tileSizeX + 2*xGhosts, tileSizeY + 2*yGhosts};
    int tileDims[2] = {tileSizeX, tileSizeY};
    int startElements[2] = {xGhosts, yGhosts};
    MPI_Type_create_subarray(2, fullDims, tileDims, startElements,
              MPI_ORDER_FORTRAN, fieldElementalTypeMPI, &ghostedTileShape2D);
    MPI_Type_commit(&ghostedTileShape2D);


    // Set up some possible hints
    MPI_Info_create(&ioHints);
    MPI_Info_set(ioHints, "collective_buffering", "true");
    char blockSize[64];
    sprintf(blockSize, "%ld", (((long)facetElements1D * 3) *
            tileSizeY) * sizeofFieldElementalType);
    MPI_Info_set(ioHints, "cb_block_size", blockSize);
}


void
initSizesAndTypes(void)
{
    /////////////////////////////////////////////
    // Fundamental values
    fieldElementalTypeMPI = MPI_DOUBLE;
    facetElements1D = sFacet;
    tileSizeX = sNx;
    tileSizeY = sNy;
    xGhosts = OLx;
    yGhosts = OLy;
    /////////////////////////////////////////////

    // Derived values
    sizeofFieldElementalType = getSizeOfMPIType(fieldElementalTypeMPI);
    long int facetElements2D = ((facetElements1D) * (facetElements1D));
    facetBytes2D  = (facetElements2D * sizeofFieldElementalType);
    fieldZlevelSizeInBytes = (13*(facetBytes2D));

    facetTilesInX = ((facetElements1D)/(tileSizeX));
    facetTilesInY = ((facetElements1D)/(tileSizeY));
    tilesPerFacet = ((facetTilesInX)*(facetTilesInY));

    tileZlevelSizeInBytes = tileSizeX * tileSizeY * sizeofFieldElementalType;
    ghostedTileZlevelSizeInBytes = (tileSizeX + 2*xGhosts) *
                                   (tileSizeY + 2*yGhosts) *
                                   sizeofFieldElementalType;

    // Create the specialized type definitions
    createMPItypes();
}


void
tileIO(
  MPI_Comm comm,
  char *filename,
  MPI_Offset tileOffsetInFile,
  MPI_Datatype tileLayoutInFile,
  void *tileBuf,
  MPI_Datatype tileLayoutInMemory,
  int  writeFlag)
{
    int fileFlags;
    MPI_File fh;
    int (*MPI_IO)();

    int res,count;
    MPI_Status status;

    if (writeFlag) {
        fileFlags = MPI_MODE_WRONLY | MPI_MODE_CREATE;
        MPI_IO = MPI_File_write_all;
    } else {
        fileFlags = MPI_MODE_RDONLY;
        MPI_IO = MPI_File_read_all;
    }

    //printf("filename is %s\n",filename);

    MPI_File_open(comm, filename, fileFlags, ioHints, &fh);
    MPI_File_set_view(fh, tileOffsetInFile, fieldElementalTypeMPI,
                      tileLayoutInFile, "native", ioHints);


    // MPI_IO(fh, tileBuf, 1, tileLayoutInMemory, MPI_STATUS_IGNORE);
    res = MPI_IO(fh, tileBuf, 1, tileLayoutInMemory, &status);
    
    MPI_Get_count(&status,tileLayoutInFile,&count);

    //fprintf(stderr,"MPI: %d %d\n",res,count);

    MPI_File_close(&fh);
}


// N.B.: tileID is 1-based, not 0-based
inline int
isInSection1(int tileID)
{ return (tileID <= (7 * tilesPerFacet)); }


// N.B.: tileID is 1-based, not 0-based
long int
tileOffsetInField(int tileID)
{
    return  isInSection1(tileID) ?
((tileID -= 1),
 (((long)tileID / facetTilesInX) * tileZlevelSizeInBytes * facetTilesInX) +
 (((long)tileID % facetTilesInX) * tileSizeX * sizeofFieldElementalType))
                                 :
((tileID -= 1 + (7 * tilesPerFacet)),
 (7 * facetBytes2D) +
 ((tileID / (3*facetTilesInX)) * tileZlevelSizeInBytes * 3*facetTilesInX) +
 ((tileID % (3*facetTilesInX)) * tileSizeX * sizeofFieldElementalType));
}


void
readField(
  MPI_Comm  appComm,
  char *filename,
  MPI_Offset fieldOffsetInFile,
  void *tileBuf,
  int tileID,
  int zLevels)
{
    int writeFlag = 0;
    MPI_Comm  sectionComm = MPI_COMM_NULL;
    MPI_Datatype  section1_ioShape, section2_ioShape;
    MPI_Datatype  inMemoryShape, tileShape, ghostedTileShape;

    int inSection1 = isInSection1(tileID);
    MPI_Offset tileOffsetInFile = fieldOffsetInFile + tileOffsetInField(tileID);


    // Create a type with the "shape" of a tile in memory,
    // with the given number of z-levels.
    MPI_Type_hvector(zLevels, 1, tileZlevelSizeInBytes,
                     tileShape2D, &tileShape);
    MPI_Type_commit(&tileShape);

    // Create a type with the "shape" of a tile in memory,
    // with ghost-cells, with the given number of z-levels.
    MPI_Type_hvector(zLevels, 1, ghostedTileZlevelSizeInBytes,
                     ghostedTileShape2D, &ghostedTileShape);
    MPI_Type_commit(&ghostedTileShape);


    // choose the i/o type
    inMemoryShape = tileShape;


    // Split between section1 tiles and section2 tiles.
    // If a rank has been assigned multiple tiles, it is possible that
    // some of those tiles are in section1, and some are in section2.
    // So we have to dynamically do the comm_split each time because we
    // cannot absolutely guarentee that each rank will always be on the
    // same side of the split every time.
    MPI_Comm_split(appComm, inSection1, 0, &sectionComm);

    memset(tileBuf,-1,tileZlevelSizeInBytes*zLevels);

    if (inSection1) {

        // Create a type with the "shape" of a section1 tile
        // in the file, with the given number of z-levels.
        MPI_Type_hvector(zLevels, 1, fieldZlevelSizeInBytes,
                         section1_ioShape2D, &section1_ioShape);
        MPI_Type_commit(&section1_ioShape);

        //printf("section 1: %d -> %ld (%ld + %ld)\n",tileID,tileOffsetInFile,fieldOffsetInFile,tileOffsetInField(tileID));

        // Do the i/o
        tileIO(sectionComm, filename, tileOffsetInFile, section1_ioShape,
               tileBuf, inMemoryShape, writeFlag);
        // I believe (?) this is needed to ensure consistency when writting
        if (writeFlag)  MPI_Barrier(appComm);

        MPI_Type_free(&section1_ioShape);

    } else {

        // Create a type with the "shape" of a section2 tile
        // in the file, with the given number of z-levels.
        MPI_Type_hvector(zLevels, 1, fieldZlevelSizeInBytes,
                         section2_ioShape2D, &section2_ioShape);
        MPI_Type_commit(&section2_ioShape);

        //printf("section 2: %d -> %ld (%ld + %ld)\n",tileID,tileOffsetInFile,fieldOffsetInFile,tileOffsetInField(tileID));

        // Do the i/o
        // I believe (?) this is needed to ensure consistency when writting
        if (writeFlag)  MPI_Barrier(appComm);
        tileIO(sectionComm, filename, tileOffsetInFile, section2_ioShape,
               tileBuf, inMemoryShape, writeFlag);

        MPI_Type_free(&section2_ioShape);
    }

    /*
    if (tileID==315){
      int i;
      printf("field offset: %ld   ",fieldOffsetInFile);
      printf("tile offset: %ld    ",tileOffsetInField(tileID));
      printf("zlevels: %d     ",zLevels);
      for (i=0;i<10;++i)
        printf("%f ",((double*)tileBuf)[i]);
      printf("\n");
    }
    */

    // Clean up
    MPI_Type_free(&tileShape);
    MPI_Type_free(&ghostedTileShape);
    MPI_Comm_free(&sectionComm);
}


// Fortran interface
// This uses the "usual" method for passing Fortran strings:
// the string length is passed, by value, as an extra "hidden" argument
// after the end of the normal argument list.  So for example, this
// routine would be invoked on the Fortran side like this:
//     call readField(comm, filename, offset, tilebuf, tileid, zlevels)
// This method of passing FOrtran strings is NOT defined by the Fortran
// standard, but it is the method of choice for many compilers, including
// gcc (GNU/Linux), and icc (Intel).
//
// PLEASE NOTE that the "offset" field is of type "MPI_Offset", which
// is synonymous with the Fortran type "integer(kind=MPI_OFFSET_KIND)".
// This will typically be integer*8.  But in particular it is almost
// certainly NOT of type "default integer", which means in particular
// that you CANNOT simply pass a constant (e.g. "0") as the argument,
// since that type will be of the wrong size.
void
readfield_(
  MPI_Fint  *fortranAppComm,
  char *fortranFilename,
  int  *fieldOffsetInFileInPencils,
  void *tileBuf,
  int *tileID,
  int *zLevels,
  int filenameLength)
{
    int i;
    char namebuf[filenameLength+1];
    char *filename = namebuf;

//fprintf (stderr, "In NEW readfield_\n");

    MPI_Offset fieldOffsetInFile = *fieldOffsetInFileInPencils * tileSizeX * sizeofFieldElementalType;

    // Translate the MPI communicator from a Fortran-style handle
    // into a C-style handle.
    MPI_Comm appComm = MPI_Comm_f2c(*fortranAppComm);

    // Translate the Fortran-style string into a C-style string
    //memset(filename, ' ', filenameLength));
    strncpy(filename, fortranFilename, filenameLength);
    for (i = filenameLength;  (i > 0) && (' ' == filename[i-1]);  --i) ;
    filename[i] = '\0';
    //while(' ' == *filename) ++filename;
    assert(strlen(filename) > 0);

    //fprintf(stderr,"%d ::%s:: %d %ld  \n",appComm,filename,filenameLength,fieldOffsetInFile);

    // Make the translated call
    readField(appComm, filename, fieldOffsetInFile, tileBuf, *tileID, *zLevels);
}


// For testing
void initsizesandtypes_(void) {initSizesAndTypes();}


/////////////////////////////////////////////////////////////
// Test case
#if 0

#define FILENAME  "./dataFile"
long int fieldOffsetInFile = 0;

void
doIO(MPI_Comm appComm)
{
    int sizeZ = 3;
    int tile1[sizeZ][tileSizeY][tileSizeX];
    int tile2[sizeZ][tileSizeY][tileSizeX];
    int ghostedTile[sizeZ][tileSizeY + 2*yGhosts][tileSizeX + 2*xGhosts];
    int tileID;
    int i,j,k;

    int  appCommSize, appCommRank;
    MPI_Comm_size(appComm, &appCommSize);
    MPI_Comm_rank(appComm, &appCommRank);

    assert((facetTilesInX * tileSizeX) == facetElements1D);
    assert((facetTilesInY * tileSizeY) == facetElements1D);

    // Ignore the dry tiles ("holes") for the moment
    if (facetTilesInX * facetTilesInY * 13 != appCommSize) {
        if (0 == appCommRank) {
            printf("Unexpected number of ranks: is %d, expected %ld\n",
                   appCommSize, facetTilesInX * facetTilesInY * 13);
        }
    }
    tileID = appCommRank + 1;

#if 0
    // Fill tile1 with distinguished values
    for (k = 0;  k < sizeZ;  ++k) {
        for (j = 0;  j < (tileSizeY + 2*yGhosts);  ++j) {
            for (i = 0;  i < (tileSizeX + 2*xGhosts);  ++i) {
                ghostedTile[k][j][i] = -appCommRank;
            }
        }
    }
    for (k = 0;  k < sizeZ;  ++k) {
        for (j = 0;  j < tileSizeY;  ++j) {
            for (i = 0;  i < tileSizeX;  ++i) {
                tile1[k][j][i] = appCommRank;
                ghostedTile[k][j+yGhosts][i+xGhosts] = appCommRank;
            }
        }
    }
#endif


if (0 == appCommRank) system("/bin/echo -n 'begin io: ' ; date ");
    readField(appComm, FILENAME, 0, ghostedTile, tileID, sizeZ);
if (0 == appCommRank) system("/bin/echo -n 'half: ' ; date ");
    readField(appComm, FILENAME, sizeZ*fieldZlevelSizeInBytes,
              ghostedTile, tileID, sizeZ);

#if 1
    for (k = 0;  k < sizeZ;  ++k) {
        for (j = 0;  j < tileSizeY;  ++j) {
            for (i = 0;  i < tileSizeX;  ++i) {
                int value = ghostedTile[k][j+yGhosts][i+xGhosts];
                if (value != appCommRank) {
                    printf("Fail: %d %d %d:  %d %d\n", k,j,i, value, appCommRank);
                    exit(1);
                }
            }
        }
    }
    if (0 == appCommRank) printf("Verification complete\n");
#endif

MPI_Barrier(appComm);
if (0 == appCommRank) system("/bin/echo -n 'finish: ' ; date ");

    MPI_Finalize();
}


int
main(int argc, char *argv[])
{
    MPI_Comm  appComm = MPI_COMM_NULL;

    MPI_Init(&argc, &argv);
    MPI_Comm_dup(MPI_COMM_WORLD, &appComm);

    initSizesAndTypes();
    doIO(appComm);

    MPI_Finalize();
    return 0;
}
#endif

1	//
2	// MPI IO for MITgcm
3	//
4
5	#include <stdio.h>
6	#include <string.h>
7	#include <assert.h>
8	#include <mpi.h>
9
10	#include "SIZE.h"
11
12
13	// lat-lon-cap decomposition has 13 square facets
14	// facetElements1D is typically 1080, or 2160, or 4320
15
16
17
18	//////////////////////////////////////////////////////////////
19	// These values filled in during "initSizesAndTypes()"
20	MPI_Datatype fieldElementalTypeMPI;
21	size_t sizeofFieldElementalType;
22
23	long int tileSizeX;
24	long int tileSizeY;
25	long int xGhosts;
26	long int yGhosts;
27
28	long int facetElements1D;
29	long int facetBytes2D;
30	long int facetTilesInX;
31	long int facetTilesInY;
32	long int tilesPerFacet;
33
34	long int fieldZlevelSizeInBytes;
35	long int tileZlevelSizeInBytes;
36	long int ghostedTileZlevelSizeInBytes;
37
38	// The first 7 facets all use the same style of layout; here we call
39	// them "section1". The last 6 facets share a layout style, but this
40	// style is different than section1. Here, we call them "section2".
41	MPI_Datatype section1_ioShape2D, section2_ioShape2D;
42	MPI_Datatype tileShape2D, ghostedTileShape2D;
43	MPI_Info ioHints;
44	//////////////////////////////////////////////////////////////
45
46
47
48	int
49	getSizeOfMPIType(MPI_Datatype mpi_type)
50	{
51	switch (mpi_type) {
52
53	case MPI_INT: case MPI_FLOAT: case MPI_REAL4:
54	return 4;
55	break;
56
57	case MPI_LONG_INT: case MPI_DOUBLE: case MPI_REAL8:
58	return 8;
59	break;
60
61	default:
62	assert(("unexpected mpi elemental type", 0));
63	break;
64
65	}
66	return -1;
67	}
68
69
70
71	void
72	createMPItypes(void)
73	{
74	// Create a type with the "shape" of a section1, 2D tile
75	MPI_Type_vector(tileSizeY, tileSizeX, facetElements1D,
76	fieldElementalTypeMPI, &section1_ioShape2D);
77	MPI_Type_commit(&section1_ioShape2D);
78
79	// Create a type with the "shape" of a section2, 2D tile
80	MPI_Type_vector(tileSizeY, tileSizeX, 3*facetElements1D,
81	fieldElementalTypeMPI, &section2_ioShape2D);
82	MPI_Type_commit(&section2_ioShape2D);
83
84
85	// Create a type that describes a 2D tile in memory
86	MPI_Type_vector(tileSizeY, tileSizeX, tileSizeX,
87	fieldElementalTypeMPI, &tileShape2D);
88	MPI_Type_commit(&tileShape2D);
89
90	// Create a type that describes a 2D tile in memory with ghost-cells.
91	int fullDims[2] = {tileSizeX + 2xGhosts, tileSizeY + 2yGhosts};
92	int tileDims[2] = {tileSizeX, tileSizeY};
93	int startElements[2] = {xGhosts, yGhosts};
94	MPI_Type_create_subarray(2, fullDims, tileDims, startElements,
95	MPI_ORDER_FORTRAN, fieldElementalTypeMPI, &ghostedTileShape2D);
96	MPI_Type_commit(&ghostedTileShape2D);
97
98
99	// Set up some possible hints
100	MPI_Info_create(&ioHints);
101	MPI_Info_set(ioHints, "collective_buffering", "true");
102	char blockSize[64];
103	sprintf(blockSize, "%ld", (((long)facetElements1D * 3) *
104	tileSizeY) * sizeofFieldElementalType);
105	MPI_Info_set(ioHints, "cb_block_size", blockSize);
106	}
107
108
109
110	void
111	initSizesAndTypes(void)
112	{
113	/////////////////////////////////////////////
114	// Fundamental values
115	fieldElementalTypeMPI = MPI_DOUBLE;
116	facetElements1D = sFacet;
117	tileSizeX = sNx;
118	tileSizeY = sNy;
119	xGhosts = OLx;
120	yGhosts = OLy;
121	/////////////////////////////////////////////
122
123	// Derived values
124	sizeofFieldElementalType = getSizeOfMPIType(fieldElementalTypeMPI);
125	long int facetElements2D = ((facetElements1D) * (facetElements1D));
126	facetBytes2D = (facetElements2D * sizeofFieldElementalType);
127	fieldZlevelSizeInBytes = (13*(facetBytes2D));
128
129	facetTilesInX = ((facetElements1D)/(tileSizeX));
130	facetTilesInY = ((facetElements1D)/(tileSizeY));
131	tilesPerFacet = ((facetTilesInX)*(facetTilesInY));
132
133	tileZlevelSizeInBytes = tileSizeX * tileSizeY * sizeofFieldElementalType;
134	ghostedTileZlevelSizeInBytes = (tileSizeX + 2xGhosts)
135	(tileSizeY + 2yGhosts)
136	sizeofFieldElementalType;
137
138	// Create the specialized type definitions
139	createMPItypes();
140	}
141
142
143
144	void
145	tileIO(
146	MPI_Comm comm,
147	char *filename,
148	MPI_Offset tileOffsetInFile,
149	MPI_Datatype tileLayoutInFile,
150	void *tileBuf,
151	MPI_Datatype tileLayoutInMemory,
152	int writeFlag)
153	{
154	int fileFlags;
155	MPI_File fh;
156	int (*MPI_IO)();
157
158	int res,count;
159	MPI_Status status;
160
161	if (writeFlag) {
162	fileFlags = MPI_MODE_WRONLY \| MPI_MODE_CREATE;
163	MPI_IO = MPI_File_write_all;
164	} else {
165	fileFlags = MPI_MODE_RDONLY;
166	MPI_IO = MPI_File_read_all;
167	}
168
169	//printf("filename is %s\n",filename);
170
171	MPI_File_open(comm, filename, fileFlags, ioHints, &fh);
172	MPI_File_set_view(fh, tileOffsetInFile, fieldElementalTypeMPI,
173	tileLayoutInFile, "native", ioHints);
174
175
176	// MPI_IO(fh, tileBuf, 1, tileLayoutInMemory, MPI_STATUS_IGNORE);
177	res = MPI_IO(fh, tileBuf, 1, tileLayoutInMemory, &status);
178
179	MPI_Get_count(&status,tileLayoutInFile,&count);
180
181	//fprintf(stderr,"MPI: %d %d\n",res,count);
182
183	MPI_File_close(&fh);
184	}
185
186
187
188	// N.B.: tileID is 1-based, not 0-based
189	inline int
190	isInSection1(int tileID)
191	{ return (tileID <= (7 * tilesPerFacet)); }
192
193
194	// N.B.: tileID is 1-based, not 0-based
195	long int
196	tileOffsetInField(int tileID)
197	{
198	return isInSection1(tileID) ?
199	((tileID -= 1),
200	(((long)tileID / facetTilesInX) * tileZlevelSizeInBytes * facetTilesInX) +
201	(((long)tileID % facetTilesInX) * tileSizeX * sizeofFieldElementalType))
202	:
203	((tileID -= 1 + (7 * tilesPerFacet)),
204	(7 * facetBytes2D) +
205	((tileID / (3facetTilesInX)) tileZlevelSizeInBytes * 3*facetTilesInX) +
206	((tileID % (3facetTilesInX)) tileSizeX * sizeofFieldElementalType));
207	}
208
209
210
211	void
212	readField(
213	MPI_Comm appComm,
214	char *filename,
215	MPI_Offset fieldOffsetInFile,
216	void *tileBuf,
217	int tileID,
218	int zLevels)
219	{
220	int writeFlag = 0;
221	MPI_Comm sectionComm = MPI_COMM_NULL;
222	MPI_Datatype section1_ioShape, section2_ioShape;
223	MPI_Datatype inMemoryShape, tileShape, ghostedTileShape;
224
225	int inSection1 = isInSection1(tileID);
226	MPI_Offset tileOffsetInFile = fieldOffsetInFile + tileOffsetInField(tileID);
227
228
229	// Create a type with the "shape" of a tile in memory,
230	// with the given number of z-levels.
231	MPI_Type_hvector(zLevels, 1, tileZlevelSizeInBytes,
232	tileShape2D, &tileShape);
233	MPI_Type_commit(&tileShape);
234
235	// Create a type with the "shape" of a tile in memory,
236	// with ghost-cells, with the given number of z-levels.
237	MPI_Type_hvector(zLevels, 1, ghostedTileZlevelSizeInBytes,
238	ghostedTileShape2D, &ghostedTileShape);
239	MPI_Type_commit(&ghostedTileShape);
240
241
242	// choose the i/o type
243	inMemoryShape = tileShape;
244
245
246	// Split between section1 tiles and section2 tiles.
247	// If a rank has been assigned multiple tiles, it is possible that
248	// some of those tiles are in section1, and some are in section2.
249	// So we have to dynamically do the comm_split each time because we
250	// cannot absolutely guarentee that each rank will always be on the
251	// same side of the split every time.
252	MPI_Comm_split(appComm, inSection1, 0, &sectionComm);
253
254	memset(tileBuf,-1,tileZlevelSizeInBytes*zLevels);
255
256	if (inSection1) {
257
258	// Create a type with the "shape" of a section1 tile
259	// in the file, with the given number of z-levels.
260	MPI_Type_hvector(zLevels, 1, fieldZlevelSizeInBytes,
261	section1_ioShape2D, &section1_ioShape);
262	MPI_Type_commit(&section1_ioShape);
263
264	//printf("section 1: %d -> %ld (%ld + %ld)\n",tileID,tileOffsetInFile,fieldOffsetInFile,tileOffsetInField(tileID));
265
266	// Do the i/o
267	tileIO(sectionComm, filename, tileOffsetInFile, section1_ioShape,
268	tileBuf, inMemoryShape, writeFlag);
269	// I believe (?) this is needed to ensure consistency when writting
270	if (writeFlag) MPI_Barrier(appComm);
271
272	MPI_Type_free(&section1_ioShape);
273
274	} else {
275
276	// Create a type with the "shape" of a section2 tile
277	// in the file, with the given number of z-levels.
278	MPI_Type_hvector(zLevels, 1, fieldZlevelSizeInBytes,
279	section2_ioShape2D, &section2_ioShape);
280	MPI_Type_commit(&section2_ioShape);
281
282	//printf("section 2: %d -> %ld (%ld + %ld)\n",tileID,tileOffsetInFile,fieldOffsetInFile,tileOffsetInField(tileID));
283
284	// Do the i/o
285	// I believe (?) this is needed to ensure consistency when writting
286	if (writeFlag) MPI_Barrier(appComm);
287	tileIO(sectionComm, filename, tileOffsetInFile, section2_ioShape,
288	tileBuf, inMemoryShape, writeFlag);
289
290	MPI_Type_free(&section2_ioShape);
291	}
292
293	/*
294	if (tileID==315){
295	int i;
296	printf("field offset: %ld ",fieldOffsetInFile);
297	printf("tile offset: %ld ",tileOffsetInField(tileID));
298	printf("zlevels: %d ",zLevels);
299	for (i=0;i<10;++i)
300	printf("%f ",((double*)tileBuf)[i]);
301	printf("\n");
302	}
303	*/
304
305	// Clean up
306	MPI_Type_free(&tileShape);
307	MPI_Type_free(&ghostedTileShape);
308	MPI_Comm_free(&sectionComm);
309	}
310
311
312
313	// Fortran interface
314	// This uses the "usual" method for passing Fortran strings:
315	// the string length is passed, by value, as an extra "hidden" argument
316	// after the end of the normal argument list. So for example, this
317	// routine would be invoked on the Fortran side like this:
318	// call readField(comm, filename, offset, tilebuf, tileid, zlevels)
319	// This method of passing FOrtran strings is NOT defined by the Fortran
320	// standard, but it is the method of choice for many compilers, including
321	// gcc (GNU/Linux), and icc (Intel).
322	//
323	// PLEASE NOTE that the "offset" field is of type "MPI_Offset", which
324	// is synonymous with the Fortran type "integer(kind=MPI_OFFSET_KIND)".
325	// This will typically be integer*8. But in particular it is almost
326	// certainly NOT of type "default integer", which means in particular
327	// that you CANNOT simply pass a constant (e.g. "0") as the argument,
328	// since that type will be of the wrong size.
329	void
330	readfield_(
331	MPI_Fint *fortranAppComm,
332	char *fortranFilename,
333	int *fieldOffsetInFileInPencils,
334	void *tileBuf,
335	int *tileID,
336	int *zLevels,
337	int filenameLength)
338	{
339	int i;
340	char namebuf[filenameLength+1];
341	char *filename = namebuf;
342
343	//fprintf (stderr, "In NEW readfield_\n");
344
345	MPI_Offset fieldOffsetInFile = fieldOffsetInFileInPencils tileSizeX * sizeofFieldElementalType;
346
347	// Translate the MPI communicator from a Fortran-style handle
348	// into a C-style handle.
349	MPI_Comm appComm = MPI_Comm_f2c(*fortranAppComm);
350
351	// Translate the Fortran-style string into a C-style string
352	//memset(filename, ' ', filenameLength));
353	strncpy(filename, fortranFilename, filenameLength);
354	for (i = filenameLength; (i > 0) && (' ' == filename[i-1]); --i) ;
355	filename[i] = '\0';
356	//while(' ' == *filename) ++filename;
357	assert(strlen(filename) > 0);
358
359	//fprintf(stderr,"%d ::%s:: %d %ld \n",appComm,filename,filenameLength,fieldOffsetInFile);
360
361	// Make the translated call
362	readField(appComm, filename, fieldOffsetInFile, tileBuf, tileID, zLevels);
363	}
364
365
366
367	// For testing
368	void initsizesandtypes_(void) {initSizesAndTypes();}
369
370
371
372
373
374	/////////////////////////////////////////////////////////////
375	// Test case
376	#if 0
377
378	#define FILENAME "./dataFile"
379	long int fieldOffsetInFile = 0;
380
381	void
382	doIO(MPI_Comm appComm)
383	{
384	int sizeZ = 3;
385	int tile1[sizeZ][tileSizeY][tileSizeX];
386	int tile2[sizeZ][tileSizeY][tileSizeX];
387	int ghostedTile[sizeZ][tileSizeY + 2yGhosts][tileSizeX + 2xGhosts];
388	int tileID;
389	int i,j,k;
390
391	int appCommSize, appCommRank;
392	MPI_Comm_size(appComm, &appCommSize);
393	MPI_Comm_rank(appComm, &appCommRank);
394
395	assert((facetTilesInX * tileSizeX) == facetElements1D);
396	assert((facetTilesInY * tileSizeY) == facetElements1D);
397
398	// Ignore the dry tiles ("holes") for the moment
399	if (facetTilesInX * facetTilesInY * 13 != appCommSize) {
400	if (0 == appCommRank) {
401	printf("Unexpected number of ranks: is %d, expected %ld\n",
402	appCommSize, facetTilesInX * facetTilesInY * 13);
403	}
404	}
405	tileID = appCommRank + 1;
406
407	#if 0
408	// Fill tile1 with distinguished values
409	for (k = 0; k < sizeZ; ++k) {
410	for (j = 0; j < (tileSizeY + 2*yGhosts); ++j) {
411	for (i = 0; i < (tileSizeX + 2*xGhosts); ++i) {
412	ghostedTile[k][j][i] = -appCommRank;
413	}
414	}
415	}
416	for (k = 0; k < sizeZ; ++k) {
417	for (j = 0; j < tileSizeY; ++j) {
418	for (i = 0; i < tileSizeX; ++i) {
419	tile1[k][j][i] = appCommRank;
420	ghostedTile[k][j+yGhosts][i+xGhosts] = appCommRank;
421	}
422	}
423	}
424	#endif
425
426
427
428	if (0 == appCommRank) system("/bin/echo -n 'begin io: ' ; date ");
429	readField(appComm, FILENAME, 0, ghostedTile, tileID, sizeZ);
430	if (0 == appCommRank) system("/bin/echo -n 'half: ' ; date ");
431	readField(appComm, FILENAME, sizeZ*fieldZlevelSizeInBytes,
432	ghostedTile, tileID, sizeZ);
433
434	#if 1
435	for (k = 0; k < sizeZ; ++k) {
436	for (j = 0; j < tileSizeY; ++j) {
437	for (i = 0; i < tileSizeX; ++i) {
438	int value = ghostedTile[k][j+yGhosts][i+xGhosts];
439	if (value != appCommRank) {
440	printf("Fail: %d %d %d: %d %d\n", k,j,i, value, appCommRank);
441	exit(1);
442	}
443	}
444	}
445	}
446	if (0 == appCommRank) printf("Verification complete\n");
447	#endif
448
449	MPI_Barrier(appComm);
450	if (0 == appCommRank) system("/bin/echo -n 'finish: ' ; date ");
451
452	MPI_Finalize();
453	}
454
455
456
457	int
458	main(int argc, char *argv[])
459	{
460	MPI_Comm appComm = MPI_COMM_NULL;
461
462	MPI_Init(&argc, &argv);
463	MPI_Comm_dup(MPI_COMM_WORLD, &appComm);
464
465	initSizesAndTypes();
466	doIO(appComm);
467
468	MPI_Finalize();
469	return 0;
470	}
471	#endif
472