概述
让MATLAB调用fortran代码写的程序,fortran的程序往往具有多个输入参数和多个输出参数,如何编写MEX文件的源代码?
利用MATLAB自带的一个示例文件进行了修改测试。
示例文件是 xtimesy.F。
这个程序是计算一个标量与一个矩阵相乘的结果。输入为一个标量和一个m×n的矩阵,输出为一个m×n的矩阵。这个示例可以作为多个输入,单个输出的情况的一个例子。需要提示的是,在对 xtimesy.F mex后,使用该函数时,应该这么用
首先给x和y赋值,然后
z=xtimesy(x,y)
如果只是输入 xtimesy(x,y),会提示
错误使用 xtimesy
One output required.
xtimesy.F源代码
#include "fintrf.h"
C======================================================================
#if 0
C
C xtimesy.F
C .F file needs to be preprocessed to generate .for equivalent
C
#endif
C
C xtimesy.f
C
C Multiple the first input by the second input
C This is a MEX file for MATLAB.
C Copyright 1984-2018 The MathWorks, Inc.
C
C======================================================================
C Gateway routine
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
C Declarations
implicit none
C mexFunction arguments:
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
C Function declarations:
mwPointer mxCreateDoubleMatrix
#if MX_HAS_INTERLEAVED_COMPLEX
mwPointer mxGetDoubles
#else
mwPointer mxGetPr
#endif
mwPointer mxGetM, mxGetN
integer mxIsNumeric
C Pointers to input/output mxArrays:
mwPointer x_pr, y_pr, z_pr
C Array information:
mwPointer m, n
mwSize size, row
C Arguments for computational routine:
real*8 x, y(100), z(100)
C-----------------------------------------------------------------------
C Check for proper number of arguments.
if (nrhs .ne. 2) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:nInput',
+ 'Two inputs required.')
elseif (nlhs .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:nOutput',
+ 'One output required.')
endif
C-----------------------------------------------------------------------
C Validate inputs
C Check to see both inputs are numeric.
if (mxIsNumeric(prhs(1)) .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:NonNumeric1',
+ 'Input # 1 is not a numeric.')
elseif (mxIsNumeric(prhs(2)) .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:NonNumeric2',
+ 'Input #2 is not a numeric array.')
endif
C-----------------------------------------------------------------------
C Check that input #1 is a scalar.
m = mxGetM(prhs(1))
n = mxGetN(prhs(1))
if(n .ne. 1 .or. m .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:NonScalar',
+ 'Input #1 is not a scalar.')
endif
C-----------------------------------------------------------------------
C Get the size of the input matrix.
m = mxGetM(prhs(2))
n = mxGetN(prhs(2))
size = m*n
if(size .gt. 100) then
call mexErrMsgIdAndTxt ('MATLAB:xtimesy:InputTooBig',
+ 'Input #2 number of elements exceeds buffer')
endif
C-----------------------------------------------------------------------
C Create matrix for the return argument.
plhs(1) = mxCreateDoubleMatrix(m,n,0)
#if MX_HAS_INTERLEAVED_COMPLEX
x_pr = mxGetDoubles(prhs(1))
y_pr = mxGetDoubles(prhs(2))
z_pr = mxGetDoubles(plhs(1))
#else
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(prhs(2))
z_pr = mxGetPr(plhs(1))
#endif
C Load the data into Fortran arrays.
row = 1
call mxCopyPtrToReal8(x_pr,x,row)
call mxCopyPtrToReal8(y_pr,y,size)
C Call the computational subroutine.
call xtimesy(x,y,z,m,n)
C Load the output into a MATLAB array.
call mxCopyReal8ToPtr(z,z_pr,size)
return
end
C-----------------------------------------------------------------------
C Computational routine
subroutine xtimesy(x,y,z,m,n)
mwSize m, n
real*8 x, y(m,n), z(m,n)
mwSize i, j
do 20 i=1,m
do 10 j=1,n
z(i,j)=x*y(i,j)
10 continue
20 continue
return
end
修改
为了测试,如何输出多个结果,我在这个源代码的基础上进行了修改。增加了一个输出 k(i,j) = 2.0 * x * y(i,j) 。为了能更清楚得看到修改的地方,对修改的地方进行了加粗。但是,。。。。
由于markdown的编辑器自动识别了代码,代码块里不能再加粗,所以留下了加粗的标签,但没有加粗。凑合看吧。
#include “fintrf.h”
C======================================================================
#if 0
C mxtimesy.F
C .F file needs to be preprocessed to generate .for equivalent
C
#endif
C======================================================================
C Gateway routine
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
C Declarations
implicit none
C mexFunction arguments:
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
C Function declarations:
mwPointer mxCreateDoubleMatrix
#if MX_HAS_INTERLEAVED_COMPLEX
mwPointer mxGetDoubles
#else
mwPointer mxGetPr
#endif
mwPointer mxGetM, mxGetN
integer mxIsNumeric
C Pointers to input/output mxArrays:
mwPointer x_pr, y_pr, z_pr , k_pr
C Array information:
mwPointer m, n
mwSize size, row
C Arguments for computational routine:
real*8 x, y(100), z(100)**, k(100)**
C-----------------------------------------------------------------------
C Check for proper number of arguments.
if (nrhs .ne. 2) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:nInput',
+ 'Two inputs required.')
elseif (nlhs .ne. **2**) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:nOutput',
+ '**Two** output required.')
endif
C-----------------------------------------------------------------------
C Validate inputs
C Check to see both inputs are numeric.
if (mxIsNumeric(prhs(1)) .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:NonNumeric1',
+ 'Input # 1 is not a numeric.')
elseif (mxIsNumeric(prhs(2)) .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:NonNumeric2',
+ 'Input #2 is not a numeric array.')
endif
C-----------------------------------------------------------------------
C Check that input #1 is a scalar.
m = mxGetM(prhs(1))
n = mxGetN(prhs(1))
if(n .ne. 1 .or. m .ne. 1) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:NonScalar',
+ 'Input #1 is not a scalar.')
endif
C-----------------------------------------------------------------------
C Get the size of the input matrix.
m = mxGetM(prhs(2))
n = mxGetN(prhs(2))
size = m*n
if(size .gt. 100) then
call mexErrMsgIdAndTxt ('MATLAB:mxtimesy:InputTooBig',
+ 'Input #2 number of elements exceeds buffer')
endif
C Create matrix for the return argument.
plhs(1) = mxCreateDoubleMatrix(m,n,0)
**plhs(2) = mxCreateDoubleMatrix(m,n,0)**
#if MX_HAS_INTERLEAVED_COMPLEX
x_pr = mxGetDoubles(prhs(1))
y_pr = mxGetDoubles(prhs(2))
**k_pr = mxGetDoubles(plhs(1))**
z_pr = mxGetDoubles(plhs(2))
#else
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(prhs(2))
**k_pr = mxGetPr(plhs(1))**
z_pr = mxGetPr(plhs(2))
#endif
C Load the data into Fortran arrays.
row = 1
call mxCopyPtrToReal8(x_pr,x,row)
call mxCopyPtrToReal8(y_pr,y,size)
C Call the computational subroutine.
call xtimesy(x,y,**k,z,**m,n)
C Load the output into a MATLAB array.
call mxCopyReal8ToPtr(z,z_pr,size)
**call mxCopyReal8ToPtr(k,k_pr,size)**
return
end
C-----------------------------------------------------------------------
C Computational routine
subroutine xtimesy(x,y,**k,z**,m,n)
mwSize m, n
real*8 x, y(m,n), z(m,n),** k(m,n)**
mwSize i, j
do 20 i=1,m
do 10 j=1,n
**k(i,j) = 2.0 * x * y(i,j)**
z(i,j) = x * y(i,j)
10 continue
20 continue
return
end