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