FORTRAN二进制文件的读写

! A fortran95 program for G95
! By WQY
program main
  implicit none
  integer::i,j,n
  parameter(n=4000)
    real*4::k(n)
  open(10,file='C:\\test\\vx.sgy',status='old',FORM='unformatted',access='direct',recl=4)
 ! rewind(10) 3200+400+240=3840,   1920*2 bytes,  960*4
do i=1,1
 do j=1,n
  read(10,rec=900+i*60+j) k(j)
 end do
end do
  close(10)
open(10,file='C:\\test\\vx.txt',status='replace',FORM='formatted',access='direct',recl=13)
!REWIND(10)
do j=1,n
   write(10,'(F10.2)',rec=j) k(j)
end do
  close(10)
end

得到的VX.TXT文件像是一张TXT表格。

program main
  implicit none
  integer::i,j,n
  parameter(n=4000)
    real*4::k(n)
  open(10,file='C:\\test\\vx.sgy',status='old',FORM='unformatted',access='direct',recl=4)
 ! rewind(10) 3200+400+240=3840,   1920*2 bytes,  960*4
do i=1,1
 do j=1,n
  read(10,rec=900+i*60+j) k(j)
 end do
end do
  close(10)
!open(10,file='C:\\test\\vx.txt',status='replace',FORM='formatted',access='direct',recl=13)
!REWIND(10)
!do j=1,n
!   write(10,'(F10.2)',rec=j) k(j) !one record one entrie
!end do
!close(10)

open(10,file='C:\\test\\vx.txt',status='replace',FORM='formatted',access='sequential',recl=13)
do j=1,n
   write(10,'(F10.2)') k(j)! one line one record
end do
close(10)
end

FORTRAN调用MATLAB的MAT文件

matlabroot/extern/examples/eng_mat 文件夹包含示例的 C/C++ 和 Fortran 源代码，用于说明如何使用 MAT 文件例程。这些示例创建独立程序。源代码对于 Windows®、macOS 和 Linux® 系统都是相同的。

要编译代码示例，请首先将文件复制到一个可写文件夹中，例如 Windows 路径上的 c:\work。

copyfile(fullfile(matlabroot,'extern','examples','eng_mat',...
'filename'), fullfile('c:','work'))

其中，filename 为源代码文件的名称。

MATLAB MEX 调用FORTRAN

Create Fortran Source MEX File

This example shows how to write a MEX file to call a Fortran subroutine, timestwo, in MATLAB® using a MATLAB matrix. You can view the complete source file here. This example uses the MATLAB Editor to write the source code and the MATLAB mex command to create the MEX function.

Fortran subroutine timestwo

The following code defines the timestwo subroutine, which multiplies an n-dimensional array, x_input, by 2, and returns the results in array, y_output.

      subroutine timestwo(y_output, x_input)
      real*8 x_input, y_output

      y_output = 2.0 * x_input
      return
      end

Create Source File

Open MATLAB Editor, create a file, and document the MEX file with the following information.

C======================================================================
C     timestwo.f
C     Computational function that takes a scalar and doubles it.
C     This is a MEX file for MATLAB.
C======================================================================

Add the Fortran header file, fintrf.h, containing the MATLAB API function declarations.

#include "fintrf.h"

Save the file on your MATLAB path, for example, in c:\work, and name it timestwo.F. The name of your MEX file is timestwo.

Create Gateway Routine

MATLAB uses the gateway routine, mexfunction, as the entry point to a Fortran subroutine. Add the following mexFunction code.

C     Gateway routine
      subroutine mexFunction(nlhs, plhs, nrhs, prhs)

C     Declarations

C     Statements

      return
      end

Add the following statement to your mexfunction subroutine to force you to declare all variables.

implicit none

Explicit type declaration is necessary for 64-bit arrays.

Declare mexfunction Arguments

To declare mxArray variables, use the MATLAB type, mwPointer. Add this code after the Declarations statement.

C     mexFunction arguments:
      mwPointer plhs(*), prhs(*)
      integer nlhs, nrhs

Declare Functions and Local Variables

• Declare the symbolic names and types of MATLAB API functions used in this MEX file.

  C     Function declarations:
        mwPointer mxGetDoubles
        mwPointer mxCreateDoubleMatrix
        integer mxIsNumeric
        mwPointer mxGetM, mxGetN

  To determine the type of a function, refer to the MATLAB API function reference documentation. For example, see the documentation for mxGetDoubles.

• Declare local variables for the mexfunction arguments.

  C     Pointers to input/output mxArrays:
        mwPointer x_ptr, y_ptr
  

• Declare matrix variables.

  C     Array information:
        mwPointer mrows, ncols
        mwSize size
  

Verify MEX File Input and Output Arguments

Verify the number of MEX file input and output arguments using the nrhs and nlhs arguments. Add these statements to the mexfunction code block.

C     Check for proper number of arguments. 
      if(nrhs .ne. 1) then
         call mexErrMsgIdAndTxt ('MATLAB:timestwo:nInput',
     +                           'One input required.')
      elseif(nlhs .gt. 1) then
         call mexErrMsgIdAndTxt ('MATLAB:timestwo:nOutput',
     +                           'Too many output arguments.')
      endif

Verify the input argument type using the prhs argument.

C     Check that the input is a number.
      if(mxIsNumeric(prhs(1)) .eq. 0) then
         call mexErrMsgIdAndTxt ('MATLAB:timestwo:NonNumeric',
     +                           'Input must be a number.')
      endif

Create Computational Routine

Add the timestwo code. This subroutine is your computational routine, the source code that performs the functionality you want to use in MATLAB.

C     Computational routine

      subroutine timestwo(y_output, x_input)
      real*8 x_input, y_output

      y_output = 2.0 * x_input
      return
      end

A computational routine is optional. Alternatively, you can place the code within the mexfunction function block.

Declare Variables for Computational Routine

Put the following variable declarations in mexFunction.

C     Arguments for computational routine:
      real*8  x_input, y_output

Read Input Array

To point to the input matrix data, use the mxGetDoubles function.

x_ptr = mxGetDoubles(prhs(1))

To create a Fortran array, x_input, use the mxCopyPtrToReal8 function.

C     Get the size of the input array.
      mrows = mxGetM(prhs(1))
      ncols = mxGetN(prhs(1))
      size = mrows*ncols

C     Create Fortran array from the input argument.
      call mxCopyPtrToReal8(x_ptr,x_input,size)

Prepare Output Data

To create the output argument, plhs(1), use the mxCreateDoubleMatrix function.

C     Create matrix for the return argument.
      plhs(1) = mxCreateDoubleMatrix(mrows,ncols,0)

Use the mxGetDoubles function to assign the y_ptr argument to plhs(1).

      y_ptr = mxGetDoubles(plhs(1))

Perform Calculation

Pass the arguments to timestwo.

C     Call the computational subroutine.
      call timestwo(y_output, x_input)

Copy Results to Output Argument

C     Load the data into y_ptr, which is the output to MATLAB.
      call mxCopyReal8ToPtr(y_output,y_ptr,size)     

View Complete Source File

Compare your source file with timestwo.F, located in the matlabroot/extern/examples/refbook folder. Open the file in the editor.

Build Binary MEX File

At the MATLAB command prompt, build the binary MEX file.

mex -R2018a timestwo.F

Test the MEX File

x = 99;
y = timestwo(x)

y =
   198

第 2 章 Fortran 输入/输出

本章介绍 Sun Studio Fortran 95 编译器提供的输入/输出功能。

2.1 从 Fortran 程序内部访问文件

数据通过 Fortran 逻辑单元在程序、设备或文件间进行传送。逻辑单元在 I/O 语句中用逻辑单元号来标识，逻辑单元号是从 0 到最大 4 字节整数值 (2,147,483,647) 的非负整数。

字符 * 可以作为逻辑单元标识符出现。当星号出现在 READ 语句中时，它代表标准输入文件；当它出现在 WRITE 或 PRINT 语句中时，则代表标准输出文件。

Fortran 逻辑单元可通过 OPEN 语句与特定的命名文件相关联。另外，在程序开始执行时，某些预连接单元会自动与特定文件相关联。

2.1.1 访问命名文件

OPEN 语句的 FILE= 说明符在运行时建立逻辑单元到命名物理文件的关联。该文件可以是预先就有的，也可以由程序创建。

OPEN 语句中的 FILE= 说明符可以指定一个简单文件名 (FILE=’myfile.out’)，也可以指定一个前面带有绝对或相对目录路径的文件名 (FILE=’../Amber/Qproj/myfile.out’)。另外，说明符还可以是字符常量、变量或字符表达式。

可以使用库例程将命令行参数和环境变量以字符变量形式送入程序，用作 OPEN 语句中的文件名。

以下示例 (GetFilNam.f) 展示了一种由键入的名称来构建绝对路径文件名的方法。程序使用库例程 GETENV、LNBLNK 和 GETCWD 返回 $HOME 环境变量的值，查找字符串中最后的非空格字符，确定当前工作目录：

CHARACTER F*128, FN*128, FULLNAME*128 PRINT*, ’ENTER FILE NAME:’ READ *, F FN = FULLNAME( F ) PRINT *, ’PATH IS: ’,FN END CHARACTER*128 FUNCTION FULLNAME( NAME ) CHARACTER NAME*(*), PREFIX*128 C This assumes C shell. C Leave absolute path names unchanged. C If name starts with ’~/’, replace tilde with home C directory; otherwise prefix relative path name with C path to current directory. IF ( NAME(1:1) .EQ. ’/’ ) THEN FULLNAME = NAME ELSE IF ( NAME(1:2) .EQ. ’~/’ ) THEN CALL GETENV( ’HOME’, PREFIX ) FULLNAME = PREFIX(:LNBLNK(PREFIX)) // 1 NAME(2:LNBLNK(NAME)) ELSE CALL GETCWD( PREFIX ) FULLNAME = PREFIX(:LNBLNK(PREFIX)) // 1 ’/’ // NAME(:LNBLNK(NAME)) ENDIF RETURN END

编译并运行 GetFilNam.f，结果如下：

demo% pwd /home/users/auser/subdir demo% f95 -o getfil GetFilNam.f demo% getfil ENTER FILE NAME: getfil PATH IS: /home/users/auser/subdir/atest.f demo%

2.1.4 向程序传递文件名中对这些例程进行了详细的介绍。有关详细信息，请参见 getarg(3F)、getcwd(3F) 和 getenv(3F) 的相应手册页条目；这些内容以及其他有用的库例程在《Fortran 库参考》中也有介绍。

2.1.2 不用文件名打开文件

OPEN 语句不需要指定名称；运行时系统会根据几个惯例提供文件名。

2.1.2.1 打开作为临时文件

在 OPEN 语句中指定 STATUS=’SCRATCH’ 会打开一个名称形式为 tmp.FAAAxnnnnn 的文件，其中 nnnnn 用当前进程 ID 替换，AAA 是一个包含三个字符的字符串，x 是一个字母；AAA 和 x 可确保文件名唯一。该文件在程序终止或执行 CLOSE 语句时被删除。在 FORTRAN 77 兼容模式 (-f77) 下编译时，可以在 CLOSE 语句中指定 STATUS=’KEEP’ 来保留这个临时文件。（此为非标准扩展。）

2.1.2.2 已打开

如果文件已被程序打开，可以使用后续的 OPEN 语句更改文件的某些特性；例如 BLANK 和 FORM。此时，只需指定文件的逻辑单元号以及要更改的参数。

2.1.2.3 预连接或隐式命名单元

程序执行开始时，会自动将三个单元号与特定的标准 I/O 文件相关联。这些预连接单元是标准输入、标准输出和标准错误：

• 标准输入是逻辑单元 5

• 标准输出是逻辑单元 6

• 标准错误是逻辑单元 0

通常，标准输入是从工作站键盘接收输入；标准输出和标准错误是在工作站屏幕上显示输出。

在其他所有情况下，如果在 OPEN 语句中指定了逻辑单元号而未在 FILE= 后指定任何名称，文件将以 fort.n 形式的名称打开，其中 n 为逻辑单元号。

2.1.3 不使用 OPEN 语句打开文件

在假定使用缺省惯例的情况下，并非必须使用 OPEN 语句。如果逻辑单元上的第一个操作是 I/O 语句，而不是 OPEN 或 INQUIRE，则会引用文件 fort.n，其中 n 为逻辑单元号（0、5 和 6 除外，它们有特殊意义）。

这些文件无需在程序执行前就存在。如果对文件的第一个操作不是 OPEN 或 INQUIRE 语句，则会创建这些文件。

示例：以下代码中，如果 WRITE 是该单元上的第一个输入/输出操作，则会创建文件 fort.25：

demo% cat TestUnit.f IU=25 WRITE( IU, ’(I4)’ ) IU END demo%

上述程序将打开文件 fort.25，并将一条格式化记录写入该文件：

demo% f95 -o testunit TestUnit.f demo% testunit demo% cat fort.25 25 demo%

2.1.4 向程序传递文件名

文件系统没有任何自动工具可将 Fortran 程序中的逻辑单元号与物理文件相关联。

