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	dimitri	1.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	dimitri	1.2	#include "SIZE.h"
11
12	dimitri	1.1
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	dimitri	1.2	facetElements1D = sFacet;
117			tileSizeX = sNx;
118			tileSizeY = sNy;
119			xGhosts = OLx;
120			yGhosts = OLy;
121	dimitri	1.1	/////////////////////////////////////////////
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	dimitri	1.2	//fprintf (stderr, "In NEW readfield_\n");
344
345	dimitri	1.1	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