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