但是，有几种令人满意的方式可将文件名传递给 Fortran 程序。

2.1.4.1 通过运行时参数和 GETARG

可以使用库例程 getarg(3F) 在运行时将命令行参数读入字符变量。参数会被解释为文件名并在 OPEN 语句的 FILE= 说明符中使用：

demo% cat testarg.f CHARACTER outfile*40 C Get first arg as output file name for unit 51 CALL getarg(1,outfile) OPEN(51,FILE=outfile) WRITE(51,*) ’Writing to file: ’, outfile END demo% f95 -o tstarg testarg.f demo% tstarg AnyFileName demo% cat AnyFileName Writing to file: AnyFileName demo%

2.1.4.2 通过环境变量和 GETENV

同样，可以使用库例程 getenv(3F) 在运行时将任何环境变量的值读入字符变量，该字符变量随后被解释为文件名：

demo% cat testenv.f CHARACTER outfile*40 C Get $OUTFILE as output file name for unit 51 CALL getenv(’OUTFILE’,outfile) OPEN(51,FILE=outfile) WRITE(51,*) ’Writing to file: ’, outfile END demo% f95 -o tstenv testenv.f demo% setenv OUTFILE EnvFileName demo% tstenv demo% cat EnvFileName Writing to file: EnvFileName demo%

使用 getarg 或 getenv 时，应该注意前导或尾随的空格。（Fortran 95 程序可以使用内函数 TRIM 或更早的 FORTRAN 77 库例程 LNBLNK()）在本章开头的示例中，可以随 FULLNAME 函数的代码行编写更加灵活的代码来接受相对路径名。

2.1.4.3 命令行 I/O 重定向和管道

将物理文件与程序的逻辑单元号相关联的另一种方式是重定向或管道输送预连接的标准 I/O 文件。重定向或管道在运行时执行命令中使用。

采用这种方式，读取标准输入（单元 5）和写入标准输出（单元 6）或标准错误（单元 0）的程序可以通过重定向（在命令行中使用 <、>、>>、>&、|、|&、2>、2>&1），读取或写入其他任何命名文件。

参见下表：

表 2–1 csh/sh/ksh 命令行重定向和管道

操作	使用 C Shell	使用 Bourne 或 Korn Shell
标准输入－从 mydata 读取	myprog < mydata	myprog < mydata
标准输出－写入（覆写）myoutput	myprog > myoutput	myprog > myoutput
标准输出－写入/追加至 myoutput	myprog >> myoutput	myprog >> myoutput
将标准错误重定向至文件	myprog >& errorfile	myprog 2> errorfile
将标准输出通过管道输送至另一程序的输入	myprog1 | myprog2	myprog1 | myprog2
将标准错误和输出通过管道输送至另一程序	myprog1 |& myprog2	myprog1 2>&1 | myprog2

有关命令行重定向和管道的详细信息，请参见 csh、ksh 和 sh 手册页。

2.2 直接 I/O

直接或随机 I/O 允许通过记录号直接访问文件。记录号在写入记录时分配。与顺序 I/O 不同，直接 I/O 记录可以按任何顺序读写。但是，在直接访问文件中，所有记录必须具有相同的固定长度。直接访问文件用文件的 OPEN 语句中的 ACCESS=’DIRECT’ 说明符声明。

直接访问文件中的逻辑记录是字节字符串，串长度由 OPEN 语句的 RECL= 说明符指定。READ 和 WRITE 语句指定的逻辑记录不能大于定义的记录大小。（记录大小以字节单位指定。）允许更短的记录。直接写入非格式化数据将使记录的未填写部分仍保持未定义。直接写入格式化数据将使未填写的记录用空格进行填充。

直接访问 READ 和 WRITE 语句另外还有一个参数 REC=n，用来指定要读取或写入的记录号。

示例：直接访问，非格式化：

OPEN( 2, FILE=’data.db’, ACCESS=’DIRECT’, RECL=200, & FORM=’UNFORMATTED’, ERR=90 ) READ( 2, REC=13, ERR=30 ) X, Y

本程序以直接访问、非格式化 I/O、记录固定长度为 200 字节的方式打开一个文件，然后将第十三条记录读入 X 和 Y。

示例：直接访问，格式化：

OPEN( 2, FILE=’inven.db’, ACCESS=’DIRECT’, RECL=200, & FORM=’FORMATTED’, ERR=90 ) READ( 2, FMT=’(I10,F10.3)’, REC=13, ERR=30 ) X, Y

本程序以直接访问、格式化 I/O、记录固定长度为 200 字节的方式打开一个文件。然后读取第十三条记录，并以 (I10,F10.3) 格式对其进行转换。

对于格式化文件，所写记录的大小由 FORMAT 语句确定。在上述示例中，FORMAT 语句所定义的记录大小为 20 个字符或字节。如果列表中的数据总量大于 FORMAT 语句中指定的记录大小，则可通过单条格式化写指令写入一条以上的记录。在此种情况下，会为随后的每一条记录赋予连续的记录号。

示例：直接访问、格式化、多记录写入：

OPEN( 21, ACCESS=’DIRECT’, RECL=200, FORM=’FORMATTED’) WRITE(21,’(10F10.3)’,REC=11) (X(J),J=1,100)

写入直接访问单元 21 的写指令会创建 10 条记录，每条记录 10 个元素（因为格式指定每条记录 10 个元素），这些记录从 11 到 20 进行编号。

2.3 二进制 I/O

Sun Studio Fortran 95 扩展了 OPEN 语句，允许声明“二进制”I/O 文件。

使用 FORM=’BINARY’ 打开文件与使用 FORM=’UNFORMATTED’ 具有大致相同的效果，所不同的是文件中没有嵌入记录长度。如果没有此数据，则无法知道一条记录的开始或结束位置。因此，无法对 FORM=’BINARY’ 文件执行 BACKSPACE 操作，这是因为不知道要退格到什么位置。对 ’BINARY’ 文件执行 READ 操作时，将按需要读取尽可能多的数据来填充输入列表中的变量。

• WRITE 语句：以二进制的形式将数据写入文件，并按输出列表中指定的数量传输字节。

• READ 语句：将数据读取到输入列表中的变量，并传输该列表所要求数量的字节。因为文件中没有记录标记，所以不进行“记录结束”错误检测。检测到的唯一错误是“文件结束”或异常系统错误。

• INQUIRE 语句：在使用 FORM="BINARY" 打开的文件中，INQUIRE 返回：FORM="BINARY" ACCESS="SEQUENTIAL" DIRECT="NO" FORMATTED="NO" UNFORMATTED="YES"。RECL= 和 NEXTREC= 没有定义。

• BACKSPACE 语句：不允许使用－返回一个错误。

• ENDFILE 语句：在当前位置照常截断文件。

• REWIND 语句：将文件照常重新定位到数据的开头。

2.4 流 I/O

f95 中实现了 Fortran 2003 标准“流”I/O 的新方案 。流 I/O 访问将数据文件视作连续的字节序列，用从 1 开始的正整数来寻址。可用 OPEN 语句中的 ACCESS=’STREAM’ 说明符来声明流 I/O 文件。字节地址文件定位要求 READ 或 WRITE 语句中有 POS=scalar_integer_expression 说明符。INQUIRE 语句接受 ACCESS=’STREAM’、说明符 STREAM=scalar_character_variable 和 POS=scalar_integer_variable。

流 I/O 在与 C 程序创建或读取的文件进行互操作时非常有用，如下例所示：

Fortran 95 program reads files created by C fwrite() program reader integer:: a(1024), i, result open(file="test", unit=8, access="stream",form="unformatted") ! read all of a read(8) a do i = 1,1024 if (a(i) .ne. i-1) print *,’error at ’, i enddo ! read the file backward do i = 1024,1,-1 read(8, pos=(i-1)*4+1) result if (result .ne. i-1) print *,’error at ’, i enddo close(8) end C program writes to a file #include <stdio.h> int binary_data[1024]; /* Create a file with 1024 32-bit integers */ int main(void) { int i; FILE *fp; for (i = 0; i < 1024; ++i) binary_data[i] = i; fp = fopen("test", "w"); fwrite(binary_data, sizeof(binary_data), 1, fp); fclose(fp); }

C 程序使用 C fwrite() 将 1024 个 32 位整数写入文件中。Fortran 95 读取程序以数组方式一次读取这些数据，然后再在文件中从后往前分别读取它们。第二条 read 语句中的 pos= 说明符说明位置是用字节表示的，从字节 1 开始（这一点与 C 相反，在 C 中，位置从字节 0 开始）。

2.5 内部文件

内部文件是 CHARACTER 类型的对象，如变量、子串、数组、数组元素或结构化记录的字段。内部文件 READ 可以来自常量字符串。内部文件 I/O 通过由一个字符对象向另一数据对象传送和转换数据，模拟格式化 READ 和 WRITE 语句。不执行任何文件 I/O。

使用内部文件时：

• 出现在 WRITE 语句中的是接收数据的字符对象的名称而非单元号。在 READ 语句中，出现的是字符对象源的名称而非单元号。

• 常量、变量或子串对象构成文件中的单条记录。

• 使用数组对象，每个数组元素对应于一条记录。

• 内部文件上的直接 I/O。（Fortran 95 标准只包括内部文件上的顺序格式化 I/O。）除了不能更改文件中的记录数之外，这一点与外部文件上的直接 I/O 相似。此时，记录是字符串数组的单个元素。这项非标准扩展仅在用 -f77 标志编译的 FORTRAN 77 兼容模式下可用。

• 每一顺序 READ 或 WRITE 语句均始于内部文件的开头。

示例：从内部文件（仅有一条记录）中以顺序、格式化方式进行读取：

demo% cat intern1.f CHARACTER X*80 READ( *, ’(A)’ ) X READ( X, ’(I3,I4)’ ) N1, N2 ! This codeline reads the internal file X WRITE( *, * ) N1, N2 END demo% f95 -o tstintern intern1.f demo% tstintern 12 99 12 99 demo%

示例：从内部文件（三条记录）中以顺序、格式化方式进行读取：

demo% cat intern2.f CHARACTER LINE(4)*16 DATA LINE(1) / ’ 81 81 ’ / DATA LINE(2) / ’ 82 82 ’ / DATA LINE(3) / ’ 83 83 ’ / DATA LINE(4) / ’ 84 84 ’ / READ( LINE,’(2I4)’) I,J,K,L,M,N PRINT *, I, J, K, L, M, N END demo% f95 intern2.f demo% a.out 81 81 82 82 83 83 demo%

示例：在 -f77 兼容模式下，从内部文件（一条记录）中以直接访问方式进行读取：

demo% cat intern3.f CHARACTER LINE(4)*16 DATA LINE(1) / ’ 81 81 ’ / DATA LINE(2) / ’ 82 82 ’ / DATA LINE(3) / ’ 83 83 ’ / DATA LINE(4) / ’ 84 84 ’ / READ ( LINE, FMT=20, REC=3 ) M, N 20 FORMAT( I4, I4 ) PRINT *, M, N END demo% f95 -f77 intern3.f demo% a.out 83 83 demo%

2.6 大端字节序和小端字节序平台之间的二进制 I/O

在 SPARC 和 x86 平台间移动时，新的编译器标志 -xfilebyteorder 支持二进制 I/O 文件。标志标识未格式化 I/O 文件的字节顺序和字节对齐。

例如，

-xfilebyteorder=little4:%all,big16:20

将指定所有文件（以 "SCRATCH" 打开的文件除外）包含 4 字节边界对齐（例如 32 位 x86）的“小端字节序”数据，Fortran 单元 20 除外，该单元是 64 位“大端字节序”文件（例如 64位 SPARC V9）。

有关详细信息，请参见 f95(1) 手册页或《Fortran 用户指南》。

2.7 传统 I/O 注意事项

Fortran 95 及传统 Fortran 77 程序在 I/O 上是兼容的。包含 f77 和 f95 混合编译代码的可执行文件可以同时从程序的 f77 和 f95 部分对同一单元执行 I/O 操作。

但是，Fortran 95 还提供了一些附加功能：

• ADVANCE=’NO’ 允许进行非提前式 I/O 操作，如下所示：

write(*,’(a)’,ADVANCE=’NO’) ’Enter size= ’ read(*,*) n

• NAMELIST 输入功能：

  • f95 允许输入时在组名前使用 $ 或 &。Fortran 95 标准只接受 &，并且这是 NAMELIST 写入语句的输出内容。

  • f95 接受 $ 作为输入组的终止符号，除非组中的最后一个数据项为 CHARACTER（在这种情况下，$ 被视为输入数据）。

  • f95 允许 NAMELIST 输入开始于记录的第一列。

  正如 f77 所做的那样，f95 承认并实现了 ENCODE 和 DECODE。

有关 f95 和 f77 间的 Fortran 95 I/O 扩展及兼容性方面的其他信息，请参见《Fortran 用户指南》。

地震数据格式读写

module weisegy
implicit none

    contains
