Read and process a Fortran90 binary file

Hi,

I have a Fortran90 code where a variable called 'stream' is defined such that it stores certain variable every nstep:

open(19,file='stream',form='UNFORMATTED',position='append')
if(mod(NSTEP,100).eq.0)then
    write(19)NSTEP,U,W,zeta
    endif

How can I properly read the file in Matlab for postprocessing?

FYI, the fortran code is built with double precision.

Thank you

[James Tursa]

Unformatted Fortran can contain beginning and end of line markers in the file. You need to give us a small sample file to work with to verify this, and also tell us the type and dimensions of the variables NSTEP, U, W, and zeta.

[Walter Roberson]

The FORTRAN runtime system embeds the record boundaries in the data by inserting an INTEGER*4 byte count at the beginning and end of each unformatted sequential record during an unformatted sequential WRITE. The trailing byte count enables BACKSPACE to operate on records.

https://docs.oracle.com/cd/E19957-01/805-4939/6j4m0vnc4/index.html

[Carola Forlini]

I am attaching an extract of the file.

NSTEP is an integer = 2000

U, W, zeta are matrices = 32x512

Thank you

Answer 1

[Les Beckham]

It would help to have a sample file to experiment with.

Nevertheless, assuming that the write call that you have shown is the only thing that writes to this file in your Fortran code, I would suggest something like this using fread and reshape:

fid = fopen('your_file_name.ext', 'rb');
A = fread(fid);
data = reshape(A, 4, numel(A)/4)';
nstep = A(:,1);
U = A(:,2);
W = A(:,3);
zeta = A(:,4);

Note that you might run into endian-ness issues if your Fortran code runs on a different architecture than your Matlab is running on. See the documentation for the machinefmt option to fread if necessary.

[Walter Roberson]

form='UNFORMATTED' has record length markers before and after each record.

[Walter Roberson]

If you start with fread(fid, 1, 'integer*4') and discard that, and thereafter fread() with a skip field of 8 then it should be properly positioned to fread() more.

Unfortunately we are not told the size of NSTEP, U, W, zeta

[Carola Forlini]

Hi, I provided a sample file and the dimensions of NSTEP, U, W, zeta in the post above.

Thank you.

[Les Beckham]

Thanks for adding the example.

Can you clarify what you mean by this?

NSTEP is an integer = 2000

Especially the "= 2000" part. Is that just an example of a possible value? What is the size of this integer (number of bits)?

Also, are U, W, and zeta double precision (64 bit) floating point 2-d arrays with the given size?

[Carola Forlini]

I apologize for the misunderstanding.The size of NSTEP is 32 bit.

U, W, zeta are single precision (32 bit) of size N1×N3, with N1=32 and N3=512.

Answer 2

[Les Beckham]

In light of the additional variable size information and sample file provided by the OP, and the information from Walter about Fortran inserting record boundary byte counts at the beginning and end of each record written with an unformatted sequential WRITE, here is a fresh answer.

Note that this makes an assumption about the order the data is written (column major order). If this assumption is not correct, the inner for loop will need to be modified accordingly.

Hopefully this is close enough to get you where you need to go.

filenames = unzip('stream.zip')

filenames = 1x1 cell array
    {'stream'}

fid = fopen(filenames{1}, 'rb');

% start/stop byte counts (4 bytes each) + NSTEP (4 bytes) + U, W, zeta (3*32*512*4 bytes)
recordSize = 2*4 + 4 + 3*32*512*4; 

fseek(fid, 0, 'eof'); % seek to the end of the file
nBytes = ftell(fid); % find number of bytes in the file
frewind(fid); % go back to the beginning of the file

nRecords = nBytes / recordSize;

% preallocate
nStep = zeros(nRecords, 1);      % column vector 
U     = zeros(3, 512, nRecords); % 3D arrays
W     = zeros(3, 512, nRecords);
zeta  = zeros(3, 512, nRecords);

for iRecord = 1:nRecords
    fseek(fid, 4, 'cof'); % skip begin record byte count field
    nStep(iRecord) = fread(fid, 1, 'integer*4');
    for iRow = 1:3
        U(iRow, :, iRecord)    = fread(fid, 512, 'real*4');
        W(iRow, :, iRecord)    = fread(fid, 512, 'real*4'); 
        zeta(iRow, :, iRecord) = fread(fid, 512, 'real*4');
    end
    fseek(fid, 4, 'cof'); % skip end record byte count field
end

fclose(fid); % close the file

whos nStep U W zeta % check variable sizes

  Name        Size                 Bytes  Class     Attributes

  U           3x512x21            258048  double              
  W           3x512x21            258048  double              
  nStep      21x1                    168  double              
  zeta        3x512x21            258048  double              

plot(nStep); % make a plot see if nStep makes sense
grid on

Well, I don't know what your NSTEP should look like, but I would have expected a monotonically increasing count. So, either I made a mistake, or there is some additional information needed. It does appear that the size information you gave is either correct, or just coincidentally came out to be an integer divisor of the file size.

[Voss]

Reasonable approach, but I'm not sure why you assume 3-by-512 matrices and read them one row at a time. If instead 32x512 is used and each read at once, then the nStep comes out correct (i.e., [0 100 200 ... 2000]).

filenames = unzip('stream.zip')

filenames = 1x1 cell array
    {'stream'}

fid = fopen(filenames{1}, 'rb');

% start/stop byte counts (4 bytes each) + NSTEP (4 bytes) + U, W, zeta (3*32*512*4 bytes)
recordSize = 2*4 + 4 + 3*32*512*4;

fseek(fid, 0, 'eof'); % seek to the end of the file
nBytes = ftell(fid); % find number of bytes in the file
frewind(fid); % go back to the beginning of the file

nRecords = nBytes / recordSize;

% preallocate
nStep = zeros(nRecords, 1);       % column vector 
U     = zeros(32, 512, nRecords); % 3D arrays
W     = zeros(32, 512, nRecords);
zeta  = zeros(32, 512, nRecords);

for iRecord = 1:nRecords
    fseek(fid, 4, 'cof'); % skip begin record byte count field
    nStep(iRecord) = fread(fid, 1, 'uint32');
    U(:, :, iRecord)    = fread(fid, [32 512], 'single');
    W(:, :, iRecord)    = fread(fid, [32 512], 'single');
    zeta(:, :, iRecord) = fread(fid, [32 512], 'single');
    fseek(fid, 4, 'cof'); % skip end record byte count field
end

% make sure the entire file was read
assert(isequal(ftell(fid),nBytes))

fclose(fid); % close the file

whos nStep U W zeta % check variable sizes

  Name        Size                  Bytes  Class     Attributes

  U          32x512x21            2752512  double              
  W          32x512x21            2752512  double              
  nStep      21x1                     168  double              
  zeta       32x512x21            2752512  double              

plot(nStep,'-o'); % make a plot see if nStep makes sense
grid on

[Carola Forlini]

Thank you for your reply! this is exactly how I was trying to exctract the data!

[Les Beckham]

Good catch, Voss. Thanks. I don't know why I switched from using 32 in the number of records calculation, to using 3 in the section where I did the preallocation and the reads. I should have spent more time reviewing my answer before submitting it, as it was obvious that something was wrong.

[James Tursa]

"... this makes an assumption about the order the data is written (column major order) ..."

Fortran is column-major storage order, same as MATLAB. When you specify just the variable name in a read/write statement, that is the element order for the read/write.