Reading 16 byte Fortran floats into Python from a file

Question

I am writing sequences of floating points from Fortran into binary files and want to read them in Python. All works fine with single and double precision floats (kind=4) and (kind=8), but when I attempt to move to a real(kind=16) variable type, suddenly things no longer work (data array filled with zeros). I read here: python read 16 bytes long double from binary file That a workaround is needed for the np.fromfile function. I implemented the suggested change but am still not getting the right result. Mwe Python and Fortran code is below. I tried this both with Python 2.7 + Numpy 1.8 and Python 3.4 + Numpy 1.14 to the same effect. I also checked and the generated file seems to have the right amount of data (480 bytes for 30 floats with 16 bytes). Any help will be welcome!

Python reader:

import numpy as np inputfilename = "fortranData.bin" dp = 8 nVals = 30 with open(inputfilename, 'rb') as f: # both work fine for 4 or 8 byte floats (32 and 64) # but not for 16 byte floats (np.float128) data = np.fromfile(f, dtype=np.float64) data = np.frombuffer(f.read(dp*nVals), dtype=np.float64, count=nVals) print(data)

Fortran writer (compiled with gfortran 4.8.6 on Ubuntu 14.04)

program test implicit none integer, parameter :: dp=8 integer, parameter :: outFileUnit=51, nPts=10, nDim=3 character(len=100) :: binaryFile="fortranData.bin" real(kind=dp), dimension(nPts,nDim) :: data integer :: i, j ! - generate data - do i=1, nPts do j=1, nDim data(i,j) = (i-1)*nDim + j enddo enddo ! - open file and write data - open(unit=outFileUnit, file=binaryFile, form='unformatted', access='direct', & status='replace', recl=nDim*sizeof(data(1,1))) do i=1, nPts write(outFileUnit, rec=i) data(i,:) enddo end program test

EDIT: as requested, here's a few lines of the data at the top of the fortranData.bin file:

$od fortranData.bin 0000000 000000 000000 000000 000000 000000 000000 000000 037777 0000020 000000 000000 000000 000000 000000 000000 000000 040000 0000040 000000 000000 000000 000000 000000 000000 100000 040000 0000060 000000 000000 000000 000000 000000 000000 000000 040001 0000100 000000 000000 000000 000000 000000 000000 040000 040001 0000120 000000 000000 000000 000000 000000 000000 100000 040001

Could you show 16 bytes of the file as test data so we can try it? — user555045
– user555045, Commented Sep 29, 2018 at 11:54
kind=a_number is not always the number of bytes. See stackoverflow.com/questions/838310/fortran-90-kind-parameter — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Sep 29, 2018 at 20:40
Right, thanks for pointing that out, I wasn't aware of different compilers treating it differently, I've only been testing this with gfortran. Still, since I'm using sizeof(data(1,1), adapting the code to a different compiler would just be a matter of adjusting the dp parameter I guess. Still, better practice would be to use the REAL64 etc. types from the new standard pointed out in the link you mentioned. Thanks! — Artur
– Artur, Commented Sep 30, 2018 at 7:49

user555045 · Accepted Answer · 2018-09-29 12:41:20Z

It seems that numpy does not support the IEEE Quad format, so far as I know it just needs to be converted manually. For example, if you read the file in chunks of 16 bytes, then the chunks can be converted like this (poorly tested though)

def rawQuadToDouble(raw): asint = int.from_bytes(raw, byteorder='little') sign = (-1.0) ** (asint >> 127); exponent = ((asint >> 112) & 0x7FFF) - 16383; significand = (asint & ((1 << 112) - 1)) | (1 << 112) return sign * significand * 2.0 ** (exponent - 112)

Fantastic, works great on the test piece of data. Thanks a lot!

Collectives™ on Stack Overflow

Reading 16 byte Fortran floats into Python from a file

1 Answer 1

1 Comment

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Linked

Related