!weiyiwen 2021/7
  subroutine weiwrite_segy(sisvx,sisvy,nt,nrec,DELTAT)
  use SGYmodule
  implicit none
  integer(kind=4) nt,nrec
 real(kind=4) DELTAT
 real(kind=4) sisvx(nt,nrec)
 real(kind=4) sisvy(nt,nrec)
  integer(kind=4) irec,it
  character(len=128) filename
    integer(kind=1)::SgyFormat=1
  real(kind=4)::tmpwei(nrec,nt)
  integer(kind=2),dimension(1800):: REEL !3600 BYTES
  integer(kind=2)::head(nrec,120)

! X component
reel=0
head=0
reel(1601:1606)=0
reel(1607)=nrec
reel(1609)=int(deltat*10**6)!weiyiwen
reel(1610)=int(deltat*10**6)
reel(1611)=nt!weiyiwen
head(:,58)=nt
head(:,59)=int(deltat*10**6)
HEAD(:,11)=1
HEAD(:,12)=1
head(:,5)=1
  sgyformat=1
  reel(1613)=5!IEEE 4BYTES!reel(1613)=1 !IBM WEIYIWEN 2021/7/7
!  write(filename,"('Vx_file_',i3.3,'.sgy')") nrec
  filename='vx.sgy'

   do irec=1,nrec
    HEAD(irec,1)=IREC
    head(irec,3)=irec
    head(irec,7)=irec
    head(irec,13)=irec
     do it=1,nt
    tmpwei(irec,it)=sngl(sisvx(it,irec)*10**3)
    enddo
   enddo

call    WriteSgy( FileName,SgyFormat,nrec,NT,tmpwei,REEL,Head )
!Y component
   do irec=1,nrec
     do it=1,nt
    tmpwei(irec,it)=sngl(sisvy(it,irec)*10**3)
    enddo
   enddo
!    write(file_name,"('Vy_file_',i3.3,'.sgy')") irec
  filename='vy.sgy'
call writeSgy( FileName,SgyFormat,nrec,NT,tmpwei,REEL,Head )

  end subroutine weiwrite_segy
end module weisegy
 

!C| Module        : SGY文件读写模块(简版)
!C| Last modified :10-22-2019                                            |
!C| Written by    : LAUQRU                                          |
!C| 主要功能：   1、支持IBM格式、IEEE格式、无卷头PC格式                  |
!C|               2、自动判断文件格式                                     |
!C|               3、自动读取采样长度、采样率及道数                       |
!C+-----------------------------------------------------------------------+
!C+-----------------------------------------------------------------------+

module SGYmodule

implicit none

    contains

!C+-----------------------------------------------------------------------+

!C| Subroutine    : GetFileLength                                         |

!C| Last modified : 10-22-2019                                           |

!C| Written by    : LAUQRU                                           |

!C| 主要功能：   1、获得文件大小                                         |

!C+-----------------------------------------------------------------------+

subroutine  GetFileLength( FileName,Length )

!    USE DFPORT

    implicit none

    character(len=128)::FileName

    integer(kind=4)::Length

    integer(kind=4) ierr, statb(13)

!    ierr = STAT ( FileName, statb )
    inquire(file=FILENAME,size=length)!file by
    !  WRITE(*,*) LENGTH
 !   if ( ierr .ne. 0 ) stop 'stat: error'

!    Length = statb(8)
!   if(ierr.eq.-1) stop

end subroutine

!C+-----------------------------------------------------------------------+

!C| Subroutine    : PREREAD                                               |

!C| Last modified : 10-22-2017                                            |

!C| Written by    : LAUQRU                                           |

!C| 主要功能：   1、预读SGY文件，获取相关参数                            |

!C+-----------------------------------------------------------------------+

!C| SgyFormat：0--IBM,1--IEEE,2--无卷头PC                                 |

!C+-----------------------------------------------------------------------+

subroutine  PREREAD( FileName,SgyFormat,NTR,NT,DT )

    implicit none

    character(len=128)::FileName

    integer:: ios,mark_pos=1

    character(len=128):: iom='OK'

    integer(kind=4)::Length,NTR,NT,DT,tmp

    integer(kind=1)::SgyFormat

    integer(kind=2),dimension(1800):: REEL !3600字节

    integer(kind=2),dimension(120):: HEAD  !240字节

    call GetFileLength(FileName, Length)

 !   WRITE(*,*)Length

!==============================================================================

!  打开SEG-Y文件

!==============================================================================

    open(70,file=FileName,       &

       access='stream',          &    !<< 流访问,即按字节访问文件

       form='unformatted',       &

       CONVERT='big_ENDIAN',     &    !<< IEEE Big Endian  ^^ Little_ENDIAN

       iostat=ios,               &

       iomsg=iom)

    read(70)REEL

    dt = reel(1609)!weiyiwen

    nt = reel(1611)!weiyiwen
  !  print *,dt,nt

    tmp = mod(Length - 3600, nt * 4 +240)

    read(70)HEAD

    if((tmp .eq.0 ).and. (HEAD(58) .eq. nt) .and.( head(59) .eq.dt).and.(dt.gt.0)) then

        SgyFormat = 0

        NTR = (Length - 3600) / (nt * 4 +240)

        write(*,*)"IBM Format"

    else

        close(70)

        mark_pos=1

        open(70,file=FileName,    &

        access='stream',          &    !<< 流访问,即按字节访问文件

        form='unformatted',       &

        CONVERT='Little_ENDIAN',  &    !<< IEEE Big Endian  ^^ Little_ENDIAN

        iostat=ios,               &

        iomsg=iom)

        rewind(unit=70)

        read(70)REEL

        dt = reel(1609)

        nt = reel(1611)

        tmp = mod(Length - 3600, nt * 4 +240)

        read(70)HEAD

        if(tmp .eq.0 .and. HEAD(58) .eq. nt .and. head(59) .eq.dt) then

             SgyFormat = 1

             NTR = (Length - 3600) / (nt * 4 +240)

    !        write(*,*)"IEEE Format"

        else

            rewind(unit=70)

            read(70)HEAD

            dt = head(59)!weiyiwen

            nt = HEAD(58)!weiyiwen

            tmp = mod(Length, nt * 4 +240)

            if(tmp .eq.0 ) then

                 NTR = Length/ (nt * 4 +240)

                 SgyFormat = 2

          !      write(*,*)"IEEE 无卷头Format"

            end if

        end if

    end if

    close(70)

end subroutine

!C+-----------------------------------------------------------------------+

!C| Subroutine    : READSgy                                               |

!C| Last modified : 10-22-2019                                            |

!C| Written by    : LAUQRU                                          |

!C+-----------------------------------------------------------------------+

!C| SgyFormat：0--IBM,1--IEEE,2--无卷头PC                                 |

!C+-----------------------------------------------------------------------+

subroutine READSgy( FileName,SgyFormat,NTR,NT,date,REEL,SgyHead )

    implicit none

    character(len=128)::FileName

    integer(kind=4)::NTR,NT,i

    integer(kind=1)::SgyFormat

    integer:: ios

    character(len=128):: iom='OK'

    integer*2,dimension(1800):: REEL !3600字节

!    integer(kind=2),dimension(120):: HEAD  !240字节

    real*4::date(NTR,NT)

    integer(kind=2)::SgyHead(NTR,120)

    if(SgyFormat .EQ. 0) then !IBM格式文件的读取

        open(70,file=FileName,    &

        access='stream',          &    !<< 流访问,即按字节访问文件

        form='unformatted',       &

        CONVERT='BIG_ENDIAN',     &    !<< IEEE Big Endian  ^^ Little_ENDIAN

        iostat=ios,               &

        iomsg=iom)

        read(70)REEL

        do i=1,Ntr

            read(70)SgyHead(i,:),date(i,:)

        end do

        close(70)

    else if(SgyFormat .EQ. 1) then !IEEE文件的读取

        open(70,file=FileName,    &

        access='stream',          &    !<< 流访问,即按字节访问文件

        form='unformatted',       &

        CONVERT='Little_ENDIAN',     &    !<< IEEE Big Endian  ^^ Little_ENDIAN

        iostat=ios,               &

        iomsg=iom)

        read(70)REEL

        do i=1,Ntr

            read(70)SgyHead(i,:),date(i,:)

        end do

        close(70)

    else if(SgyFormat .EQ. 2) then  !无卷头PC格式的读取

        open(70,file=FileName,    &

        access='stream',          &    !<< 流访问,即按字节访问文件

        form='unformatted',       &

        CONVERT='Little_ENDIAN',     &    !<< IEEE Big Endian  ^^ Little_ENDIAN

        iostat=ios,               &

        iomsg=iom)

        do i=1,Ntr

            read(70)SgyHead(i,:),date(i,:)

        end do

        close(70)

    end if

end subroutine

!C+-----------------------------------------------------------------------+

!C| Subroutine    : WriteSgy                                              |

!C| Last modified : 10-22-2019                                            |

!C| Written by    : LAUQRU                                           |

!C| 主要功能：   1、写SGY文件                                            |

!C+-----------------------------------------------------------------------+

!C| SgyFormat：0--IBM,1--IEEE,2--无卷头PC                                 |

!C+-----------------------------------------------------------------------+

subroutine  WriteSgy( FileName,SgyFormat,NTR,NT,date,REEL,SgyHead )

    implicit none

    character(len=128)::FileName

    integer(kind=4)::NTR,NT,i

    integer(kind=1)::SgyFormat

    integer:: ios

    character(len=128):: iom='OK'

    integer(kind=2),dimension(1800):: REEL !3600字节

!    integer(kind=2),dimension(120):: HEAD  !240字节

    real*4::date(NTR,NT)

    integer(kind=2)::SgyHead(NTR,120)

    if(SgyFormat .EQ. 0) then !IBM

        open(70,file=FileName,    &

        access='stream',          &

        form='unformatted',       &

        CONVERT='BIG_ENDIAN',     &

        iostat=ios,               &

        iomsg=iom)

        write(70) REEL

        do i=1,Ntr

            write(70) SgyHead(i,:),date(i,:)

        end do

        close(70)

    else if(SgyFormat .EQ. 1) then

        open(70,file=FileName,    &

        access='stream',          &

        form='unformatted',       &

        CONVERT='Little_ENDIAN',     &

        iostat=ios,               &

        iomsg=iom)

        write(70)REEL

        do i=1,Ntr

            write(70)SgyHead(i,:),date(i,:)

        end do

        close(70)

    else if(SgyFormat .EQ. 2) then  !无卷头PC格式的读取

        open(70,file=FileName,    &

        access='stream',          &

        form='unformatted',       &

        CONVERT='Little_ENDIAN',     &

        iostat=ios,               &

        iomsg=iom)

        do i=1,Ntr

            write(70)SgyHead(i,:),date(i,:)

        end do

        close(70)

    end if

end subroutine

end module
 

BIGENDIAN与little_endian的字节交换

program hello
implicit none!将两个2bytes合并成一个4bytes的数据
integer(kind=2)::x=3   !0000000000000011,高位变成低位而且倒序16BIT
integer(kind=2)::y=1!0000000000000001，低位编程高位而且倒序16BIT
integer(kind=4)::z1(2)
LOGICAL::BIG_ENDIAN=.TRUE.

REAL*4::F
print *,x,y
IF(BIG_ENDIAN.EQV..TRUE.) THEN
z1(1)=IOR(LSHIFT(x,8),RSHIFT(X,8))
Z1(2)=LSHIFT(IOR(LSHIFT(y,8),RSHIFT(Y,8)),16)
!F=REAL(z1,KIND=4)
F=IOR(Z1(1),Z1(2))
else
z1(1)=Y
Z1(2)=LSHIFT(X,16)
!F=REAL(z1,KIND=4)
F=IOR(Z1(1),Z1(2))
ENDIF
PRINT *,F
end program

Binary (unformatted) files are 100% compatible among all fortran compilers and all c compilers on all hardware platforms. You just have to know how to read them! Some conventions will make this trivial as opposed to a hardship. I regularly exchange data among Sun, Alpha, Athlon, Intel, G4, IBM SP, you name it, and have figured out how to make it work seamlessly -- mostly by using netCDF.

• ALWAYS use direct access for binary output. This avoids currupting the file with headers on each record and allows you to read the data in arbitrary chunks, not just the chunks of the writing program. With sequential access, you not only need to know how big the data array is, but how it was written. Also, sharing between c and fortran causes real headaches with sequential access.

  Direct access is awkward for saving mixed data to a single file. That's probably a bad idea anyway -- how are you going to remember what's in the file? Use a self-describing format like ascii or netCDF for mixed-data files. Or use separate files for each data type with filenames that tell you what's in them.

  Note that on some compilers (eg Dec Fortran) you need to specify the record length is measured in bytes (eg -assume byterecl). Thus (in f77):

        program write
        integer m,n, sizeofreal
        parameter (m=10, n=10, sizeofreal=4)
        real x(m,n)
  
  c     insert some code to fill x
  
  c     open direct access unformatted file with records sized for the whole array
  
        open (unit=1,
       $     file='out.r4',
       $     form='unformatted',
       $     access='direct',
       $     recl=m*n*sizeofreal)
  
  c     write the whole thing as one record
  
        write (1,rec=1) x
  
        end
  
        program read
  
        integer m,n
        parameter (m=10, n=10)
        real x(m,n)
  
  c     open direct access unformatted file with records sized for a single row
  c     note how we can read in the data differently from how it was written.
  c     You just need to make recl the smallest thing you will read and adjust
  c     your read statements accordingly
  
        open (unit=1,
       $     file='out.r4',
       $     form='unformatted',
       $     access='direct',
       $     recl=m*4)
  
  c     read in each row as a single record
  
        do j = 1, n
           read (1,rec=j) (x(i,j), i=1,m)
        end do
  
        end
  
  

