llc_4320/code-async/readtile_mpiio.c

//
// MPI IO for MITgcm
//

#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <mpi.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);
}


// Somehow we acquire this info at runtime
void
initSizesAndTypes(void)
{
    /////////////////////////////////////////////
    // Fundamental values
    fieldElementalTypeMPI = MPI_DOUBLE;
    facetElements1D = 4320;
    tileSizeX = 72;
    tileSizeY = 72;
    xGhosts = 8;
    yGhosts = 8;
    /////////////////////////////////////////////

    // 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;

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