xref: /phasta/phastaIO/phastaIO.cc (revision 4d60bba28c1e1f3ca80b42756ae9dcbcd5c4bc48)
1 /*
2  *
3  * This is the SyncIO library that uses MPI-IO collective functions to
4  * implement a flexible I/O checkpoint solution for a large number of
5  * processors.
6  *
7  * Previous developer: Ning Liu         (liun2@cs.rpi.edu)
8  *                     Jing Fu          (fuj@cs.rpi.edu)
9  * Current developers: Michel Rasquin   (Michel.Rasquin@colorado.edu),
10  *                     Ben Matthews     (benjamin.a.matthews@colorado.edu)
11  *
12  */
13 
14 #include <map>
15 #include <vector>
16 #include <string>
17 #include <string.h>
18 #include <ctype.h>
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include <math.h>
22 #include <sstream>
23 #include "phastaIO.h"
24 #include "mpi.h"
25 #include "phiotmrc.h"
26 #include <assert.h>
27 
28 #define VERSION_INFO_HEADER_SIZE 8192
29 #define DB_HEADER_SIZE 1024
30 #define ONE_MEGABYTE 1048576
31 #define TWO_MEGABYTE 2097152
32 #define ENDIAN_TEST_NUMBER 12180 // Troy's Zip Code!!
33 #define MAX_PHASTA_FILES 64
34 #define MAX_PHASTA_FILE_NAME_LENGTH 1024
35 #define MAX_FIELDS_NUMBER ((VERSION_INFO_HEADER_SIZE/MAX_FIELDS_NAME_LENGTH)-4) // The meta data include - MPI_IO_Tag, nFields, nFields*names of the fields, nppf
36 // -3 for MPI_IO_Tag, nFields and nppf, -4 for extra security (former nFiles)
37 #define MAX_FIELDS_NAME_LENGTH 128
38 #define DefaultMHSize (4*ONE_MEGABYTE)
39 //#define DefaultMHSize (8350) //For test
40 #define LATEST_WRITE_VERSION 1
41 #define inv1024sq 953.674316406e-9 // = 1/1024/1024
42 int MasterHeaderSize = -1;
43 
44 bool PRINT_PERF = false; // default print no perf results
45 int irank = -1; // global rank, should never be manually manipulated
46 int mysize = -1;
47 
48 // Static variables are bad but used here to store the subcommunicator and associated variables
49 // Prevent MPI_Comm_split to be called more than once, especially on BGQ with the V1R2M1 driver (leak detected in MPI_Comm_split - IBM working on it)
50 static int s_assign_local_comm = 0;
51 static MPI_Comm s_local_comm;
52 static int s_local_size = -1;
53 static int s_local_rank = -1;
54 
55 // the following defines are for debug printf
56 #define PHASTAIO_DEBUG 0 //default to not print any debugging info
57 
58 void phprintf(const char* fmt, ...) {
59   (void)fmt;
60 #if PHASTAIO_DEBUG
61   char format[1024];
62   snprintf(format, sizeof(format), "phastaIO debug: %s", fmt);
63   va_list ap;
64   va_start(ap,fmt);
65   vprintf(format,ap)
66   va_end(ap);
67 #endif
68 }
69 
70 void phprintf_0(const char* fmt, ...) {
71   (void)fmt;
72 #if PHASTAIO_DEBUG
73   int rank = 0;
74   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
75   if(rank == 0){
76     char format[1024];
77     snprintf(format, sizeof(format), "phastaIO debug: irank=0 %s", fmt);
78     va_list ap;
79     va_start(ap,s);
80     vprintf(format, ap);
81     va_end(ap);
82   }
83 #endif
84 }
85 
86 enum PhastaIO_Errors
87 {
88   MAX_PHASTA_FILES_EXCEEDED = -1,
89   UNABLE_TO_OPEN_FILE = -2,
90   NOT_A_MPI_FILE = -3,
91   GPID_EXCEEDED = -4,
92   DATA_TYPE_ILLEGAL = -5
93 };
94 
95 using namespace std;
96 
97 namespace{
98 
99   map<int, std::string> LastHeaderKey;
100   vector< FILE* > fileArray;
101   vector< bool > byte_order;
102   vector< int > header_type;
103   int DataSize=0;
104   bool LastHeaderNotFound = false;
105   bool Wrong_Endian = false ;
106   bool Strict_Error = false ;
107   bool binary_format = true;
108 
109   /***********************************************************************/
110   /***************** NEW PHASTA IO CODE STARTS HERE **********************/
111   /***********************************************************************/
112 
113   typedef struct
114   {
115     char filename[MAX_PHASTA_FILE_NAME_LENGTH];   /* defafults to 1024 */
116     unsigned long my_offset;
117     unsigned long next_start_address;
118     unsigned long **my_offset_table;
119     unsigned long **my_read_table;
120 
121     double * double_chunk;
122     double * read_double_chunk;
123 
124     int field_count;
125     int part_count;
126     int read_field_count;
127     int read_part_count;
128     int GPid;
129     int start_id;
130 
131     int mhsize;
132 
133     int myrank;
134     int numprocs;
135     int local_myrank;
136     int local_numprocs;
137 
138     int nppp;
139     int nPPF;
140     int nFiles;
141     int nFields;
142 
143     int * int_chunk;
144     int * read_int_chunk;
145 
146     int Wrong_Endian; /* default to false */
147     char * master_header;
148     MPI_File file_handle;
149     MPI_Comm local_comm;
150   } phastaio_file_t;
151 
152   typedef struct
153   {
154     int nppf, nfields;
155     char * masterHeader;
156   }serial_file;
157 
158   serial_file *SerialFile;
159   phastaio_file_t *PhastaIOActiveFiles[MAX_PHASTA_FILES];
160   int PhastaIONextActiveIndex = 0; /* indicates next index to allocate */
161 
162   // the caller has the responsibility to delete the returned string
163   // TODO: StringStipper("nbc value? ") returns NULL?
164   char* StringStripper( const char  istring[] ) {
165     int length = strlen( istring );
166     char* dest = (char *)malloc( length + 1 );
167     strcpy( dest, istring );
168     dest[ length ] = '\0';
169     if ( char* p = strpbrk( dest, " ") )
170       *p = '\0';
171     return dest;
172   }
173 
174   inline int cscompare( const char teststring[], const char targetstring[] ) {
175     char* s1 = const_cast<char*>(teststring);
176     char* s2 = const_cast<char*>(targetstring);
177 
178     while( *s1 == ' ') s1++;
179     while( *s2 == ' ') s2++;
180     while( ( *s1 )
181         && ( *s2 )
182         && ( *s2 != '?')
183         && ( tolower( *s1 )==tolower( *s2 ) ) ) {
184       s1++;
185       s2++;
186       while( *s1 == ' ') s1++;
187       while( *s2 == ' ') s2++;
188     }
189     if ( !( *s1 ) || ( *s1 == '?') ) return 1;
190     else return 0;
191   }
192 
193   inline void isBinary( const char iotype[] ) {
194     char* fname = StringStripper( iotype );
195     if ( cscompare( fname, "binary" ) ) binary_format = true;
196     else binary_format = false;
197     free (fname);
198 
199   }
200 
201   inline size_t typeSize( const char typestring[] ) {
202     char* ts1 = StringStripper( typestring );
203     if ( cscompare( "integer", ts1 ) ) {
204       free (ts1);
205       return sizeof(int);
206     } else if ( cscompare( "double", ts1 ) ) {
207       free (ts1);
208       return sizeof( double );
209     } else {
210       free (ts1);
211       fprintf(stderr,"unknown type : %s\n",ts1);
212       return 0;
213     }
214   }
215 
216   int readHeader(
217       FILE*       fileObject,
218         const char  phrase[],
219         int*        params,
220         int         expect ) {
221     char* text_header;
222     char* token;
223     char Line[1024] = "\0";
224     char junk;
225     bool FOUND = false ;
226     int real_length;
227     int skip_size, integer_value;
228     int rewind_count=0;
229 
230     if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) ) {
231       rewind( fileObject );
232       clearerr( fileObject );
233       rewind_count++;
234       fgets( Line, 1024, fileObject );
235     }
236 
237     while( !FOUND  && ( rewind_count < 2 ) )  {
238       if ( ( Line[0] != '\n' ) && ( real_length = strcspn( Line, "#" )) ) {
239         text_header = (char *)malloc( real_length + 1 );
240         strncpy( text_header, Line, real_length );
241         text_header[ real_length ] =static_cast<char>(NULL);
242         token = strtok ( text_header, ":" );
243         assert(token);
244         if( cscompare( phrase , token ) ) {
245           FOUND = true ;
246           token = strtok( NULL, " ,;<>" );
247           assert(token);
248           skip_size = atoi( token );
249           int i;
250           for( i=0; i < expect && ( token = strtok( NULL," ,;<>") ); i++) {
251             assert(token);
252             params[i] = atoi( token );
253           }
254           if ( i < expect ) {
255             fprintf(stderr,
256                 "Aloha Expected # of ints not found for: %s\n",phrase );
257           }
258         } else if ( cscompare(token,"byteorder magic number") ) {
259           if ( binary_format ) {
260             fread((void*)&integer_value,sizeof(int),1,fileObject);
261             fread( &junk, sizeof(char), 1 , fileObject );
262             if ( 362436 != integer_value ) Wrong_Endian = true;
263           } else{
264             fscanf(fileObject, "%d\n", &integer_value );
265           }
266         } else {
267           /* some other header, so just skip over */
268           token = strtok( NULL, " ,;<>" );
269           assert(token);
270           skip_size = atoi( token );
271           if ( binary_format)
272             fseek( fileObject, skip_size, SEEK_CUR );
273           else
274             for( int gama=0; gama < skip_size; gama++ )
275               fgets( Line, 1024, fileObject );
276         }
277         free (text_header);
278       } // end of if before while loop
279 
280       if ( !FOUND )
281         if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) ) {
282           rewind( fileObject );
283           clearerr( fileObject );
284           rewind_count++;
285           fgets( Line, 1024, fileObject );
286         }
287     }
288 
289     if ( !FOUND ) {
290       //fprintf(stderr, "Error: Could not find: %s\n", phrase);
291       if(irank == 0) printf("WARNING: Could not find: %s\n", phrase);
292       return 1;
293     }
294     return 0;
295   }
296 
297 } // end unnamed namespace
298 
299 
300 // begin of publicly visible functions
301 
302 /**
303  * This function takes a long long pointer and assign (start) phiotmrc value to it
304  */
305 void startTimer(double* start) {
306   if( !PRINT_PERF ) return;
307   MPI_Barrier(MPI_COMM_WORLD);
308   *start =  phiotmrc();
309 }
310 
311 /**
312  * This function takes a long long pointer and assign (end) phiotmrc value to it
313  */
314 void endTimer(double* end) {
315   if( !PRINT_PERF ) return;
316   *end = phiotmrc();
317   MPI_Barrier(MPI_COMM_WORLD);
318 }
319 
320 /**
321  * choose to print some performance results (or not) according to
322  * the PRINT_PERF macro
323  */
324 void printPerf(
325     const char* func_name,
326     double start,
327     double end,
328     unsigned long datasize,
329     int printdatainfo,
330     const char* extra_msg) {
331   if( !PRINT_PERF ) return;
332   unsigned long data_size = datasize;
333   double time = end - start;
334   unsigned long isizemin,isizemax,isizetot;
335   double sizemin,sizemax,sizeavg,sizetot,rate;
336   double tmin, tmax, tavg, ttot;
337 
338   MPI_Allreduce(&time, &tmin,1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
339   MPI_Allreduce(&time, &tmax,1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
340   MPI_Allreduce(&time, &ttot,1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
341   tavg = ttot/mysize;
342 
343   if(irank == 0) {
344     if ( PhastaIONextActiveIndex == 0 ) printf("** 1PFPP ");
345     else  printf("** syncIO ");
346     printf("%s(): Tmax = %f sec, Tmin = %f sec, Tavg = %f sec", func_name, tmax, tmin, tavg);
347   }
348 
349   if(printdatainfo == 1) { // if printdatainfo ==1, compute I/O rate and block size
350     MPI_Allreduce(&data_size,&isizemin,1,MPI_LONG_LONG_INT,MPI_MIN,MPI_COMM_WORLD);
351     MPI_Allreduce(&data_size,&isizemax,1,MPI_LONG_LONG_INT,MPI_MAX,MPI_COMM_WORLD);
352     MPI_Allreduce(&data_size,&isizetot,1,MPI_LONG_LONG_INT,MPI_SUM,MPI_COMM_WORLD);
353 
354     sizemin=(double)(isizemin*inv1024sq);
355     sizemax=(double)(isizemax*inv1024sq);
356     sizetot=(double)(isizetot*inv1024sq);
357     sizeavg=(double)(1.0*sizetot/mysize);
358     rate=(double)(1.0*sizetot/tmax);
359 
360     if( irank == 0) {
361       printf(", Rate = %f MB/s [%s] \n \t\t\t"
362              " block size: Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB\n",
363              rate, extra_msg, sizemin,sizemax,sizeavg,sizetot);
364     }
365   }
366   else {
367     if(irank == 0) {
368       printf(" \n");
369       //printf(" (%s) \n", extra_msg);
370     }
371   }
372 }
373 
374 /**
375  * This function is normally called at the beginning of a read operation, before
376  * init function.
377  * This function (uses rank 0) reads out nfields, nppf, master header size,
378  * endianess and allocates for masterHeader string.
379  * These values are essential for following read operations. Rank 0 will bcast
380  * these values to other ranks in the commm world
381  *
382  * If the file set is of old POSIX format, it would throw error and exit
383  */
384 void queryphmpiio(const char filename[],int *nfields, int *nppf)
385 {
386   MPI_Comm_rank(MPI_COMM_WORLD, &irank);
387   MPI_Comm_size(MPI_COMM_WORLD, &mysize);
388 
389   if(irank == 0) {
390     FILE * fileHandle;
391     char* fname = StringStripper( filename );
392 
393     fileHandle = fopen (fname,"rb");
394     if (fileHandle == NULL ) {
395       printf("\nError: File %s doesn't exist! Please check!\n",fname);
396     }
397     else {
398       SerialFile =(serial_file *)calloc( 1,  sizeof( serial_file) );
399       int meta_size_limit = VERSION_INFO_HEADER_SIZE;
400       SerialFile->masterHeader = (char *)malloc( meta_size_limit );
401       fread(SerialFile->masterHeader, 1, meta_size_limit, fileHandle);
402 
403       char read_out_tag[MAX_FIELDS_NAME_LENGTH];
404       char version[MAX_FIELDS_NAME_LENGTH/4];
405       int mhsize;
406       char * token;
407       int magic_number;
408 
409       memcpy( read_out_tag,
410           SerialFile->masterHeader,
411           MAX_FIELDS_NAME_LENGTH-1 );
412 
413       if ( cscompare ("MPI_IO_Tag",read_out_tag) ) {
414         // Test endianess ...
415         memcpy (&magic_number,
416             SerialFile->masterHeader + sizeof("MPI_IO_Tag : ")-1, //-1 sizeof returns the size of the string+1 for "\0"
417             sizeof(int) );                                        // masterheader should look like "MPI_IO_Tag : 12180 " with 12180 in binary format
418 
419         if ( magic_number != ENDIAN_TEST_NUMBER ) {
420           printf("Endian is different!\n");
421           // Will do swap later
422         }
423 
424         // test version, old version, default masterheader size is 4M
425         // newer version masterheader size is read from first line
426         memcpy(version,
427             SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/2,
428             MAX_FIELDS_NAME_LENGTH/4 - 1); //TODO: why -1?
429 
430         if( cscompare ("version",version) ) {
431           // if there is "version" tag in the file, then it is newer format
432           // read master header size from here, otherwise use default
433           // Note: if version is "1", we know mhsize is at 3/4 place...
434 
435           token = strtok(version, ":");
436           token = strtok(NULL, " ,;<>" );
437           int iversion = atoi(token);
438 
439           if( iversion == 1) {
440             memcpy( &mhsize,
441                 SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/4*3 + sizeof("mhsize : ")-1,
442                 sizeof(int));
443             if ( magic_number != ENDIAN_TEST_NUMBER )
444               SwapArrayByteOrder(&mhsize, sizeof(int), 1);
445 
446             if( mhsize > DefaultMHSize ) {
447               //if actual headersize is larger than default, let's re-read
448               free(SerialFile->masterHeader);
449               SerialFile->masterHeader = (char *)malloc(mhsize);
450               fseek(fileHandle, 0, SEEK_SET); // reset the file stream position
451               fread(SerialFile->masterHeader,1,mhsize,fileHandle);
452             }
453           }
454           //TODO: check if this is a valid int??
455           MasterHeaderSize = mhsize;
456         }
457         else { // else it's version 0's format w/o version tag, implicating MHSize=4M
458           MasterHeaderSize = DefaultMHSize;
459         }
460 
461         memcpy( read_out_tag,
462             SerialFile->masterHeader+MAX_FIELDS_NAME_LENGTH+1,
463             MAX_FIELDS_NAME_LENGTH ); //TODO: why +1
464 
465         // Read in # fields ...
466         token = strtok( read_out_tag, ":" );
467         token = strtok( NULL," ,;<>" );
468         *nfields = atoi( token );
469         if ( *nfields > MAX_FIELDS_NUMBER) {
470           printf("Error queryphmpiio: nfields is larger than MAX_FIELDS_NUMBER!\n");
471         }
472         SerialFile->nfields=*nfields; //TODO: sanity check of this int?
473 
474         memcpy( read_out_tag,
475             SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH * 2
476             + *nfields * MAX_FIELDS_NAME_LENGTH,
477             MAX_FIELDS_NAME_LENGTH);
478 
479         token = strtok( read_out_tag, ":" );
480         token = strtok( NULL," ,;<>" );
481         *nppf = atoi( token );
482         SerialFile->nppf=*nppf; //TODO: sanity check of int
483       } // end of if("MPI_IO_TAG")
484       else {
485         printf("Error queryphmpiio: The file you opened is not of syncIO new format, please check! read_out_tag = %s\n",read_out_tag);
486         exit(1);
487       }
488       fclose(fileHandle);
489       free(SerialFile->masterHeader);
490       free(SerialFile);
491     } //end of else
492     free(fname);
493   }
494 
495   // Bcast value to every one
496   MPI_Bcast( nfields, 1, MPI_INT, 0, MPI_COMM_WORLD );
497   MPI_Bcast( nppf, 1, MPI_INT, 0, MPI_COMM_WORLD );
498   MPI_Bcast( &MasterHeaderSize, 1, MPI_INT, 0, MPI_COMM_WORLD );
499   phprintf("Info queryphmpiio: myrank = %d, MasterHeaderSize = %d", irank, MasterHeaderSize);
500 }
501 
502 /**
503  * This function computes the right master header size (round to size of 2^n).
504  * This is only needed for file format version 1 in "write" mode.
505  */
506 int computeMHSize(int nfields, int nppf, int version) {
507   int mhsize=0;
508   if(version == 1) {
509     //int meta_info_size = (2+nfields+1) * MAX_FIELDS_NAME_LENGTH; // 2 is MPI_IO_TAG and nFields, the others 1 is nppf
510     int meta_info_size = VERSION_INFO_HEADER_SIZE;
511     int actual_size =  nfields * nppf * sizeof(unsigned long) + meta_info_size;
512     //printf("actual_size = %d, offset table size = %d\n", actual_size,  nfields * nppf * sizeof(long long));
513     if (actual_size > DefaultMHSize) {
514       mhsize = (int) ceil( (double) actual_size/DefaultMHSize); // it's rounded to ceiling of this value
515       mhsize *= DefaultMHSize;
516     }
517     else {
518       mhsize = DefaultMHSize;
519     }
520   } else {
521     int rank = 0;
522     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
523     if(!rank) {
524       fprintf(stderr,
525           "ERROR invalid version passed to %s... exiting\n", __func__);
526       exit(EXIT_FAILURE);
527     }
528   }
529   return mhsize;
530 }
531 
532 /**
533  * Computes correct color of a rank according to number of files.
534  */
535 extern "C" int computeColor( int myrank, int numprocs, int nfiles) {
536   int color =
537     (int)(myrank / (numprocs / nfiles));
538   return color;
539 }
540 
541 
542 /**
543  * Check the file descriptor.
544  */
545 void checkFileDescriptor(const char fctname[], int*  fileDescriptor ) {
546   if ( *fileDescriptor < 0 ) {
547     printf("Error: File descriptor = %d in %s\n",*fileDescriptor,fctname);
548     exit(1);
549   }
550 }
551 
552 /**
553  * Initialize the file struct members and allocate space for file struct
554  * buffers.
555  *
556  * Note: this function is only called when we are using new format. Old POSIX
557  * format should skip this routine and call openfile() directly instead.
558  */
559 int initphmpiio( int *nfields, int *nppf, int *nfiles, int *filehandle, const char mode[])
560 {
561   // we init irank again in case query not called (e.g. syncIO write case)
562   MPI_Comm_rank(MPI_COMM_WORLD, &irank);
563   MPI_Comm_size(MPI_COMM_WORLD, &mysize);
564 
565   phprintf("Info initphmpiio: entering function, myrank = %d, MasterHeaderSize = %d", irank, MasterHeaderSize);
566 
567   double timer_start, timer_end;
568   startTimer(&timer_start);
569 
570   char* imode = StringStripper( mode );
571 
572   // Note: if it's read, we presume query was called prior to init and
573   // MasterHeaderSize is already set to correct value from parsing header
574   // otherwise it's write then it needs some computation to be set
575   if ( cscompare( "read", imode ) ) {
576     // do nothing
577   }
578   else if( cscompare( "write", imode ) ) {
579     MasterHeaderSize =  computeMHSize(*nfields, *nppf, LATEST_WRITE_VERSION);
580   }
581   else {
582     printf("Error initphmpiio: can't recognize the mode %s", imode);
583     exit(1);
584   }
585   free ( imode );
586 
587   phprintf("Info initphmpiio: myrank = %d, MasterHeaderSize = %d", irank, MasterHeaderSize);
588 
589   int i, j;
590 
591   if( PhastaIONextActiveIndex == MAX_PHASTA_FILES ) {
592     printf("Error initphmpiio: PhastaIONextActiveIndex = MAX_PHASTA_FILES");
593     endTimer(&timer_end);
594     printPerf("initphmpiio", timer_start, timer_end, 0, 0, "");
595     return MAX_PHASTA_FILES_EXCEEDED;
596   }
597   //		else if( PhastaIONextActiveIndex == 0 )  //Hang in debug mode on Intrepid
598   //		{
599   //			for( i = 0; i < MAX_PHASTA_FILES; i++ );
600   //			{
601   //				PhastaIOActiveFiles[i] = NULL;
602   //			}
603   //		}
604 
605 
606   PhastaIOActiveFiles[PhastaIONextActiveIndex] = (phastaio_file_t *)calloc( 1,  sizeof( phastaio_file_t) );
607 
608   i = PhastaIONextActiveIndex;
609   PhastaIONextActiveIndex++;
610 
611   //PhastaIOActiveFiles[i]->next_start_address = 2*TWO_MEGABYTE;
612 
613   PhastaIOActiveFiles[i]->next_start_address = MasterHeaderSize;  // what does this mean??? TODO
614 
615   PhastaIOActiveFiles[i]->Wrong_Endian = false;
616 
617   PhastaIOActiveFiles[i]->nFields = *nfields;
618   PhastaIOActiveFiles[i]->nPPF = *nppf;
619   PhastaIOActiveFiles[i]->nFiles = *nfiles;
620   MPI_Comm_rank(MPI_COMM_WORLD, &(PhastaIOActiveFiles[i]->myrank));
621   MPI_Comm_size(MPI_COMM_WORLD, &(PhastaIOActiveFiles[i]->numprocs));
622 
623 
624   if( *nfiles > 1 ) { // split the ranks according to each mpiio file
625 
626     if ( s_assign_local_comm == 0) { // call mpi_comm_split for the first (and only) time
627 
628       if (PhastaIOActiveFiles[i]->myrank == 0) printf("Building subcommunicator\n");
629 
630       int color = computeColor(PhastaIOActiveFiles[i]->myrank, PhastaIOActiveFiles[i]->numprocs, PhastaIOActiveFiles[i]->nFiles);
631       MPI_Comm_split(MPI_COMM_WORLD,
632           color,
633           PhastaIOActiveFiles[i]->myrank,
634           &(PhastaIOActiveFiles[i]->local_comm));
635       MPI_Comm_size(PhastaIOActiveFiles[i]->local_comm,
636           &(PhastaIOActiveFiles[i]->local_numprocs));
637       MPI_Comm_rank(PhastaIOActiveFiles[i]->local_comm,
638           &(PhastaIOActiveFiles[i]->local_myrank));
639 
640       // back up now these variables so that we do not need to call comm_split again
641       s_local_comm = PhastaIOActiveFiles[i]->local_comm;
642       s_local_size = PhastaIOActiveFiles[i]->local_numprocs;
643       s_local_rank = PhastaIOActiveFiles[i]->local_myrank;
644       s_assign_local_comm = 1;
645     }
646     else { // recycle the subcommunicator
647       if (PhastaIOActiveFiles[i]->myrank == 0) printf("Recycling subcommunicator\n");
648       PhastaIOActiveFiles[i]->local_comm = s_local_comm;
649       PhastaIOActiveFiles[i]->local_numprocs = s_local_size;
650       PhastaIOActiveFiles[i]->local_myrank = s_local_rank;
651     }
652   }
653   else { // *nfiles == 1 here - no need to call mpi_comm_split here
654 
655     if (PhastaIOActiveFiles[i]->myrank == 0) printf("Bypassing subcommunicator\n");
656     PhastaIOActiveFiles[i]->local_comm = MPI_COMM_WORLD;
657     PhastaIOActiveFiles[i]->local_numprocs = PhastaIOActiveFiles[i]->numprocs;
658     PhastaIOActiveFiles[i]->local_myrank = PhastaIOActiveFiles[i]->myrank;
659 
660   }
661 
662   PhastaIOActiveFiles[i]->nppp =
663     PhastaIOActiveFiles[i]->nPPF/PhastaIOActiveFiles[i]->local_numprocs;
664 
665   PhastaIOActiveFiles[i]->start_id = PhastaIOActiveFiles[i]->nPPF *
666     (int)(PhastaIOActiveFiles[i]->myrank/PhastaIOActiveFiles[i]->local_numprocs) +
667     (PhastaIOActiveFiles[i]->local_myrank * PhastaIOActiveFiles[i]->nppp);
668 
669   PhastaIOActiveFiles[i]->my_offset_table =
670     ( unsigned long ** ) calloc( MAX_FIELDS_NUMBER , sizeof( unsigned long *) );
671 
672   PhastaIOActiveFiles[i]->my_read_table =
673     ( unsigned long ** ) calloc( MAX_FIELDS_NUMBER , sizeof( unsigned long *) );
674 
675   for (j=0; j<*nfields; j++)
676   {
677     PhastaIOActiveFiles[i]->my_offset_table[j] =
678       ( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nppp , sizeof( unsigned long) );
679 
680     PhastaIOActiveFiles[i]->my_read_table[j] =
681       ( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nppp , sizeof( unsigned long) );
682   }
683   *filehandle = i;
684 
685   PhastaIOActiveFiles[i]->master_header = (char *)calloc(MasterHeaderSize,sizeof( char ));
686   PhastaIOActiveFiles[i]->double_chunk = (double *)calloc(1,sizeof( double ));
687   PhastaIOActiveFiles[i]->int_chunk = (int *)calloc(1,sizeof( int ));
688   PhastaIOActiveFiles[i]->read_double_chunk = (double *)calloc(1,sizeof( double ));
689   PhastaIOActiveFiles[i]->read_int_chunk = (int *)calloc(1,sizeof( int ));
690 
691   // Time monitoring
692   endTimer(&timer_end);
693   printPerf("initphmpiio", timer_start, timer_end, 0, 0, "");
694 
695   phprintf_0("Info initphmpiio: quiting function");
696 
697   return i;
698 }
699 
700 /**
701  * Destruct the file struct and free buffers allocated in init function.
702  */
703 void finalizephmpiio( int *fileDescriptor )
704 {
705   double timer_start, timer_end;
706   startTimer(&timer_start);
707 
708   int i, j;
709   i = *fileDescriptor;
710   //PhastaIONextActiveIndex--;
711 
712   /* //free the offset table for this phasta file */
713   //for(j=0; j<MAX_FIELDS_NUMBER; j++) //Danger: undefined behavior for my_*_table.[j] not allocated or not initialized to NULL
714   for(j=0; j<PhastaIOActiveFiles[i]->nFields; j++)
715   {
716     free( PhastaIOActiveFiles[i]->my_offset_table[j]);
717     free( PhastaIOActiveFiles[i]->my_read_table[j]);
718   }
719   free ( PhastaIOActiveFiles[i]->my_offset_table );
720   free ( PhastaIOActiveFiles[i]->my_read_table );
721   free ( PhastaIOActiveFiles[i]->master_header );
722   free ( PhastaIOActiveFiles[i]->double_chunk );
723   free ( PhastaIOActiveFiles[i]->int_chunk );
724   free ( PhastaIOActiveFiles[i]->read_double_chunk );
725   free ( PhastaIOActiveFiles[i]->read_int_chunk );
726 
727   if( PhastaIOActiveFiles[i]->nFiles > 1 && s_assign_local_comm ) { // the comm was split
728     if (PhastaIOActiveFiles[i]->myrank == 0) printf("Freeing subcommunicator\n");
729     s_assign_local_comm = 0;
730     MPI_Comm_free(&(PhastaIOActiveFiles[i]->local_comm));
731   }
732 
733   free( PhastaIOActiveFiles[i]);
734 
735   endTimer(&timer_end);
736   printPerf("finalizempiio", timer_start, timer_end, 0, 0, "");
737 
738   PhastaIONextActiveIndex--;
739 }
740 
741 
742 /**
743  * Special init for M2N in order to create a subcommunicator for the reduced solution (requires PRINT_PERF to be false for now)
744  * Initialize the file struct members and allocate space for file struct buffers.
745  *
746  * Note: this function is only called when we are using new format. Old POSIX
747  * format should skip this routine and call openfile() directly instead.
748  */
749 int initphmpiiosub( int *nfields, int *nppf, int *nfiles, int *filehandle, const char mode[],MPI_Comm my_local_comm)
750 {
751   // we init irank again in case query not called (e.g. syncIO write case)
752 
753   MPI_Comm_rank(my_local_comm, &irank);
754   MPI_Comm_size(my_local_comm, &mysize);
755 
756   phprintf("Info initphmpiio: entering function, myrank = %d, MasterHeaderSize = %d", irank, MasterHeaderSize);
757 
758   double timer_start, timer_end;
759   startTimer(&timer_start);
760 
761   char* imode = StringStripper( mode );
762 
763   // Note: if it's read, we presume query was called prior to init and
764   // MasterHeaderSize is already set to correct value from parsing header
765   // otherwise it's write then it needs some computation to be set
766   if ( cscompare( "read", imode ) ) {
767     // do nothing
768   }
769   else if( cscompare( "write", imode ) ) {
770     MasterHeaderSize =  computeMHSize(*nfields, *nppf, LATEST_WRITE_VERSION);
771   }
772   else {
773     printf("Error initphmpiio: can't recognize the mode %s", imode);
774     exit(1);
775   }
776   free ( imode );
777 
778   phprintf("Info initphmpiio: myrank = %d, MasterHeaderSize = %d", irank, MasterHeaderSize);
779 
780   int i, j;
781 
782   if( PhastaIONextActiveIndex == MAX_PHASTA_FILES ) {
783     printf("Error initphmpiio: PhastaIONextActiveIndex = MAX_PHASTA_FILES");
784     endTimer(&timer_end);
785     printPerf("initphmpiio", timer_start, timer_end, 0, 0, "");
786     return MAX_PHASTA_FILES_EXCEEDED;
787   }
788   //		else if( PhastaIONextActiveIndex == 0 )  //Hang in debug mode on Intrepid
789   //		{
790   //			for( i = 0; i < MAX_PHASTA_FILES; i++ );
791   //			{
792   //				PhastaIOActiveFiles[i] = NULL;
793   //			}
794   //		}
795 
796 
797   PhastaIOActiveFiles[PhastaIONextActiveIndex] = (phastaio_file_t *)calloc( 1,  sizeof( phastaio_file_t) );
798 
799   i = PhastaIONextActiveIndex;
800   PhastaIONextActiveIndex++;
801 
802   //PhastaIOActiveFiles[i]->next_start_address = 2*TWO_MEGABYTE;
803 
804   PhastaIOActiveFiles[i]->next_start_address = MasterHeaderSize;  // what does this mean??? TODO
805 
806   PhastaIOActiveFiles[i]->Wrong_Endian = false;
807 
808   PhastaIOActiveFiles[i]->nFields = *nfields;
809   PhastaIOActiveFiles[i]->nPPF = *nppf;
810   PhastaIOActiveFiles[i]->nFiles = *nfiles;
811   MPI_Comm_rank(my_local_comm, &(PhastaIOActiveFiles[i]->myrank));
812   MPI_Comm_size(my_local_comm, &(PhastaIOActiveFiles[i]->numprocs));
813 
814   int color = computeColor(PhastaIOActiveFiles[i]->myrank, PhastaIOActiveFiles[i]->numprocs, PhastaIOActiveFiles[i]->nFiles);
815   MPI_Comm_split(my_local_comm,
816       color,
817       PhastaIOActiveFiles[i]->myrank,
818       &(PhastaIOActiveFiles[i]->local_comm));
819   MPI_Comm_size(PhastaIOActiveFiles[i]->local_comm,
820       &(PhastaIOActiveFiles[i]->local_numprocs));
821   MPI_Comm_rank(PhastaIOActiveFiles[i]->local_comm,
822       &(PhastaIOActiveFiles[i]->local_myrank));
823   PhastaIOActiveFiles[i]->nppp =
824     PhastaIOActiveFiles[i]->nPPF/PhastaIOActiveFiles[i]->local_numprocs;
825 
826   PhastaIOActiveFiles[i]->start_id = PhastaIOActiveFiles[i]->nPPF *
827     (int)(PhastaIOActiveFiles[i]->myrank/PhastaIOActiveFiles[i]->local_numprocs) +
828     (PhastaIOActiveFiles[i]->local_myrank * PhastaIOActiveFiles[i]->nppp);
829 
830   PhastaIOActiveFiles[i]->my_offset_table =
831     ( unsigned long ** ) calloc( MAX_FIELDS_NUMBER , sizeof( unsigned long *) );
832 
833   PhastaIOActiveFiles[i]->my_read_table =
834     ( unsigned long ** ) calloc( MAX_FIELDS_NUMBER , sizeof( unsigned long *) );
835 
836   for (j=0; j<*nfields; j++)
837   {
838     PhastaIOActiveFiles[i]->my_offset_table[j] =
839       ( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nppp , sizeof( unsigned long) );
840 
841     PhastaIOActiveFiles[i]->my_read_table[j] =
842       ( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nppp , sizeof( unsigned long) );
843   }
844   *filehandle = i;
845 
846   PhastaIOActiveFiles[i]->master_header = (char *)calloc(MasterHeaderSize,sizeof( char ));
847   PhastaIOActiveFiles[i]->double_chunk = (double *)calloc(1,sizeof( double ));
848   PhastaIOActiveFiles[i]->int_chunk = (int *)calloc(1,sizeof( int ));
849   PhastaIOActiveFiles[i]->read_double_chunk = (double *)calloc(1,sizeof( double ));
850   PhastaIOActiveFiles[i]->read_int_chunk = (int *)calloc(1,sizeof( int ));
851 
852   // Time monitoring
853   endTimer(&timer_end);
854   printPerf("initphmpiiosub", timer_start, timer_end, 0, 0, "");
855 
856   phprintf_0("Info initphmpiiosub: quiting function");
857 
858   return i;
859 }
860 
861 
862 
863 /** open file for both POSIX and MPI-IO syncIO format.
864  *
865  * If it's old POSIX format, simply call posix fopen().
866  *
867  * If it's MPI-IO foramt:
868  * in "read" mode, it builds the header table that points to the offset of
869  * fields for parts;
870  * in "write" mode, it opens the file with MPI-IO open routine.
871  */
872 void openfile(const char filename[], const char mode[], int*  fileDescriptor )
873 {
874   phprintf_0("Info: entering openfile");
875 
876   double timer_start, timer_end;
877   startTimer(&timer_start);
878 
879   if ( PhastaIONextActiveIndex == 0 )
880   {
881     FILE* file=NULL ;
882     *fileDescriptor = 0;
883     char* fname = StringStripper( filename );
884     char* imode = StringStripper( mode );
885 
886     if ( cscompare( "read", imode ) ) file = fopen(fname, "rb" );
887     else if( cscompare( "write", imode ) ) file = fopen(fname, "wb" );
888     else if( cscompare( "append", imode ) ) file = fopen(fname, "ab" );
889 
890     if ( !file ){
891       fprintf(stderr,"Error openfile: unable to open file %s",fname ) ;
892     } else {
893       fileArray.push_back( file );
894       byte_order.push_back( false );
895       header_type.push_back( sizeof(int) );
896       *fileDescriptor = fileArray.size();
897     }
898     free (fname);
899     free (imode);
900   }
901   else // else it would be parallel I/O, opposed to posix io
902   {
903     char* fname = StringStripper( filename );
904     char* imode = StringStripper( mode );
905     int rc;
906     int i = *fileDescriptor;
907     checkFileDescriptor("openfile",&i);
908     char* token;
909 
910     if ( cscompare( "read", imode ) )
911     {
912       //	      if (PhastaIOActiveFiles[i]->myrank == 0)
913       //                printf("\n **********\nRead open ... ... regular version\n");
914 
915       rc = MPI_File_open( PhastaIOActiveFiles[i]->local_comm,
916           fname,
917           MPI_MODE_RDONLY,
918           MPI_INFO_NULL,
919           &(PhastaIOActiveFiles[i]->file_handle) );
920 
921       if(rc)
922       {
923         *fileDescriptor = UNABLE_TO_OPEN_FILE;
924         int error_string_length;
925         char error_string[4096];
926         MPI_Error_string(rc, error_string, &error_string_length);
927         fprintf(stderr, "Error openfile: Unable to open file %s! MPI reports \"%s\"\n",fname,error_string);
928         endTimer(&timer_end);
929         printPerf("openfile", timer_start, timer_end, 0, 0, "");
930         return;
931       }
932 
933       MPI_Status read_tag_status;
934       char read_out_tag[MAX_FIELDS_NAME_LENGTH];
935       int j;
936       int magic_number;
937 
938       if ( PhastaIOActiveFiles[i]->local_myrank == 0 ) {
939         MPI_File_read_at( PhastaIOActiveFiles[i]->file_handle,
940             0,
941             PhastaIOActiveFiles[i]->master_header,
942             MasterHeaderSize,
943             MPI_CHAR,
944             &read_tag_status );
945       }
946 
947       MPI_Bcast( PhastaIOActiveFiles[i]->master_header,
948           MasterHeaderSize,
949           MPI_CHAR,
950           0,
951           PhastaIOActiveFiles[i]->local_comm );
952 
953       memcpy( read_out_tag,
954           PhastaIOActiveFiles[i]->master_header,
955           MAX_FIELDS_NAME_LENGTH-1 );
956 
957       if ( cscompare ("MPI_IO_Tag",read_out_tag) )
958       {
959         // Test endianess ...
960         memcpy ( &magic_number,
961             PhastaIOActiveFiles[i]->master_header+sizeof("MPI_IO_Tag : ")-1, //-1 sizeof returns the size of the string+1 for "\0"
962             sizeof(int) );                                                   // masterheader should look like "MPI_IO_Tag : 12180 " with 12180 in binary format
963 
964         if ( magic_number != ENDIAN_TEST_NUMBER )
965         {
966           PhastaIOActiveFiles[i]->Wrong_Endian = true;
967         }
968 
969         memcpy( read_out_tag,
970             PhastaIOActiveFiles[i]->master_header+MAX_FIELDS_NAME_LENGTH+1, // TODO: WHY +1???
971             MAX_FIELDS_NAME_LENGTH );
972 
973         // Read in # fields ...
974         token = strtok ( read_out_tag, ":" );
975         token = strtok( NULL," ,;<>" );
976         PhastaIOActiveFiles[i]->nFields = atoi( token );
977 
978         unsigned long **header_table;
979         header_table = ( unsigned long ** )calloc(PhastaIOActiveFiles[i]->nFields, sizeof(unsigned long *));
980 
981         for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ )
982         {
983           header_table[j]=( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nPPF , sizeof( unsigned long));
984         }
985 
986         // Read in the offset table ...
987         for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ )
988         {
989           if ( PhastaIOActiveFiles[i]->local_myrank == 0 ) {
990             memcpy( header_table[j],
991                 PhastaIOActiveFiles[i]->master_header +
992                 VERSION_INFO_HEADER_SIZE +
993                 j * PhastaIOActiveFiles[i]->nPPF * sizeof(unsigned long),
994                 PhastaIOActiveFiles[i]->nPPF * sizeof(unsigned long) );
995           }
996 
997           MPI_Scatter( header_table[j],
998               PhastaIOActiveFiles[i]->nppp,
999               MPI_LONG_LONG_INT,
1000               PhastaIOActiveFiles[i]->my_read_table[j],
1001               PhastaIOActiveFiles[i]->nppp,
1002               MPI_LONG_LONG_INT,
1003               0,
1004               PhastaIOActiveFiles[i]->local_comm );
1005 
1006           // Swap byte order if endianess is different ...
1007           if ( PhastaIOActiveFiles[i]->Wrong_Endian ) {
1008             SwapArrayByteOrder( PhastaIOActiveFiles[i]->my_read_table[j],
1009                 sizeof(unsigned long),
1010                 PhastaIOActiveFiles[i]->nppp );
1011           }
1012         }
1013 
1014         for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ ) {
1015           free ( header_table[j] );
1016         }
1017         free (header_table);
1018 
1019       } // end of if MPI_IO_TAG
1020       else //else not valid MPI file
1021       {
1022         *fileDescriptor = NOT_A_MPI_FILE;
1023         printf("Error openfile: The file %s you opened is not in syncIO new format, please check again! File descriptor = %d, MasterHeaderSize = %d, read_out_tag = %s\n",fname,*fileDescriptor,MasterHeaderSize,read_out_tag);
1024         //Printing MasterHeaderSize is useful to test a compiler bug on Intrepid BGP
1025         endTimer(&timer_end);
1026         printPerf("openfile", timer_start, timer_end, 0, 0, "");
1027         return;
1028       }
1029     } // end of if "read"
1030     else if( cscompare( "write", imode ) )
1031     {
1032       rc = MPI_File_open( PhastaIOActiveFiles[i]->local_comm,
1033           fname,
1034           MPI_MODE_WRONLY | MPI_MODE_CREATE,
1035           MPI_INFO_NULL,
1036           &(PhastaIOActiveFiles[i]->file_handle) );
1037       if(rc != MPI_SUCCESS)
1038       {
1039         *fileDescriptor = UNABLE_TO_OPEN_FILE;
1040         int error_string_length;
1041         char error_string[4096];
1042         MPI_Error_string(rc, error_string, &error_string_length);
1043         fprintf(stderr, "Error openfile: Unable to open file %s! MPI reports \"%s\"\n",fname,error_string);
1044         return;
1045       }
1046     } // end of if "write"
1047     free (fname);
1048     free (imode);
1049   } // end of if FileIndex != 0
1050 
1051   endTimer(&timer_end);
1052   printPerf("openfile", timer_start, timer_end, 0, 0, "");
1053 }
1054 
1055 /** close file for both POSIX and MPI-IO syncIO format.
1056  *
1057  * If it's old POSIX format, simply call posix fclose().
1058  *
1059  * If it's MPI-IO foramt:
1060  * in "read" mode, it simply close file with MPI-IO close routine.
1061  * in "write" mode, rank 0 in each group will re-assemble the master header and
1062  * offset table and write to the beginning of file, then close the file.
1063  */
1064 void closefile( int* fileDescriptor, const char mode[] )
1065 {
1066   double timer_start, timer_end;
1067   startTimer(&timer_start);
1068 
1069   int i = *fileDescriptor;
1070   checkFileDescriptor("closefile",&i);
1071 
1072   if ( PhastaIONextActiveIndex == 0 ) {
1073     char* imode = StringStripper( mode );
1074 
1075     if( cscompare( "write", imode )
1076         || cscompare( "append", imode ) ) {
1077       fflush( fileArray[ *fileDescriptor - 1 ] );
1078     }
1079 
1080     fclose( fileArray[ *fileDescriptor - 1 ] );
1081     free (imode);
1082   }
1083   else {
1084     char* imode = StringStripper( mode );
1085 
1086     //write master header here:
1087     if ( cscompare( "write", imode ) ) {
1088       //	      if ( PhastaIOActiveFiles[i]->nPPF * PhastaIOActiveFiles[i]->nFields < 2*ONE_MEGABYTE/8 )  //SHOULD BE CHECKED
1089       //		MasterHeaderSize = 4*ONE_MEGABYTE;
1090       //	      else
1091       //		MasterHeaderSize = 4*ONE_MEGABYTE + PhastaIOActiveFiles[i]->nPPF * PhastaIOActiveFiles[i]->nFields * 8 - 2*ONE_MEGABYTE;
1092 
1093       MasterHeaderSize = computeMHSize( PhastaIOActiveFiles[i]->nFields, PhastaIOActiveFiles[i]->nPPF, LATEST_WRITE_VERSION);
1094       phprintf_0("Info closefile: myrank = %d, MasterHeaderSize = %d\n", PhastaIOActiveFiles[i]->myrank, MasterHeaderSize);
1095 
1096       MPI_Status write_header_status;
1097       char mpi_tag[MAX_FIELDS_NAME_LENGTH];
1098       char version[MAX_FIELDS_NAME_LENGTH/4];
1099       char mhsize[MAX_FIELDS_NAME_LENGTH/4];
1100       int magic_number = ENDIAN_TEST_NUMBER;
1101 
1102       if ( PhastaIOActiveFiles[i]->local_myrank == 0 )
1103       {
1104         bzero((void*)mpi_tag,MAX_FIELDS_NAME_LENGTH);
1105         sprintf(mpi_tag, "MPI_IO_Tag : ");
1106         memcpy(PhastaIOActiveFiles[i]->master_header,
1107             mpi_tag,
1108             MAX_FIELDS_NAME_LENGTH);
1109 
1110         bzero((void*)version,MAX_FIELDS_NAME_LENGTH/4);
1111         // this version is "1", print version in ASCII
1112         sprintf(version, "version : %d",1);
1113         memcpy(PhastaIOActiveFiles[i]->master_header + MAX_FIELDS_NAME_LENGTH/2,
1114             version,
1115             MAX_FIELDS_NAME_LENGTH/4);
1116 
1117         // master header size is computed using the formula above
1118         bzero((void*)mhsize,MAX_FIELDS_NAME_LENGTH/4);
1119         sprintf(mhsize, "mhsize : ");
1120         memcpy(PhastaIOActiveFiles[i]->master_header + MAX_FIELDS_NAME_LENGTH/4*3,
1121             mhsize,
1122             MAX_FIELDS_NAME_LENGTH/4);
1123 
1124         bzero((void*)mpi_tag,MAX_FIELDS_NAME_LENGTH);
1125         sprintf(mpi_tag,
1126             "\nnFields : %d\n",
1127             PhastaIOActiveFiles[i]->nFields);
1128         memcpy(PhastaIOActiveFiles[i]->master_header+MAX_FIELDS_NAME_LENGTH,
1129             mpi_tag,
1130             MAX_FIELDS_NAME_LENGTH);
1131 
1132         bzero((void*)mpi_tag,MAX_FIELDS_NAME_LENGTH);
1133         sprintf(mpi_tag, "\nnPPF : %d\n", PhastaIOActiveFiles[i]->nPPF);
1134         memcpy( PhastaIOActiveFiles[i]->master_header+
1135             PhastaIOActiveFiles[i]->nFields *
1136             MAX_FIELDS_NAME_LENGTH +
1137             MAX_FIELDS_NAME_LENGTH * 2,
1138             mpi_tag,
1139             MAX_FIELDS_NAME_LENGTH);
1140 
1141         memcpy( PhastaIOActiveFiles[i]->master_header+sizeof("MPI_IO_Tag : ")-1, //-1 sizeof returns the size of the string+1 for "\0"
1142             &magic_number,
1143             sizeof(int));
1144 
1145         memcpy( PhastaIOActiveFiles[i]->master_header+sizeof("mhsize : ") -1 + MAX_FIELDS_NAME_LENGTH/4*3,
1146             &MasterHeaderSize,
1147             sizeof(int));
1148       }
1149 
1150       int j = 0;
1151       unsigned long **header_table;
1152       header_table = ( unsigned long ** )calloc(PhastaIOActiveFiles[i]->nFields, sizeof(unsigned long *));
1153 
1154       for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ ) {
1155         header_table[j]=( unsigned long * ) calloc( PhastaIOActiveFiles[i]->nPPF , sizeof( unsigned long));
1156       }
1157 
1158       //if( irank == 0 ) printf("gonna mpi_gather, myrank = %d\n", irank);
1159       for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ ) {
1160         MPI_Gather( PhastaIOActiveFiles[i]->my_offset_table[j],
1161             PhastaIOActiveFiles[i]->nppp,
1162             MPI_LONG_LONG_INT,
1163             header_table[j],
1164             PhastaIOActiveFiles[i]->nppp,
1165             MPI_LONG_LONG_INT,
1166             0,
1167             PhastaIOActiveFiles[i]->local_comm );
1168       }
1169 
1170       if ( PhastaIOActiveFiles[i]->local_myrank == 0 ) {
1171 
1172         //if( irank == 0 ) printf("gonna memcpy for every procs, myrank = %d\n", irank);
1173         for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ ) {
1174           memcpy ( PhastaIOActiveFiles[i]->master_header +
1175               VERSION_INFO_HEADER_SIZE +
1176               j * PhastaIOActiveFiles[i]->nPPF * sizeof(unsigned long),
1177               header_table[j],
1178               PhastaIOActiveFiles[i]->nPPF * sizeof(unsigned long) );
1179         }
1180 
1181         //if( irank == 0 ) printf("gonna file_write_at(), myrank = %d\n", irank);
1182         MPI_File_write_at( PhastaIOActiveFiles[i]->file_handle,
1183             0,
1184             PhastaIOActiveFiles[i]->master_header,
1185             MasterHeaderSize,
1186             MPI_CHAR,
1187             &write_header_status );
1188       }
1189 
1190       ////free(PhastaIOActiveFiles[i]->master_header);
1191 
1192       for ( j = 0; j < PhastaIOActiveFiles[i]->nFields; j++ ) {
1193         free ( header_table[j] );
1194       }
1195       free (header_table);
1196     }
1197 
1198     //if( irank == 0 ) printf("gonna file_close(), myrank = %d\n", irank);
1199     MPI_File_close( &( PhastaIOActiveFiles[i]->file_handle ) );
1200     free ( imode );
1201   }
1202 
1203   endTimer(&timer_end);
1204   printPerf("closefile_", timer_start, timer_end, 0, 0, "");
1205 }
1206 
1207 int readHeader( FILE* f, const char phrase[],
1208     int* params, int numParams, const char* iotype) {
1209   isBinary(iotype);
1210   return readHeader(f,phrase,params,numParams);
1211 }
1212 
1213 void readheader(
1214     int* fileDescriptor,
1215     const  char keyphrase[],
1216     void* valueArray,
1217     int*  nItems,
1218     const char  datatype[],
1219     const char  iotype[] )
1220 {
1221   double timer_start, timer_end;
1222 
1223   startTimer(&timer_start);
1224 
1225   int i = *fileDescriptor;
1226   checkFileDescriptor("readheader",&i);
1227 
1228   if ( PhastaIONextActiveIndex == 0 ) {
1229     int filePtr = *fileDescriptor - 1;
1230     FILE* fileObject;
1231     int* valueListInt;
1232 
1233     if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
1234       fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
1235       fprintf(stderr,"openfile_ function has to be called before \n") ;
1236       fprintf(stderr,"acessing the file\n ") ;
1237       fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
1238       endTimer(&timer_end);
1239       printPerf("readheader", timer_start, timer_end, 0, 0, "");
1240       return;
1241     }
1242 
1243     LastHeaderKey[filePtr] = std::string(keyphrase);
1244     LastHeaderNotFound = false;
1245 
1246     fileObject = fileArray[ filePtr ] ;
1247     Wrong_Endian = byte_order[ filePtr ];
1248 
1249     isBinary( iotype );
1250     typeSize( datatype );   //redundant call, just avoid a compiler warning.
1251 
1252     // right now we are making the assumption that we will only write integers
1253     // on the header line.
1254 
1255     valueListInt = static_cast< int* >( valueArray );
1256     int ierr = readHeader( fileObject ,
1257         keyphrase,
1258         valueListInt,
1259         *nItems ) ;
1260 
1261     byte_order[ filePtr ] = Wrong_Endian ;
1262 
1263     if ( ierr ) LastHeaderNotFound = true;
1264 
1265     //return ; // don't return, go to the end to print perf
1266   }
1267   else {
1268     int* valueListInt;
1269     valueListInt = static_cast <int*>(valueArray);
1270     char* token = NULL;
1271     bool FOUND = false ;
1272     isBinary( iotype );
1273 
1274     MPI_Status read_offset_status;
1275     char read_out_tag[MAX_FIELDS_NAME_LENGTH];
1276     memset(read_out_tag, '\0', MAX_FIELDS_NAME_LENGTH);
1277     char readouttag[MAX_FIELDS_NUMBER][MAX_FIELDS_NAME_LENGTH];
1278     int j;
1279 
1280     int string_length = strlen( keyphrase );
1281     char* buffer = (char*) malloc ( string_length+1 );
1282 
1283     strcpy ( buffer, keyphrase );
1284     buffer[ string_length ] = '\0';
1285 
1286     char* st2 = strtok ( buffer, "@" );
1287     st2 = strtok (NULL, "@");
1288     PhastaIOActiveFiles[i]->GPid = atoi(st2);
1289     if ( char* p = strpbrk(buffer, "@") )
1290       *p = '\0';
1291 
1292     // Check if the user has input the right GPid
1293     if ( ( PhastaIOActiveFiles[i]->GPid <=
1294           PhastaIOActiveFiles[i]->myrank *
1295           PhastaIOActiveFiles[i]->nppp )||
1296         ( PhastaIOActiveFiles[i]->GPid >
1297           ( PhastaIOActiveFiles[i]->myrank + 1 ) *
1298           PhastaIOActiveFiles[i]->nppp ) )
1299     {
1300       *fileDescriptor = NOT_A_MPI_FILE;
1301       printf("Error readheader: The file is not in syncIO new format, please check! myrank = %d, GPid = %d, nppp = %d, keyphrase = %s\n", PhastaIOActiveFiles[i]->myrank, PhastaIOActiveFiles[i]->GPid, PhastaIOActiveFiles[i]->nppp, keyphrase);
1302       // It is possible atoi could not generate a clear integer from st2 because of additional garbage character in keyphrase
1303       endTimer(&timer_end);
1304       printPerf("readheader", timer_start, timer_end, 0, 0, "");
1305       return;
1306     }
1307 
1308     // Find the field we want ...
1309     //for ( j = 0; j<MAX_FIELDS_NUMBER; j++ )
1310     for ( j = 0; j<PhastaIOActiveFiles[i]->nFields; j++ )
1311     {
1312       memcpy( readouttag[j],
1313           PhastaIOActiveFiles[i]->master_header + j*MAX_FIELDS_NAME_LENGTH+MAX_FIELDS_NAME_LENGTH*2+1,
1314           MAX_FIELDS_NAME_LENGTH-1 );
1315     }
1316 
1317     for ( j = 0; j<PhastaIOActiveFiles[i]->nFields; j++ )
1318     {
1319       token = strtok ( readouttag[j], ":" );
1320 
1321       //if ( cscompare( buffer, token ) )
1322       if ( cscompare( token , buffer ) && cscompare( buffer, token ) )
1323         // This double comparison is required for the field "number of nodes" and all the other fields that start with "number of nodes" (i.g. number of nodes in the mesh").
1324         // Would be safer to rename "number of nodes" by "number of nodes in the part" so that the name are completely unique. But much more work to do that (Nspre, phParAdapt, etc).
1325         // Since the field name are unique in SyncIO (as it includes part ID), this should be safe and there should be no issue with the "?" trailing character.
1326       {
1327         PhastaIOActiveFiles[i]->read_field_count = j;
1328         FOUND = true;
1329         //printf("buffer: %s | token: %s | j: %d\n",buffer,token,j);
1330         break;
1331       }
1332     }
1333     free(buffer);
1334 
1335     if (!FOUND)
1336     {
1337       //if(irank==0) printf("Warning readheader: Not found %s \n",keyphrase); //PhastaIOActiveFiles[i]->myrank is certainly initialized here.
1338       if(PhastaIOActiveFiles[i]->myrank == 0) printf("WARNING readheader: Not found %s\n",keyphrase);
1339       endTimer(&timer_end);
1340       printPerf("readheader", timer_start, timer_end, 0, 0, "");
1341       return;
1342     }
1343 
1344     // Find the part we want ...
1345     PhastaIOActiveFiles[i]->read_part_count = PhastaIOActiveFiles[i]->GPid -
1346       PhastaIOActiveFiles[i]->myrank * PhastaIOActiveFiles[i]->nppp - 1;
1347 
1348     PhastaIOActiveFiles[i]->my_offset =
1349       PhastaIOActiveFiles[i]->my_read_table[PhastaIOActiveFiles[i]->read_field_count][PhastaIOActiveFiles[i]->read_part_count];
1350 
1351     // 	  printf("****Rank %d offset is %d\n",PhastaIOActiveFiles[i]->myrank,PhastaIOActiveFiles[i]->my_offset);
1352 
1353     // Read each datablock header here ...
1354 
1355     MPI_File_read_at_all( PhastaIOActiveFiles[i]->file_handle,
1356         PhastaIOActiveFiles[i]->my_offset+1,
1357         read_out_tag,
1358         MAX_FIELDS_NAME_LENGTH-1,
1359         MPI_CHAR,
1360         &read_offset_status );
1361     token = strtok ( read_out_tag, ":" );
1362 
1363     // 	  printf("&&&&Rank %d read_out_tag is %s\n",PhastaIOActiveFiles[i]->myrank,read_out_tag);
1364 
1365     if( cscompare( keyphrase , token ) ) //No need to compare also token with keyphrase like above. We should already have the right one. Otherwise there is a problem.
1366     {
1367       FOUND = true ;
1368       token = strtok( NULL, " ,;<>" );
1369       for( j=0; j < *nItems && ( token = strtok( NULL," ,;<>") ); j++ )
1370         valueListInt[j] = atoi( token );
1371 
1372       if ( j < *nItems )
1373       {
1374         fprintf( stderr, "Expected # of ints not found for: %s\n", keyphrase );
1375       }
1376     }
1377     else {
1378       //if(irank==0)
1379       if(PhastaIOActiveFiles[i]->myrank == 0)
1380         // If we enter this if, there is a problem with the name of some fields
1381       {
1382         printf("Error readheader: Unexpected mismatch between keyphrase = %s and token = %s\n",keyphrase,token);
1383       }
1384     }
1385   }
1386 
1387   endTimer(&timer_end);
1388   printPerf("readheader", timer_start, timer_end, 0, 0, "");
1389 
1390 }
1391 
1392 void readDataBlock(
1393     FILE* fileObject,
1394     void* valueArray,
1395     int nItems,
1396     const char  datatype[],
1397     const char  iotype[] )
1398 {
1399   isBinary(iotype);
1400   size_t type_size = typeSize( datatype );
1401   if ( binary_format ) {
1402     char junk = '\0';
1403     fread( valueArray, type_size, nItems, fileObject );
1404     fread( &junk, sizeof(char), 1 , fileObject );
1405     if ( Wrong_Endian ) SwapArrayByteOrder( valueArray, type_size, nItems );
1406   } else {
1407     char* ts1 = StringStripper( datatype );
1408     if ( cscompare( "integer", ts1 ) ) {
1409       for( int n=0; n < nItems ; n++ )
1410         fscanf(fileObject, "%d\n",(int*)((int*)valueArray+n) );
1411     } else if ( cscompare( "double", ts1 ) ) {
1412       for( int n=0; n < nItems ; n++ )
1413         fscanf(fileObject, "%lf\n",(double*)((double*)valueArray+n) );
1414     }
1415     free (ts1);
1416   }
1417 }
1418 
1419 void readdatablock(
1420     int*  fileDescriptor,
1421     const char keyphrase[],
1422     void* valueArray,
1423     int*  nItems,
1424     const char  datatype[],
1425     const char  iotype[] )
1426 {
1427   //if(irank == 0) printf("entering readdatablock()\n");
1428   unsigned long data_size = 0;
1429   double timer_start, timer_end;
1430   startTimer(&timer_start);
1431 
1432   int i = *fileDescriptor;
1433   checkFileDescriptor("readdatablock",&i);
1434 
1435   if ( PhastaIONextActiveIndex == 0 ) {
1436     int filePtr = *fileDescriptor - 1;
1437     FILE* fileObject;
1438 
1439     if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
1440       fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
1441       fprintf(stderr,"openfile_ function has to be called before\n") ;
1442       fprintf(stderr,"acessing the file\n ") ;
1443       fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
1444       endTimer(&timer_end);
1445       printPerf("readdatablock", timer_start, timer_end, 0, 0, "");
1446       return;
1447     }
1448 
1449     // error check..
1450     // since we require that a consistant header always preceed the data block
1451     // let us check to see that it is actually the case.
1452 
1453     if ( ! cscompare( LastHeaderKey[ filePtr ].c_str(), keyphrase ) ) {
1454       fprintf(stderr, "Header not consistant with data block\n");
1455       fprintf(stderr, "Header: %s\n", LastHeaderKey[ filePtr ].c_str() );
1456       fprintf(stderr, "DataBlock: %s\n ", keyphrase );
1457       fprintf(stderr, "Please recheck read sequence \n");
1458       if( Strict_Error ) {
1459         fprintf(stderr, "fatal error: cannot continue, returning out of call\n");
1460         endTimer(&timer_end);
1461         printPerf("readdatablock", timer_start, timer_end, 0, 0, "");
1462         return;
1463       }
1464     }
1465 
1466     if ( LastHeaderNotFound ) {
1467       endTimer(&timer_end);
1468       printPerf("readdatablock", timer_start, timer_end, 0, 0, "");
1469       return;
1470     }
1471     fileObject = fileArray[ filePtr ];
1472     Wrong_Endian = byte_order[ filePtr ];
1473     LastHeaderKey.erase(filePtr);
1474     readDataBlock(fileObject,valueArray,*nItems,datatype,iotype);
1475 
1476     //return;
1477   }
1478   else {
1479     // 	  printf("read data block\n");
1480     MPI_Status read_data_status;
1481     size_t type_size = typeSize( datatype );
1482     int nUnits = *nItems;
1483     isBinary( iotype );
1484 
1485     // read datablock then
1486     //MR CHANGE
1487     //             if ( cscompare ( datatype, "double"))
1488     char* ts2 = StringStripper( datatype );
1489     if ( cscompare ( "double" , ts2))
1490       //MR CHANGE END
1491     {
1492 
1493       MPI_File_read_at_all_begin( PhastaIOActiveFiles[i]->file_handle,
1494           PhastaIOActiveFiles[i]->my_offset + DB_HEADER_SIZE,
1495           valueArray,
1496           nUnits,
1497           MPI_DOUBLE );
1498       MPI_File_read_at_all_end( PhastaIOActiveFiles[i]->file_handle,
1499           valueArray,
1500           &read_data_status );
1501       data_size=8*nUnits;
1502 
1503     }
1504     //MR CHANGE
1505     //             else if ( cscompare ( datatype, "integer"))
1506     else if ( cscompare ( "integer" , ts2))
1507       //MR CHANGE END
1508     {
1509       MPI_File_read_at_all_begin(PhastaIOActiveFiles[i]->file_handle,
1510           PhastaIOActiveFiles[i]->my_offset + DB_HEADER_SIZE,
1511           valueArray,
1512           nUnits,
1513           MPI_INT );
1514       MPI_File_read_at_all_end( PhastaIOActiveFiles[i]->file_handle,
1515           valueArray,
1516           &read_data_status );
1517       data_size=4*nUnits;
1518     }
1519     else
1520     {
1521       *fileDescriptor = DATA_TYPE_ILLEGAL;
1522       printf("readdatablock - DATA_TYPE_ILLEGAL - %s\n",datatype);
1523       endTimer(&timer_end);
1524       printPerf("readdatablock", timer_start, timer_end, 0, 0, "");
1525       return;
1526     }
1527     free(ts2);
1528 
1529 
1530     // 	  printf("%d Read finishe\n",PhastaIOActiveFiles[i]->myrank);
1531 
1532     // Swap data byte order if endianess is different ...
1533     if ( PhastaIOActiveFiles[i]->Wrong_Endian )
1534     {
1535       SwapArrayByteOrder( valueArray, type_size, nUnits );
1536     }
1537   }
1538 
1539   endTimer(&timer_end);
1540   char extra_msg[1024];
1541   memset(extra_msg, '\0', 1024);
1542   char* key = StringStripper(keyphrase);
1543   sprintf(extra_msg, " field is %s ", key);
1544   printPerf("readdatablock", timer_start, timer_end, data_size, 1, extra_msg);
1545   free(key);
1546 
1547 }
1548 
1549 void writeHeader(
1550     FILE* f,
1551     const char keyphrase[],
1552     const void* valueArray,
1553     const int nItems,
1554     const int ndataItems,
1555     const char datatype[],
1556     const char iotype[])
1557 {
1558   isBinary( iotype );
1559 
1560   const int _newline =
1561     ( ndataItems > 0 ) ? sizeof( char ) : 0;
1562   int size_of_nextblock =
1563     ( binary_format ) ? typeSize(datatype) * ndataItems + _newline : ndataItems;
1564 
1565   fprintf( f, "%s : < %d > ", keyphrase, size_of_nextblock );
1566   for( int i = 0; i < nItems; i++ )
1567     fprintf( f, "%d ", *((int*)((int*)valueArray+i)));
1568   fprintf( f, "\n");
1569 }
1570 
1571 void writeheader(
1572     const int* fileDescriptor,
1573     const char keyphrase[],
1574     const void* valueArray,
1575     const int* nItems,
1576     const int* ndataItems,
1577     const char datatype[],
1578     const char iotype[])
1579 {
1580 
1581   //if(irank == 0) printf("entering writeheader()\n");
1582 
1583   double timer_start, timer_end;
1584   startTimer(&timer_start);
1585 
1586   int i = *fileDescriptor;
1587   checkFileDescriptor("writeheader",&i);
1588 
1589   if ( PhastaIONextActiveIndex == 0 ) {
1590     int filePtr = *fileDescriptor - 1;
1591     if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
1592       fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
1593       fprintf(stderr,"openfile_ function has to be called before \n") ;
1594       fprintf(stderr,"acessing the file\n ") ;
1595       fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
1596       endTimer(&timer_end);
1597       printPerf("writeheader", timer_start, timer_end, 0, 0, "");
1598       return;
1599     }
1600 
1601     LastHeaderKey[filePtr] = std::string(keyphrase);
1602     DataSize = *ndataItems;
1603     FILE* fileObject = fileArray[ filePtr ] ;
1604     header_type[ filePtr ] = typeSize( datatype );
1605     writeHeader(fileObject,keyphrase,valueArray,*nItems,
1606         *ndataItems,datatype,iotype);
1607   }
1608   else { // else it's parallel I/O
1609     DataSize = *ndataItems;
1610     size_t type_size = typeSize( datatype );
1611     isBinary( iotype );
1612     char mpi_tag[MAX_FIELDS_NAME_LENGTH];
1613 
1614     int string_length = strlen( keyphrase );
1615     char* buffer = (char*) malloc ( string_length+1 );
1616 
1617     strcpy ( buffer, keyphrase);
1618     buffer[ string_length ] = '\0';
1619 
1620     char* st2 = strtok ( buffer, "@" );
1621     st2 = strtok (NULL, "@");
1622     PhastaIOActiveFiles[i]->GPid = atoi(st2);
1623 
1624     if ( char* p = strpbrk(buffer, "@") )
1625       *p = '\0';
1626 
1627     bzero((void*)mpi_tag,MAX_FIELDS_NAME_LENGTH);
1628     sprintf(mpi_tag, "\n%s : %d\n", buffer, PhastaIOActiveFiles[i]->field_count);
1629     unsigned long offset_value;
1630 
1631     int temp = *ndataItems;
1632     unsigned long number_of_items = (unsigned long)temp;
1633     MPI_Barrier(PhastaIOActiveFiles[i]->local_comm);
1634 
1635     MPI_Scan( &number_of_items,
1636         &offset_value,
1637         1,
1638         MPI_LONG_LONG_INT,
1639         MPI_SUM,
1640         PhastaIOActiveFiles[i]->local_comm );
1641 
1642     offset_value = (offset_value - number_of_items) * type_size;
1643 
1644     offset_value += PhastaIOActiveFiles[i]->local_myrank *
1645       DB_HEADER_SIZE +
1646       PhastaIOActiveFiles[i]->next_start_address;
1647     // This offset is the starting address of each datablock header...
1648     PhastaIOActiveFiles[i]->my_offset = offset_value;
1649 
1650     // Write in my offset table ...
1651     PhastaIOActiveFiles[i]->my_offset_table[PhastaIOActiveFiles[i]->field_count][PhastaIOActiveFiles[i]->part_count] =
1652       PhastaIOActiveFiles[i]->my_offset;
1653 
1654     // Update the next-start-address ...
1655     PhastaIOActiveFiles[i]->next_start_address = offset_value +
1656       number_of_items * type_size +
1657       DB_HEADER_SIZE;
1658     MPI_Bcast( &(PhastaIOActiveFiles[i]->next_start_address),
1659         1,
1660         MPI_LONG_LONG_INT,
1661         PhastaIOActiveFiles[i]->local_numprocs-1,
1662         PhastaIOActiveFiles[i]->local_comm );
1663 
1664     // Prepare datablock header ...
1665     int _newline = (*ndataItems>0)?sizeof(char):0;
1666     unsigned int size_of_nextblock = type_size * (*ndataItems) + _newline;
1667 
1668     //char datablock_header[255];
1669     //bzero((void*)datablock_header,255);
1670     char datablock_header[DB_HEADER_SIZE];
1671     bzero((void*)datablock_header,DB_HEADER_SIZE);
1672 
1673     PhastaIOActiveFiles[i]->GPid = PhastaIOActiveFiles[i]->nppp*PhastaIOActiveFiles[i]->myrank+PhastaIOActiveFiles[i]->part_count;
1674     sprintf( datablock_header,
1675         "\n%s : < %u >",
1676         keyphrase,
1677         size_of_nextblock );
1678 
1679     for ( int j = 0; j < *nItems; j++ )
1680     {
1681       sprintf( datablock_header,
1682           "%s %d ",
1683           datablock_header,
1684           *((int*)((int*)valueArray+j)));
1685     }
1686     sprintf( datablock_header,
1687         "%s\n ",
1688         datablock_header );
1689 
1690     // Write datablock header ...
1691     //MR CHANGE
1692     // 	if ( cscompare(datatype,"double") )
1693     char* ts1 = StringStripper( datatype );
1694     if ( cscompare("double",ts1) )
1695       //MR CHANGE END
1696     {
1697       free ( PhastaIOActiveFiles[i]->double_chunk );
1698       PhastaIOActiveFiles[i]->double_chunk = ( double * )malloc( (sizeof( double )*number_of_items+ DB_HEADER_SIZE));
1699 
1700       double * aa = ( double * )datablock_header;
1701       memcpy(PhastaIOActiveFiles[i]->double_chunk, aa, DB_HEADER_SIZE);
1702     }
1703     //MR CHANGE
1704     // 	if  ( cscompare(datatype,"integer") )
1705     else if ( cscompare("integer",ts1) )
1706       //MR CHANGE END
1707     {
1708       free ( PhastaIOActiveFiles[i]->int_chunk );
1709       PhastaIOActiveFiles[i]->int_chunk = ( int * )malloc( (sizeof( int )*number_of_items+ DB_HEADER_SIZE));
1710 
1711       int * aa = ( int * )datablock_header;
1712       memcpy(PhastaIOActiveFiles[i]->int_chunk, aa, DB_HEADER_SIZE);
1713     }
1714     else {
1715       //             *fileDescriptor = DATA_TYPE_ILLEGAL;
1716       printf("writeheader - DATA_TYPE_ILLEGAL - %s\n",datatype);
1717       endTimer(&timer_end);
1718       printPerf("writeheader", timer_start, timer_end, 0, 0, "");
1719       return;
1720     }
1721     free(ts1);
1722 
1723     PhastaIOActiveFiles[i]->part_count++;
1724     if ( PhastaIOActiveFiles[i]->part_count == PhastaIOActiveFiles[i]->nppp ) {
1725       //A new field will be written
1726       if ( PhastaIOActiveFiles[i]->local_myrank == 0 ) {
1727         memcpy( PhastaIOActiveFiles[i]->master_header +
1728             PhastaIOActiveFiles[i]->field_count *
1729             MAX_FIELDS_NAME_LENGTH +
1730             MAX_FIELDS_NAME_LENGTH * 2,
1731             mpi_tag,
1732             MAX_FIELDS_NAME_LENGTH-1);
1733       }
1734       PhastaIOActiveFiles[i]->field_count++;
1735       PhastaIOActiveFiles[i]->part_count=0;
1736     }
1737     free(buffer);
1738   }
1739 
1740   endTimer(&timer_end);
1741   printPerf("writeheader", timer_start, timer_end, 0, 0, "");
1742 }
1743 
1744 void writeDataBlock(
1745     FILE* f,
1746     const void* valueArray,
1747     const int   nItems,
1748     const char  datatype[],
1749     const char  iotype[]  )
1750 {
1751   isBinary( iotype );
1752   size_t type_size = typeSize( datatype );
1753   if ( binary_format ) {
1754     fwrite( valueArray, type_size, nItems, f );
1755     fprintf( f,"\n");
1756   } else {
1757     char* ts1 = StringStripper( datatype );
1758     if ( cscompare( "integer", ts1 ) ) {
1759       const int* vals = (int*) valueArray;
1760       for( int n=0; n < nItems ; n++ )
1761         fprintf(f,"%d\n",vals[n]);
1762     } else if ( cscompare( "double", ts1 ) ) {
1763       const double* vals = (double*) valueArray;
1764       for( int n=0; n < nItems ; n++ )
1765         fprintf(f,"%f\n",vals[n]);
1766     }
1767     free (ts1);
1768   }
1769 }
1770 
1771 void writedatablock(
1772     const int* fileDescriptor,
1773     const char keyphrase[],
1774     const void* valueArray,
1775     const int* nItems,
1776     const char datatype[],
1777     const char iotype[] )
1778 {
1779   //if(irank == 0) printf("entering writedatablock()\n");
1780 
1781   unsigned long data_size = 0;
1782   double timer_start, timer_end;
1783   startTimer(&timer_start);
1784 
1785   int i = *fileDescriptor;
1786   checkFileDescriptor("writedatablock",&i);
1787 
1788   if ( PhastaIONextActiveIndex == 0 ) {
1789     int filePtr = *fileDescriptor - 1;
1790 
1791     if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
1792       fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
1793       fprintf(stderr,"openfile_ function has to be called before \n") ;
1794       fprintf(stderr,"acessing the file\n ") ;
1795       fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
1796       endTimer(&timer_end);
1797       printPerf("writedatablock", timer_start, timer_end, 0, 0, "");
1798       return;
1799     }
1800     // since we require that a consistant header always preceed the data block
1801     // let us check to see that it is actually the case.
1802 
1803     if ( ! cscompare( LastHeaderKey[ filePtr ].c_str(), keyphrase ) ) {
1804       fprintf(stderr, "Header not consistant with data block\n");
1805       fprintf(stderr, "Header: %s\n", LastHeaderKey[ filePtr ].c_str() );
1806       fprintf(stderr, "DataBlock: %s\n ", keyphrase );
1807       fprintf(stderr, "Please recheck write sequence \n");
1808       if( Strict_Error ) {
1809         fprintf(stderr, "fatal error: cannot continue, returning out of call\n");
1810         endTimer(&timer_end);
1811         printPerf("writedatablock", timer_start, timer_end, 0, 0, "");
1812         return;
1813       }
1814     }
1815 
1816     FILE* fileObject =  fileArray[ filePtr ] ;
1817     size_t type_size=typeSize( datatype );
1818     isBinary( iotype );
1819 
1820     LastHeaderKey.erase(filePtr);
1821 
1822     if ( header_type[filePtr] != (int)type_size ) {
1823       fprintf(stderr,"header and datablock differ on typeof data in the block for\n");
1824       fprintf(stderr,"keyphrase : %s\n", keyphrase);
1825       if( Strict_Error ) {
1826         fprintf(stderr,"fatal error: cannot continue, returning out of call\n" );
1827         endTimer(&timer_end);
1828         printPerf("writedatablock", timer_start, timer_end, 0, 0, "");
1829         return;
1830       }
1831     }
1832 
1833     int nUnits = *nItems;
1834 
1835     if ( nUnits != DataSize ) {
1836       fprintf(stderr,"header and datablock differ on number of data items for\n");
1837       fprintf(stderr,"keyphrase : %s\n", keyphrase);
1838       if( Strict_Error ) {
1839         fprintf(stderr,"fatal error: cannot continue, returning out of call\n" );
1840         endTimer(&timer_end);
1841         printPerf("writedatablock", timer_start, timer_end, 0, 0, "");
1842         return;
1843       }
1844     }
1845     writeDataBlock(fileObject,valueArray,*nItems,datatype,iotype);
1846   }
1847   else {  // syncIO case
1848     MPI_Status write_data_status;
1849     isBinary( iotype );
1850     int nUnits = *nItems;
1851 
1852     //MR CHANGE
1853     // 	if ( cscompare(datatype,"double") )
1854     char* ts1 = StringStripper( datatype );
1855     if ( cscompare("double",ts1) )
1856       //MR CHANGE END
1857     {
1858       memcpy((PhastaIOActiveFiles[i]->double_chunk+DB_HEADER_SIZE/sizeof(double)), valueArray, nUnits*sizeof(double));
1859       MPI_File_write_at_all_begin( PhastaIOActiveFiles[i]->file_handle,
1860           PhastaIOActiveFiles[i]->my_offset,
1861           PhastaIOActiveFiles[i]->double_chunk,
1862           //BLOCK_SIZE/sizeof(double),
1863           nUnits+DB_HEADER_SIZE/sizeof(double),
1864           MPI_DOUBLE );
1865       MPI_File_write_at_all_end( PhastaIOActiveFiles[i]->file_handle,
1866           PhastaIOActiveFiles[i]->double_chunk,
1867           &write_data_status );
1868       data_size=8*nUnits;
1869     }
1870     //MR CHANGE
1871     // 	else if ( cscompare ( datatype, "integer"))
1872     else if ( cscompare("integer",ts1) )
1873       //MR CHANGE END
1874     {
1875       memcpy((PhastaIOActiveFiles[i]->int_chunk+DB_HEADER_SIZE/sizeof(int)), valueArray, nUnits*sizeof(int));
1876       MPI_File_write_at_all_begin( PhastaIOActiveFiles[i]->file_handle,
1877           PhastaIOActiveFiles[i]->my_offset,
1878           PhastaIOActiveFiles[i]->int_chunk,
1879           nUnits+DB_HEADER_SIZE/sizeof(int),
1880           MPI_INT );
1881       MPI_File_write_at_all_end( PhastaIOActiveFiles[i]->file_handle,
1882           PhastaIOActiveFiles[i]->int_chunk,
1883           &write_data_status );
1884       data_size=4*nUnits;
1885     }
1886     else {
1887       printf("Error: writedatablock - DATA_TYPE_ILLEGAL - %s\n",datatype);
1888       endTimer(&timer_end);
1889       printPerf("writedatablock", timer_start, timer_end, 0, 0, "");
1890       return;
1891     }
1892     free(ts1);
1893   }
1894 
1895   endTimer(&timer_end);
1896   char extra_msg[1024];
1897   memset(extra_msg, '\0', 1024);
1898   char* key = StringStripper(keyphrase);
1899   sprintf(extra_msg, " field is %s ", key);
1900   printPerf("writedatablock", timer_start, timer_end, data_size, 1, extra_msg);
1901   free(key);
1902 
1903 }
1904 
1905 void SwapArrayByteOrder( void* array, int nbytes, int nItems )
1906 {
1907   /* This swaps the byte order for the array of nItems each
1908      of size nbytes , This will be called only locally  */
1909   int i,j;
1910   unsigned char* ucDst = (unsigned char*)array;
1911   for(i=0; i < nItems; i++) {
1912     for(j=0; j < (nbytes/2); j++)
1913       std::swap( ucDst[j] , ucDst[(nbytes - 1) - j] );
1914     ucDst += nbytes;
1915   }
1916 }
1917 
1918 void writestring( int* fileDescriptor, const char inString[] )
1919 {
1920   int filePtr = *fileDescriptor - 1;
1921   FILE* fileObject = fileArray[filePtr] ;
1922   fprintf(fileObject,"%s",inString );
1923   return;
1924 }
1925 
1926 void Gather_Headers( int* fileDescriptor, vector< string >& headers )
1927 {
1928   FILE* fileObject;
1929   char Line[1024];
1930 
1931   fileObject = fileArray[ (*fileDescriptor)-1 ];
1932 
1933   while( !feof(fileObject) ) {
1934     fgets( Line, 1024, fileObject);
1935     if ( Line[0] == '#' ) {
1936       headers.push_back( Line );
1937     } else {
1938       break;
1939     }
1940   }
1941   rewind( fileObject );
1942   clearerr( fileObject );
1943 }
1944 
1945 void isWrong( void ) {
1946   (Wrong_Endian) ? fprintf(stdout,"YES\n") : fprintf(stdout,"NO\n");
1947 }
1948 
1949 void togglestrictmode( void ) { Strict_Error = !Strict_Error; }
1950 
1951 int isLittleEndian( void )
1952 {
1953   // this function returns a 1 if the current running architecture is
1954   // LittleEndian Byte Ordered, else it returns a zero
1955   union  {
1956     long a;
1957     char c[sizeof( long )];
1958   } endianUnion;
1959   endianUnion.a = 1 ;
1960   if ( endianUnion.c[sizeof(long)-1] != 1 ) return 1 ;
1961   else return 0;
1962 }
1963 
1964 namespace PHASTA {
1965   const char* const PhastaIO_traits<int>::type_string = "integer";
1966   const char* const PhastaIO_traits<double>::type_string = "double";
1967 }
1968