• The above will work when run on any two machines that have the same byte order convention. Intel/AMD and Alpha (Dec/Compaq/HP) are so-called "little-endian". Mac, Sun, HP (risc), SGI, and most everything else are "big endian". Big endian used to be standard, but now with the rise of intel beowulf clusters, you get a lot of litte-endian. Two solutions:

  1. Use a compiler flag to force big endian (or little endian if you prefer) so all your platforms work the same. This only works if you control both writing and reading. Eg, in Dec fortran, this is -convert big_endian

  2. Detect the byte order of your machine and swap as needed. Again, you need to know what convention your input file is to know if it is different from the machine you are reading from (usually easy to figure out since the data will be garbage if you have it wrong). The following routines will effect a byte swap for 2-byte integers and 4-byte reals:

           subroutine endian(litend)
     
     c     checks if this is a little endian machine
     c     returns litend=.true. if it is, litend=.false. if not
     
           integer*1 j(2)
           integer*2 i
           equivalence (i,j)
           logical litend
     
           i = 1
           if (j(1).eq.1) then
              litend = .true.
           else
              litend = .false.
           end if
     
           end
     
           subroutine byteswapi2(k)
     
     c     does a byteswap on integer2 number
     
           integer*1 ii(2), jj(2)
           integer*2 i, j, k
           equivalence (i,ii)
           equivalence (j,jj)
     
           i = k
     
           jj(1) = ii(2)
           jj(2) = ii(1)
     
           k = j
     
           end
     
     
           subroutine byteswapr4(r)
     
     c     does a byteswap on real*4 number
     
           integer*1 ii(4), jj(4)
           real*4 r, s, t
           equivalence (s,ii)
           equivalence (t,jj)
     
           s = r
     
           jj(1) = ii(4)
           jj(2) = ii(3)
           jj(3) = ii(2)
           jj(4) = ii(1)
     
           r = t
     
           end
     
     

     When using NAGware f95, you need the -dusty flag to compile the above, so there is evidently something non-standard or obsolete about it. There is probably a more elegant solution, but this works.

• If someone does send you a sequential-access file, corrupting the data with garbage headers, you'll need the code they used to write the file and hope your machine uses the same record headers. You can still do the byte swapping. If record headers are different, open as direct access with recl=4 and experiment to find out how many bytes (records on read) you need to skip over to get to the data. Tell them never to do that again.

• The best solution, by far, is to use a platform-independent file format like netCDF or HDF -- or even gzipped ascii (use system calls to gunzip on the fly).

/***************************************************************************** 
*                                                                               
*                         NCSA HDF version 3.10                                 
*                               July 1, 1990                                    
*                                                                               
* NCSA HDF Version 3.10 source code and documentation are in the public         
* domain.  Specifically, we give to the public domain all rights for future                                                                          
* licensing of the source code, all resale rights, and all publishing rights.                                                                         
*                                                                               
* We ask, but do not require, that the following message be included in all                                                                             
* derived works:                                                                
*                                                                               
* Portions developed at the National Center for Supercomputing Applications at                                                                 
* the University of Illinois at Urbana-Champaign.                               
*                                                                               
* THE UNIVERSITY OF ILLINOIS GIVES NO WARRANTY, EXPRESSED OR IMPLIED, FOR THE                                                                         
* SOFTWARE AND/OR DOCUMENTATION PROVIDED, INCLUDING, WITHOUT LIMITATION,                                                                     
* WARRANTY OF MERCHANTABILITY AND WARRANTY OF FITNESS FOR A PARTICULAR PURPOSE                                                                         
*                                                                               
*****************************************************************************/  
/* ------------------------------------------------------------------- *        
 * procedure: convert_float_ibm_to_ieee32(float ibm,float ieee,int *n) *
 * ------------------------------------------------------------------- *        
 * PROCESSOR    - High C Compiler R2.1n                                *        
 *                                                                     *        
 * DEPENDENCES  - NONE                                                 *        
 *                                                                     *        
 * ATTRIBUTES   - REENTERANT                                           *        
 *                                                                     *        
 * ENTRY POINT  - _cfsi32                                              *        
 *                                                                     *        
 * STATUS       - NEW:           22 February 1990                      *        
 *                LAST REVISION: 19 July     1990                      *        
 *                                                                     *        
 *                            Val I. Garger, Technology Integration    *        
 *                            Group, CNSF, Cornell University          *        
 *                                                                     *        
 *                               vig@eagle.cnsf.cornell.edu            *        
 *-------------------------------------------------------------------- *        
 *                                                                     *        
 *  COPYRIGHT -  VAL GARGER, CORNELL NATIONAL SUPERCOMPUTER FACILITY,  *        
 *               (JUNE 1990) CORNELL UNIVERSITY, ITHACA, NY.           *        
 *               CONTAINS RESTRICTED MATERIALS OF CORNELL UNIVERSITY,  *        
 *               (C) COPYRIGHT CORNELL UNIVERSITY 1990                 *        
 *                                                                     *        
 * ------------------------------------------------------------------- *        
 * Change log:                                                         *        
 * ------------------------------------------------------------------- *        
 * Convert floating point, 32-bit IBM to 32-bit IEEE standard.         *        
 *                                                                     *        
 *   Example: call cfsi32(ibm, ieee, n)                                *        
 *                                                                     *        
 *  input: ibm    Array of IBM floating point numbers, REAL*4 values.  *        
 *         n      Number of elements in ibm to convert, integer.       *        
 * output: ieee   Array of 32-bit IEEE floating point numbers,         *        
 *                  single  precision.                                 *        
 *                                                                     *        
 * Format (bits, left to right):            |    Exponent bias:        *        
 *              sign   exponent   mantissa  |                          *        
 *  IBM           1      7           24     |     64 hex               *        
 *  IEEE          1      8           23     |     127                  *        
 *                                          |                          *        
 * Usage notes:                                                        *        
 * 1. Data could be converted "inplace".                               *        
 * 2. IBM values that do not conform to IEEE standard are converted to *        
 *    either infinite IEEE values (positive or negative)  or to zero.  *        
 * 3. Non-normalized with zero exponent values are kept intact.        *        
 * 4. Conversion does not incur the loss of mantissa accuracy.         *        
 * =================================================================== *        
 */                                                                             
#include <stdint.h>
#define   exp   0x7F000000                                                      
#define   sign  0x80000000                                                      
#define   tiss  0x00FFFFFF                                                      
#define   etis  0x007FFFFF                                                      
#define   nrm   0x00F00000                                                      
                                                                                
#pragma linkage (cfsi32, fortran)                                               

int convert_float_ibm_to_ieee32(int ibm[], int ieee[], int* n)
{
  int status = 0;
  int32_t j, ibs, ibe, ibt, it, k;
  union { int32_t i; float r; } u;

  for(j = 0; j < *n; j++) {
    ibs = ibm[j] & sign;
    ibe = ibm[j] & exp ;
    ibt = ibm[j] & tiss;
    if (ibt == 0) {
      ibe = 0 ;
    } else {
      if ( (ibe != 0) && (ibt & nrm) == 0 ) {
        u.i = ibm[j] ;
        u.r = u.r + 0e0 ;
        ibe = u.i & exp ;
        ibt = u.i & tiss ;
      }
      /* mantissa */
      it = ibt << 8;
      for (k = 0; (k < 5) && (it >= 0); k++ ) {
        it = it << 1;
      }
      if ( k < 4 ) {
        ibt = (it >> 8) & etis;
        ibe = (ibe >> 22) - 256 + 127 - k - 1;
        if (ibe < 0) {
          ibe = ibt = 0;
        }
        if (ibe >= 255) {
         ibe = 255; ibt = 0;
        }
        ibe = ibe << 23;
       }
    }
    ieee[j] = ibs | ibe | ibt;
    if (ibe == 255<<23) {
      status = -1;
    }
  }
  return status;
}

C++ CODE

/*******************************************************************************
Conversion between IBM and IEEE floating point 

Reference:
https://en.wikipedia.org/wiki/IBM_hexadecimal_floating-point
https://www.crewes.org/Documents/ResearchReports/2017/CRR201725.pdf
*******************************************************************************/
#include <cmath>
#include <iostream>

float ibm2ieee_intuitive(unsigned int ibm)
{
    char bias = 64;
    char sign = (ibm & 0x80000000) >> 31;
    char expn = (ibm & 0x7F000000) >> 24;
    expn = expn - bias;

    int frac = (ibm << 8);

    int shift = 0;
    do
    {
        frac = frac << 4;
        shift++;
    } while (frac != 0);

    int fracNoPadding = (ibm << 8);
    fracNoPadding = fracNoPadding >> (sizeof(unsigned int) * 8 - shift * 4);

    float rst = pow(-1, sign) * fracNoPadding * pow(16, expn - shift);
    return rst;
}

float ibm2ieee(unsigned int ibm)
{
    // copy sign
    unsigned int signBit = (ibm & 0x80000000);

    // fraction
    int frac = (ibm & 0x00FFFFFF);
    int shift = 0;
    do
    {
        // According to IBM standard, the fraction part is normalized by shifting 
        // 4bits at a time, until the leftmost digit is 0
        // Meaning the first bit bust be 0, so a do-while is correct
        frac = frac << 1;
        shift++;
    } while ((frac & 0x00800000) != 0x00800000);

    // Handle implicit 1
    frac = frac & 0x007FFFFF;
    shift++;

    // exponent
    char bias16 = 64;
    char expn16 = (ibm & 0x7f000000) >> 24;
    expn16 = expn16 - bias16;

    char expn2 = 4 * expn16 - shift;
    char bias2 = 127;
    expn2 = expn2 + bias2;
    unsigned int expn = expn2;
    expn = (expn & 0x000000FF) << 23;

    unsigned int buffer = signBit | expn | frac;
    return reinterpret_cast<float&>(buffer);
}

unsigned int ieee2ibm(float ieee)
{
    unsigned int buffer = reinterpret_cast<unsigned int&>(ieee);

    // Get the sign bit
    unsigned int signBit = (buffer) & 0x80000000;
    
    // exponent
    int base2Expn = ((buffer) & 0x7F800000) >> 23;
    int ieeeBias = 127;
    base2Expn = base2Expn - ieeeBias; // remove ieee bias

    // fraction
    unsigned int frac = (buffer) & 0x007FFFFF;
    // add explicit 1
    frac = (frac | 0x00800000);

    // normalize
    // Because the first bit of the fraction in IEEE must be 1, 
    // right shift 4 bit will finish the normalization
    frac = frac >> 3;
    base2Expn = base2Expn + 4;

    // convert from base 2 to base 16
    unsigned int remainder = base2Expn % 4;
    unsigned base16Expn = (base2Expn - remainder) / 4;
    char ibmBias = 64; 
    base16Expn = base16Expn + ibmBias; // add ibm bias

    // multiply the remainder to the fraction 
    frac = frac << (remainder);

    unsigned int ibm = signBit | (base16Expn << 24) | frac;
    return ibm;
}

int main()
{
    unsigned int ibmFloat[6];
    float ieeeFloat[6];
    // Case 1: Wiki example
    // IEEE -118.625, 0xC2ED4000, 11000010 11101101 01000000 00000000
    // IBM            0xC276A000, 11000010 01110110 10100000 00000000
    ibmFloat[0] = 0xC276A000; 
    ieeeFloat[0] = -118.625;

    // Case 2: 
    // IEEE 0.15625
    // IBM 0x40280000
    ibmFloat[1] = 0x40280000;
    ieeeFloat[1] = 0.15625;
    
    // Test ibm to ieee
    for (int i = 4; i < sizeof(ibmFloat) / sizeof(ibmFloat[0]); i++)
    {
        unsigned int t = ibmFloat[i];
        float f1 = ibm2ieee_intuitive(t);
        float f2 = ibm2ieee(t);

        if (f1 != ieeeFloat[i])
            std::cout << "Test case " << i << " failed for ibm2ieee_intuitive()" << std::endl;
        if (f2 != ieeeFloat[i])
            std::cout << "Test case " << i << " failed for ibm2ieee()" << std::endl;
    }

    // Test ieee to ibm
    for (int i = 0; i < sizeof(ieeeFloat) / sizeof(ieeeFloat[0]); i++)
    {
        unsigned int ui1 = ieee2ibm(ieeeFloat[i]);
        if(ui1 != ibmFloat[i])
            std::cout << "Test case " << i << " failed for ieee2ibm()" << std::endl;
    }

    return 0;
}