Writing Matlab MEX files using Fortran 90 / 95 (and the Intel Fortran Compiler) (Update: This text refers primarily to Matlab 6.5 (aka R13) on Linux. The new Matlab 7 introduces some serious problems with non-GNU compilers due to throwing C++ exceptions. Read on below for details.) The Mathworks do not officially support Fortran 90 / 95, but thanks to people like Benjamin Barrowes who wrote Matlab2FMex and the people at NAGWare it is easy to figure out how to write your MEX files in Fortran 95 without using the awkward mxCopy-functions. Lately, I also found this page which describes in detail all three possible ways of writing MEX files entirely in Fortran. However, of these three ways, two make a local copy of the data which is less than optimal when dealing with big data sets. The trick is to declare the pointer you get from mxGetPr as an integer pointer and pass it along with the array dimensions to a subroutine which expects a double precision array. Until now this works for me without any problems (Linux w/ Intel Compiler), but note that this procedure may not work with every compiler or on any platform. I may also add that this procedure does NOT work with the F programming language as it forbids passing "wrong" arguments to subroutines. First, we will have to explicitly define the used Matlab API functions, at least those who have a return value: module mexf90 interface function mxGetPr(pm) integer,pointer :: mxGetPr integer :: pm end function mxGetPr function mxGetM(pm) integer :: mxGetM,pm end function mxGetM function mxGetN(pm) integer :: mxGetN,pm end function mxGetN function mxCreateDoubleMatrix(m,n,type) integer :: m,n,type,mxCreateDoubleMatrix end function mxCreateDoubleMatrix function mxGetScalar(pm) integer :: pm double precision :: mxGetScalar end function mxGetScalar end interface end module mexf90 Just save this file under the name "mexf90.f90" and compile it. (For those unfamiliar with modules: you must prevent the compiler from also linking the file - this is usually done with the "-c" switch.) This module contains only the API functions we need, but adding further functions is straightforward. You may also click here to obtain a module file with more function definitions. Now the actual function, consisting of the subroutine "mexFunction" (the gateway to Matlab) and an example subroutine "scalarMult" which multiplies the matrix A with a scalar "alpha" and returns the result in matrix B. The comments should make clear how these functions work: subroutine mexFunction(nlhs, plhs, nrhs, prhs) use mexf90 ! API function definitions implicit none integer, intent(in) :: nlhs, nrhs integer, intent(in), dimension(*) :: prhs integer, intent(out), dimension(*) :: plhs integer :: m,n integer, pointer :: A,B ! These are the pointers to the matrix data double precision :: alpha ! Check input arguments if(nrhs /= 2) then call mexErrMsgTxt('Function requires two input arguments.'); end if ! Get data and size of the input matrix A = mxGetPr(prhs(1)) m = mxGetM(prhs(1)) n = mxGetN(prhs(1)) ! Get scalar value alpha = mxGetScalar(prhs(2)) ! Create output matrix plhs(1) = mxCreateDoubleMatrix(m,n,0) B = mxGetPr(plhs(1)) ! Call subroutine for multiplication call scalarMult(A,B,alpha,m,n) end subroutine mexFunction subroutine scalarMult(A,B,alpha,m,n) implicit none integer :: m,n double precision :: alpha ! Now define the matrices with the actual data type and dimension double precision, dimension(m,n) :: A,B ! Now do something useful... B = alpha * A end subroutine scalarMult For more complex functions it is a bit tedious to pass the dimensions along with the data pointers - if someone finds a better solution, please let me know! The Intel Fortran Compiler for Linux is free for non-commercial use, but you probably already know that... These are the relevant parts in the mexopts.sh file for MEXing the above file with the IFC FC='ifort' FFLAGS='-fPIC -u -w95 -warn all' FLIBS='$RPATH $MLIBS -L/opt/intel_fc_80/lib -lm' FOPTIMFLAGS='-O3' FDEBUGFLAGS='-g' LD='icc' LDFLAGS='-pthread -shared' LDOPTIMFLAGS='$FOPTIMFLAGS' LDDEBUGFLAGS='-g' As you can see, I do not use any map-files as that didn't work for me. It is important to note that linking with "ifort" will usually NOT work, the above file uses the Intel C++ Compiler "icc" instead. If you do not have icc installed, the GNU compiler "g++" should also work. If you also want to use functions like "print", e.g. for debugging purposes (mexprintf does not allow format strings in Fortran...), you will also have to link "ifcore", but you should remove that for the release version of your code. If you get error messages about unresolved symbols when executing the mex-file, make sure that the Intel libraries are scanned by ldconfig or put them into your LD_LIBRARY_PATH before starting Matlab. If you get an error message saying that __MAIN couldn't be found, check again that you are NOT linking with ifort. If you still have problems with unresolved symbols, try using the "-static-libcxa" switch. Now happy MEXing! Oh, and also check out LAPACK95, incredibly useful for porting M-files to MEX. And Matran looks very promising, although I didn't use it yet. Update: Problems with Matlab 7 (aka R14) The new Matlab 7 now uses C++ exceptions in library functions like mexErrMsgTxt. As Matlab was built on Linux with the GNU compiler, this breaks practically all other compilers which are unable to correctly propagate these exceptions - including the Intel Fortran Compiler (and also Ocaml, by the way...). (A sidemark: the Intel C Compiler (icc) will usually work as long as you always compile with "-Kc++"). This means, as soon as a Fortran MEX file compiled with the Intel Fortran Compiler calls mexErrMsgTxt, the complete Matlab session will abort and return to the shell, as there is no handler to catch the exception. I doubt that there exists an easy solution to this problem, as the implementation of exception handling in C++ differs between different compilers. Thus, if you want to write your files entirely in Fortran 95 and you want to use the IFC, don't call mexErrMsgTxt or your session will abort. If you have to use mexErrMsgTxt, you are for now forced to use a GNU Fortran compiler. Although Fortran and C don't "know" exception handling, the GNU C and Fortran compilers have a switch "-fexceptions" which includes code to correctly propagate these exceptions. For Fortran 95 you can use either g95 or GNU Fortran 95, both work for creating Matlab MEX files as long as you compile with "-fexceptions". However, both compilers are not yet stable, but they are making rapid progress recently.
最近，好几个人问到公共块出错的问题，都是出于一个原因，特在这里把老贴整理一下供大家参考。 问题： 20924057 Compiling Fortran... D:\FRAME2D.FOR D:\FRAME2D.FOR(314) : Warning: Because of COMMON, the alignment of object is inconsistent with its type COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6)--------------------------^ D:\FRAME2D.FOR(314) : Warning: Because of COMMON, the alignment of object is inconsistent with its type COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6) ---------------------------^ 解答： 这是老Fortran的规则。看你的公共语句： COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6),SF(6),T(3,3) 根据程序的隐含规则，LM(6),ND一共是7个整型数，每个整型数4个字节，共28个字节，不是8字节的整倍数。紧跟着是双精度的ASA(6,6)，每个数是8个字节，可是ASA(1,1)不能从8字节整倍数位置开始存储，这就是the alignment of object is inconsistent with its type 的意思。 对于警告错误，可以忽略，一般不致于影响执行结果。如果要改的话，有几种办法： 1）在公共语句中，把双精度数放在前边，整形数跟在后边； 2）在ASA()前插一个整型变量（哪怕是没用的），用来占用4个字节，以使得后面的双精度数可以从8字节整数倍位置开始存储。 准则：作为好的编程习惯，建议在公共块中，把实型变量排在前边，把整型数据放在后边，就不会有对位不整的错误！ 仅供参考。 补充：（05.06.01） 公共块不是用堆栈实现的，是内存中的一段连续存储的数据。 按照Fortran的规定，当读取双精度数据时，总是假定前面的数据长度是双精度数字节长度（8个字节）的整数倍。 对于本例，ASA(1,1)从第29个字节开始存放8个字节；可是读取的时候，要从第33个字节（28不是8的倍数，32是28之后最小的8的倍数）开始读入8个字节，这就是定位（alignment）错误。 所以F90之后不提倡用公共块共享数据，而可用更为灵活的module来代替公共块共享数据。 公共块是过时的语言功能
转自：编程爱好者论坛 Fortran区的一篇置顶文章 最近，好几个人问到公共块出错的问题，都是出于一个原因，特在这里把老贴整理一下供大家参考。 问题： 20924057 Compiling Fortran... D:\FRAME2D.FOR D:\FRAME2D.FOR(314) : Warning: Because of COMMON, the alignment of object is inconsistent with its type COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6) --------------------------^ D:\FRAME2D.FOR(314) : Warning: Because of COMMON, the alignment of object is inconsistent with its type COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6) ---------------------------^ 解答： 这是老Fortran的规则。看你的公共语句： COMMON /EM/LM(6),ND,ASA(6,6),RF(6),SA(6,6),SF(6),T(3,3) 根据程序的隐含规则，LM(6),ND一共是7个整型数，每个整型数4个字节，共28个字节，不是8字节的整倍数。紧跟着是双精度的ASA(6,6)，每个数是8个字节，可是ASA(1,1)不能从8字节整倍数位置开始存储，这就是the alignment of object is inconsistent with its type 的意思。 对于警告错误，可以忽略，一般不致于影响执行结果。如果要改的话，有几种办法： 1）在公共语句中，把双精度数放在前边，整形数跟在后边； 2）在ASA()前插一个整型变量（哪怕是没用的），用来占用4个字节，以使得后面的双精度数可以从8字节整数倍位置开始存储。 准则：作为好的编程习惯，建议在公共块中，把实型变量排在前边，把整型数据放在后边，就不会有对位不整的错误！ 仅供参考。 补充：（05.06.01） 公共块不是用堆栈实现的，是内存中的一段连续存储的数据。 按照Fortran的规定，当读取双精度数据时，总是假定前面的数据长度是双精度数字节长度（8个字节）的整数倍。 对于本例，ASA(1,1)从第29个字节开始存放8个字节；可是读取的时候，要从第33个字节（28不是8的倍数，32是28之后最小的8的倍数）开始读入8个字节，这就是定位（alignment）错误。 所以F90之后不提倡用公共块共享数据，而可用更为灵活的module来代替公共块共享数据。 公共块是过时的语言功能
如何在fortran中读写文件时不换行（上传供大家学习） 转自ifelseif的博客 如何在fortran中读写文件时不换行？这是个极简单又极复杂的问题，简单到只要一个字符，复杂到翻破了好几本语法书也没找见。fortran中默认一条read或者write结束之后就换一行，但是读和写还有些不太一样。 读文件时，read之后如果写了一个数组，就像这样： read(10,*)Y(1:n) 整整一行数就全都读到数组里了。但是如果用write，写到文件中却不是这个样子，会给你一个超级长的文件然后每行只有一个数。有一个选项叫ADVANCE='YES'/'NO'，可以控制输入输出语句完了之后要不要换行，默认是'YES'，很不幸，在intel的fortran中这个选项只对read起作用，write依旧不行。 在fortran的输入输出中，/表示换行，那么\表示什么意思呢，就是不换行。这是我在网上逛了老半天才看到的，为了防止忘记，写到博客里面，立此存照。 下面是一段fortran代码样例，要处理的数据20个数就会换一行，一般来讲最后一行是不满20个数的，需要用个同余判断一下 PROGRAM MAIN IMPLICIT NONE INTEGER I,J,NY,A,B REAL X,Y(60),Z(60) OPEN(UNIT=10,FILE='SX-RIVER.TXT') OPEN(UNIT=11,FILE='SX-OUTPUT.TXT') DO I=1,9 READ(10,*) END DO DO I=1,372 !!!!!!!!! DATA INPUT !!!!!!!!!!!!!!!! READ(10,(18X,I2,1X,F7.3))NY,X CALL MOD(20,NY,A,B) DO J=1,A READ(10,(20(1X,F7.2)))Y(20*(J-1)+1:20*J) END DO DO J=1,B READ(10,(1X,F7.2),ADVANCE='NO')Y(20*A+J) END DO READ(10,*) DO J=1,A READ(10,(20(1X,F7.2)))Z(20*(J-1)+1:20*J) END DO DO J=1,B READ(10,(1X,F7.2),ADVANCE='NO')Z(20*A+J) END DO READ(10,*) !!!!!!!!! DATA OUTPUT !!!!!!!!!!!!!!! !WRITE(11,(I4,1X,F7.3))NY,X WRITE(11,(1X,F7.2,\))Y(1:NY) WRITE(11,*) WRITE(11,(1X,F7.2,\))Z(1:NY) WRITE(11,*) END DO CLOSE(10) CLOSE(11) STOP END PROGRAM MAIN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE FOR CONGRUENCE(TONGYU) ! ! Y=A*X+B ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE MOD(X,Y,A,B) IMPLICIT NONE INTEGER X,Y,A,B,TMP A=1 TMP=X DO WHILE(Y.GT.TMP) A=A+1 TMP=TMP+X END DO A=A-1 B=X-(TMP-Y) RETURN END SUBROUTINE MOD
有一个问题困扰我许多年了：同样是Fortran语言用Unformatted格式生成的数据文件，在32位和64位机器上大小不同。关于这一问题，在 这里 描述道： Update (Jan 2008): It would appear that on 64 bit machines the 2 header elements are each written as 4 bytes instead of 2 bytes each. 但也没有给出解决方案。 终于我想到一个另类做法，可以通用地读取这样的文件。 假设原来的文件是这样生成的： real*8 u_tmp(k1:kmax,k1:kmax,0:n/16-1) open(10,file=filename,form='unformatted') write(10)(((u_tmp(i,j,k),i=k1,k2),j=k1,k2),k=0,n/16-1) close(10) 如果原来文件是在64位机上生成的，那么不需要做特殊处理。否则，假设生成的文件是u.dat，则另外新建一个4字节的文件，比如叫tmp。执行 cat z u.dat u_new.dat 将4个字节插到u.dat的文件头，使得新文件的偏移量和原文件相同。 接着，用下面的程序读取 open(10,file=dat(ifn),access='direct', recl=8) do k = 0, n/16-1 do j = k1, k2 do i = k1, k2 read(10, rec = k*n*n+(j-k1)*n+i-k1+2) u_tmp(i,j,k) enddo enddo enddo 这仅仅是针对real*8的读法，但其他格式的做法也类似
我今天弄到了一本书，名字是：Sample page from NUMERICAL RECIPES IN FORTRAN 77 THE ART OF SCIENTIFIC COMPUTING 1100多页啊，刚才看了一下。这本书太牛了，几乎面面俱到。谁要是把这本书一个一个地练习一遍，我觉得吧：不但是fortran编程对他小菜一碟，就是有限元以后也不是什么难的东西了。有点夸张啊，太高兴了，原谅原谅！ 呵呵，找了一个晚上，最后是有点巧取豪夺弄过来的。谁要的话，可以说一声。大约7M，最好是留QQ邮箱吧！ 下面给出本书的目录吧： 目录 分卷压缩： part1 分卷压缩： part2 我压缩的时候没注意，这里下载之后注意要把扩展名.zip改为.rar
标签: Fortran
