Commit 249de949 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Impemented reading of C binary files with doubles

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@123 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent c6985522
......@@ -8,11 +8,15 @@ MODULE mpdef
! precision constants
INTRINSIC :: selected_real_kind
INTRINSIC :: selected_int_kind
INTEGER, PARAMETER :: mpi = selected_int_kind(9) !> 4 byte integer
INTEGER, PARAMETER :: mpl = selected_int_kind(18) !> 8 byte integer
INTEGER, PARAMETER :: mps = selected_real_kind(6, 37) !> 4 byte float
INTEGER, PARAMETER :: mpd = selected_real_kind(15, 307) !> 8 byte float
INTEGER, PARAMETER :: mpq = selected_real_kind(33, 4931) !> 16 byte float, gcc needs libquadmath
INTEGER, PARAMETER :: mpi4 = selected_int_kind(9) !> 4 byte integer
INTEGER, PARAMETER :: mpi8 = selected_int_kind(18) !> 8 byte integer
INTEGER, PARAMETER :: mpr4 = selected_real_kind(6, 37) !> 4 byte float
INTEGER, PARAMETER :: mpr8 = selected_real_kind(15, 307) !> 8 byte float
INTEGER, PARAMETER :: mpr16 = selected_real_kind(33, 4931) !> 16 byte float, gcc needs libquadmath INTEGER, PARAMETER :: mpi = selected_int_kind(9) !> 4 byte integer
INTEGER, PARAMETER :: mpi = mpi4 !> integer
INTEGER, PARAMETER :: mpl = mpi8 !> long integer
INTEGER, PARAMETER :: mps = mpr4 !> single precision
INTEGER, PARAMETER :: mpd = mpr8 !> double precision
!> list items from steering file
TYPE listItem
INTEGER(mpi) :: label
......
......@@ -86,7 +86,8 @@ MODULE mpmod
INTEGER(mpi) :: nagbn !< max number of global paramters per record
INTEGER(mpi) :: nalcn !< max number of local paramters per record
INTEGER(mpi) :: naeqn !< max number of equations (measurements) per record
INTEGER(mpi) :: nrec !< (current) record number
INTEGER(mpi) :: nrec !< number of records read
INTEGER(mpi) :: nrecd !< number of records read containing doubles
REAL(mps) :: dflim !< convergence limit
INTEGER(mpi), DIMENSION(0:3) :: nrejec !< rejected events
REAL(mps), DIMENSION(0:8) :: times !< cpu time counters
......@@ -172,7 +173,8 @@ MODULE mpmod
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: readBufferInfo !< buffer management (per thread)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: readBufferPointer !< pointer to used buffers
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: readBufferDataI !< integer data
REAL(mps), DIMENSION(:), ALLOCATABLE :: readBufferDataF !< float data
REAL(mpr4), DIMENSION(:), ALLOCATABLE :: readBufferDataF !< float data
REAL(mpr8), DIMENSION(:), ALLOCATABLE :: readBufferDataD !< double data
! global parameter usage in record
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: globalIndexUsage !< indices of global par in record
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: backIndexUsage !< list of global par in record
......@@ -188,7 +190,7 @@ MODULE mpmod
REAL(mpd), DIMENSION(:), ALLOCATABLE::vzru !< local fit 'border solution'
REAL(mpd), DIMENSION(:), ALLOCATABLE::scdiag !< local fit workspace (D)
INTEGER(mpi), DIMENSION(:), ALLOCATABLE:: scflag !< local fit workspace (I)
REAL(mps), DIMENSION(:), ALLOCATABLE :: localCorrections !< local fit corrections (to residuals)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: localCorrections !< local fit corrections (to residuals)
REAL(mpd), DIMENSION(:), ALLOCATABLE :: localGlobalMatrix !< matrix correlating local and global par
! update of global matrix
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: writeBufferInfo !< write buffer management (per thread)
......
This diff is collapsed.
/**
* \file
*
......@@ -22,9 +21,13 @@
* - update on July 14th, 2008
* - update on October 29th, 2008: return for file number in \c readC()
*
* Last major update on April 24th, 2012 by C.Kleinwort:
* Major updates on April 24th, 2012 by C.Kleinwort:
* - skip records larger than buffer size (to determine max record length)
* - dynamic allocation of file pointer list (no hard-coded max number of files)
*
* Last major update on February 26th, 2014 by C.Kleinwort:
* - implement reading of records containing doubles (instead of floats)
* indicated by negative record length.
*/
#ifdef USE_SHIFT_RFIO
......@@ -52,30 +55,28 @@ FILE **files; ///< pointer to list of pointers to opened binary files
unsigned int maxNumFiles; ///< max number of files
unsigned int numAllFiles; ///< number of opened files
/*______________________________________________________________*/
/// Initialises the 'global' variables used for file handling.
/**
* \param[in] nFiles Maximal number of C binary files to use.
*/
void initC(int *nFiles)
{
maxNumFiles = *nFiles;
void initC(int *nFiles) {
maxNumFiles = *nFiles;
#ifdef USE_ZLIB
printf(" initC: using zlib version %s\n",ZLIB_VERSION);
files = (gzFile **) malloc(sizeof(gzFile *)*maxNumFiles);
printf(" initC: using zlib version %s\n",ZLIB_VERSION);
files = (gzFile **) malloc(sizeof(gzFile *)*maxNumFiles);
#else
files = (FILE **) malloc(sizeof(FILE *)*maxNumFiles);
files = (FILE **) malloc(sizeof(FILE *) * maxNumFiles);
#endif
{
int i = 0;
for ( ; i < maxNumFiles; ++i) {
files[i] = 0;
}
}
numAllFiles = 0;
{
int i = 0;
for (; i < maxNumFiles; ++i) {
files[i] = 0;
}
}
numAllFiles = 0;
}
FCALLSCSUB1(initC,INITC,initc,PINT)
FCALLSCSUB1( initC, INITC, initc, PINT)
/*______________________________________________________________*/
......@@ -94,19 +95,19 @@ FCALLSCSUB1(initC,INITC,initc,PINT)
/**
* \param[in] nFileIn File number (1 .. maxNumFiles)
*/
void resetC(int *nFileIn)
{
int fileIndex = *nFileIn-1; /* index of current file */
if (fileIndex < 0) return; /* no file opened at all... */
void resetC(int *nFileIn) {
int fileIndex = *nFileIn - 1; /* index of current file */
if (fileIndex < 0)
return; /* no file opened at all... */
#ifdef USE_ZLIB
gzrewind(files[fileIndex]);
gzrewind(files[fileIndex]);
#else
/* rewind(files[fileIndex]); Does not work with rfio, so call: */
fseek(files[fileIndex], 0L, SEEK_SET);
clearerr(files[fileIndex]); /* These two should be the same as rewind... */
/* rewind(files[fileIndex]); Does not work with rfio, so call: */
fseek(files[fileIndex], 0L, SEEK_SET);
clearerr(files[fileIndex]); /* These two should be the same as rewind... */
#endif
}
FCALLSCSUB1(resetC,RESETC,resetc,PINT)
FCALLSCSUB1( resetC, RESETC, resetc, PINT)
/*______________________________________________________________*/
/// Open file.
......@@ -120,39 +121,41 @@ void openC(const char *fileName, int *errorFlag)
* * 3: if file opened, but with error (can that happen?)
*/
{
/* No return value since to be called as subroutine from fortran */
/* No return value since to be called as subroutine from fortran */
if (!errorFlag) return; /* 'printout' error? */
if (!errorFlag)
return; /* 'printout' error? */
if (numAllFiles >= maxNumFiles) {
*errorFlag = 1;
} else {
if (numAllFiles >= maxNumFiles) {
*errorFlag = 1;
} else {
#ifdef USE_ZLIB
files[numAllFiles] = gzopen(fileName, "rb");
if (!files[numAllFiles]) {
*errorFlag = 2;
} else
files[numAllFiles] = gzopen(fileName, "rb");
if (!files[numAllFiles]) {
*errorFlag = 2;
} else
#else
files[numAllFiles] = fopen(fileName, "rb");
if (!files[numAllFiles]) {
*errorFlag = 2;
} else if (ferror(files[numAllFiles])) {
fclose(files[numAllFiles]);
files[numAllFiles] = 0;
*errorFlag = 3;
} else
files[numAllFiles] = fopen(fileName, "rb");
if (!files[numAllFiles]) {
*errorFlag = 2;
} else if (ferror(files[numAllFiles])) {
fclose(files[numAllFiles]);
files[numAllFiles] = 0;
*errorFlag = 3;
} else
#endif
{
++numAllFiles; /* We have one more opened file! */
*errorFlag = 0;
}
}
{
++numAllFiles; /* We have one more opened file! */
*errorFlag = 0;
}
}
}
FCALLSCSUB2(openC,OPENC,openc,STRING,PINT)
FCALLSCSUB2( openC, OPENC, openc, STRING, PINT)
/*______________________________________________________________*/
/// Read record from file.
/**
* \param[out] bufferDouble read buffer for doubles
* \param[out] bufferFloat read buffer for floats
* \param[out] bufferInt read buffer for integers
* \param[in,out] lengthBuffers in: buffer length, out: number of floats/ints in records
......@@ -164,168 +167,227 @@ FCALLSCSUB2(openC,OPENC,openc,STRING,PINT)
* * -4: given buffers too short for record
* * -8: problem with stream or EOF reading floats
* * -16: problem with stream or EOF reading ints
* * -32: problem with stream or EOF reading doubles
* * =0: reached end of file (or read empty record?!)
* * >0: number of words (floats + integers) read and stored in buffers
* * =4: found floats
* * =8: found doubles
*/
void readC(float *bufferFloat, int *bufferInt, int *lengthBuffers,
int *nFileIn, int *errorFlag)
{
/* No return value since to be called as subroutine from fortran,
negative *errorFlag are errors, otherwise fine.
*nFileIn: number of the file the record is read from,
starting from 1 (not 0)
*/
void readC(double *bufferDouble, float *bufferFloat, int *bufferInt,
int *lengthBuffers, int *nFileIn, int *errorFlag) {
/* No return value since to be called as subroutine from fortran,
negative *errorFlag are errors, otherwise fine.
if (!errorFlag) return;
*errorFlag = 0;
int fileIndex = *nFileIn-1; /* index of current file */
if (fileIndex < 0) return; /* no file opened at all... */
if (!bufferFloat || !bufferInt || !lengthBuffers) {
*errorFlag = -1;
return;
}
*nFileIn: number of the file the record is read from,
starting from 1 (not 0)
*/
int doublePrec = 0;
/* read length of 'record' */
int recordLength = 0; /* becomes number of words following in file */
if (!errorFlag)
return;
*errorFlag = 0;
int fileIndex = *nFileIn - 1; /* index of current file */
if (fileIndex < 0)
return; /* no file opened at all... */
if (!bufferFloat || !bufferInt || !lengthBuffers) {
*errorFlag = -1;
return;
}
/* read length of 'record' */
int recordLength = 0; /* becomes number of words following in file */
#ifdef USE_ZLIB
int nCheckR = gzread(files[fileIndex], &recordLength, sizeof(recordLength));
if (gzeof(files[fileIndex])) {
gzrewind(files[fileIndex]);
*errorFlag = 0; /* Means EOF of file. */
return;
}
int nCheckR = gzread(files[fileIndex], &recordLength, sizeof(recordLength));
if (gzeof(files[fileIndex])) {
gzrewind(files[fileIndex]);
*errorFlag = 0; /* Means EOF of file. */
return;
}
if (recordLength<0) {
doublePrec = 1;
recordLength = -recordLength;
}
if (sizeof(recordLength) != nCheckR) {
printf("readC: problem reading length of record file %d\n", fileIndex);
*errorFlag = -2;
return;
}
if (sizeof(recordLength) != nCheckR) {
printf("readC: problem reading length of record file %d\n", fileIndex);
*errorFlag = -2;
return;
}
if (recordLength/2 > *lengthBuffers) {
/* printf("readC: given buffers too short (%d, need > %d)\n", *lengthBuffers,
recordLength/2); */
/* skip floats */
int i=0;
if (doublePrec) {
for (; i< recordLength/2; ++i)
{
int nCheckD = gzread(files[fileIndex], bufferDouble, sizeof(bufferDouble[0]));
if (nCheckD != sizeof(bufferDouble[0])) {
printf("readC: problem with stream or EOF skipping doubles\n");
*errorFlag = -32;
return;
}
}
} else {
for (; i< recordLength/2; ++i)
{
int nCheckF = gzread(files[fileIndex], bufferFloat, sizeof(bufferFloat[0]));
if (nCheckF != sizeof(bufferFloat[0])) {
printf("readC: problem with stream or EOF skipping floats\n");
*errorFlag = -8;
return;
}
}
}
i=0;
/* skip ints */
for (; i< recordLength/2; ++i)
{
int nCheckI = gzread(files[fileIndex], bufferInt, sizeof(bufferInt[0]));
if (nCheckI != sizeof(bufferInt[0])) {
printf("readC: problem with stream or EOF skipping ints\n");
*errorFlag = -16;
return;
}
}
if (recordLength/2 > *lengthBuffers) {
/* printf("readC: given buffers too short (%d, need > %d)\n", *lengthBuffers,
recordLength/2); */
/* skip floats */
int i=0;
for ( ; i< recordLength/2; ++i)
{
int nCheckF = gzread(files[fileIndex], bufferFloat, sizeof(bufferFloat[0]));
if (nCheckF != sizeof(bufferFloat[0])) {
printf("readC: problem with stream or EOF skiping floats\n");
*errorFlag = -8;
return;
}
}
i=0;
/* skip ints */
for ( ; i< recordLength/2; ++i)
{
int nCheckI = gzread(files[fileIndex], bufferInt, sizeof(bufferInt[0]));
if (nCheckI != sizeof(bufferInt[0])) {
printf("readC: problem with stream or EOF skiping ints\n");
*errorFlag = -16;
return;
}
}
*errorFlag = -4;
*lengthBuffers = recordLength/2;
return;
} else {
*lengthBuffers = recordLength/2;
}
*errorFlag = -4;
*lengthBuffers = recordLength/2;
return;
} else {
*lengthBuffers = recordLength/2;
}
/* read floats (i.e. derivatives + value + sigma) */
int nCheckF = gzread(files[fileIndex], bufferFloat, *lengthBuffers*4);
if (nCheckF != *lengthBuffers*4) {
printf("readC: problem with stream or EOF reading floats\n");
*errorFlag = -8;
return;
}
/* read floats (i.e. derivatives + value + sigma) */
if (doublePrec) {
int nCheckD = gzread(files[fileIndex], bufferDouble, *lengthBuffers*8);
if (nCheckD != *lengthBuffers*8) {
printf("readC: problem with stream or EOF reading doubles\n");
*errorFlag = -32;
return;
}
} else {
int nCheckF = gzread(files[fileIndex], bufferFloat, *lengthBuffers*4);
if (nCheckF != *lengthBuffers*4) {
printf("readC: problem with stream or EOF reading floats\n");
*errorFlag = -8;
return;
}
int i=0;
for (; i< recordLength/2; ++i) bufferDouble[i] = (double) bufferFloat[i];
}
/* read ints (i.e. parameter lables) */
int nCheckI = gzread(files[fileIndex], bufferInt, *lengthBuffers*4);
if (nCheckI != *lengthBuffers*4) {
printf("readC: problem with stream or EOF reading ints\n");
*errorFlag = -16;
return;
}
/* read ints (i.e. parameter labels) */
int nCheckI = gzread(files[fileIndex], bufferInt, *lengthBuffers*4);
if (nCheckI != *lengthBuffers*4) {
printf("readC: problem with stream or EOF reading ints\n");
*errorFlag = -16;
return;
}
#else
size_t nCheckR = fread(&recordLength, sizeof(recordLength), 1,
files[fileIndex]);
if (feof(files[fileIndex])) {
/* rewind(files[fileIndex]); Does not work with rfio, so call: */
fseek(files[fileIndex], 0L, SEEK_SET);
clearerr(files[fileIndex]); /* These two should be the same as rewind... */
*errorFlag = 0; /* Means EOF of file. */
return;
}
size_t nCheckR = fread(&recordLength, sizeof(recordLength), 1,
files[fileIndex]);
if (feof(files[fileIndex])) {
/* rewind(files[fileIndex]); Does not work with rfio, so call: */
fseek(files[fileIndex], 0L, SEEK_SET);
clearerr(files[fileIndex]); /* These two should be the same as rewind... */
*errorFlag = 0; /* Means EOF of file. */
return;
}
if (1 != nCheckR || ferror(files[fileIndex])) {
printf("readC: problem reading length of record, file %d\n",
fileIndex);
*errorFlag = -2;
return;
}
if (1 != nCheckR || ferror(files[fileIndex])) {
printf("readC: problem reading length of record, file %d\n", fileIndex);
*errorFlag = -2;
return;
}
if (recordLength/2 > *lengthBuffers) {
/* printf("readC: given buffers too short (%d, need > %d)\n", *lengthBuffers,
recordLength/2); */
/* skip floats */
int i=0;
for ( ; i< recordLength/2; ++i)
{
size_t nCheckF = fread(bufferFloat, sizeof(bufferFloat[0]), 1,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckF != *lengthBuffers) {
printf("readC: problem with stream or EOF skiping floats\n");
*errorFlag = -8;
return;
}
}
i=0;
/* skip ints */
for ( ; i< recordLength/2; ++i)
{
size_t nCheckI = fread(bufferInt, sizeof(bufferInt[0]), 1,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckI != *lengthBuffers) {
printf("readC: problem with stream or EOF skiping ints\n");
*errorFlag = -16;
return;
}
}
*errorFlag = -4;
*lengthBuffers = recordLength/2;
return;
} else {
*lengthBuffers = recordLength/2;
}
if (recordLength < 0) {
doublePrec = 1;
recordLength = -recordLength;
}
if (recordLength / 2 > *lengthBuffers) {
/* printf("readC: given buffers too short (%d, need > %d)\n", *lengthBuffers,
recordLength/2); */
/* skip floats */
int i = 0;
if (doublePrec) {
for (; i < recordLength / 2; ++i) {
size_t nCheckD = fread(bufferDouble, sizeof(bufferDouble[0]), 1,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckD != *lengthBuffers) {
printf(
"readC: problem with stream or EOF skipping doubles\n");
*errorFlag = -32;
return;
}
}
} else {
for (; i < recordLength / 2; ++i) {
size_t nCheckF = fread(bufferFloat, sizeof(bufferFloat[0]), 1,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckF != *lengthBuffers) {
printf(
"readC: problem with stream or EOF skipping floats\n");
*errorFlag = -8;
return;
}
}
}
i = 0;
/* skip ints */
for (; i < recordLength / 2; ++i) {
size_t nCheckI = fread(bufferInt, sizeof(bufferInt[0]), 1,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckI != *lengthBuffers) {
printf("readC: problem with stream or EOF skiping ints\n");
*errorFlag = -16;
return;
}
}
/* read floats (i.e. derivatives + value + sigma) */
size_t nCheckF = fread(bufferFloat, sizeof(bufferFloat[0]), *lengthBuffers,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckF != *lengthBuffers) {
printf("readC: problem with stream or EOF reading floats\n");
*errorFlag = -8;
return;
}
*errorFlag = -4;
*lengthBuffers = recordLength / 2;
return;
} else {
*lengthBuffers = recordLength / 2;
}
/* read ints (i.e. parameter lables) */
size_t nCheckI = fread(bufferInt, sizeof(bufferInt[0]), *lengthBuffers,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckI != *lengthBuffers) {
printf("readC: problem with stream or EOF reading ints\n");
*errorFlag = -16;
return;
}
/* read floats (i.e. derivatives + value + sigma) */
if (doublePrec) {
size_t nCheckD = fread(bufferDouble, sizeof(bufferDouble[0]),
*lengthBuffers, files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckD != *lengthBuffers) {
printf("readC: problem with stream or EOF reading doubles\n");
*errorFlag = -32;
return;
}
} else {
size_t nCheckF = fread(bufferFloat, sizeof(bufferFloat[0]),
*lengthBuffers, files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckF != *lengthBuffers) {
printf("readC: problem with stream or EOF reading floats\n");
*errorFlag = -8;
return;
}
int i = 0;
for (; i < recordLength / 2; ++i)
bufferDouble[i] = (double) bufferFloat[i];
}
/* read ints (i.e. parameter labels) */
size_t nCheckI = fread(bufferInt, sizeof(bufferInt[0]), *lengthBuffers,
files[fileIndex]);
if (ferror(files[fileIndex]) || feof(files[fileIndex])
|| nCheckI != *lengthBuffers) {
printf("readC: problem with stream or EOF reading ints\n");
*errorFlag = -16;
return;
}
#endif
*errorFlag = *lengthBuffers;
}
FCALLSCSUB5(readC,READC,readc,PFLOAT,PINT,PINT,PINT,PINT)
*errorFlag = 4 * (doublePrec + 1);
}
FCALLSCSUB6(readC,READC,readc,PDOUBLE,PFLOAT,PINT,PINT,PINT,PINT)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment