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