USING FORTRAN90 CODE IN MATLAB MEX-FILES

 

As the standard for Fortran programming slowly advances in the direction of Fortran90, I have had the problem of creating mex-files using the fortran90 coding standard. Herefore I have developed the following method of creating and compiling the mex-files. As a help I will give an example. The program will compute the inverse of a 2x2 matrix, together with the determinant. Furthermore an optional verbose flag can be defined. The code itself is not very useful, but it contains all the possible difficulties one can expect to have during developing of mex-files.
 

STEP 1.

First write the matlab source file: a .m file  I will not elaborate on that, beautiful coding examples are available on the Internet, on your local machine and in the matlab manual

EXAMPLE

function [InvMat,Determ]=inverse(Mat,verbose);

%INVERSE compute inverse of 2x2 matrix
%
%SYNTAX: [InvMat,Determ]=inverse(Mat,verbose);
%
%INPUT:  Mat     =   2x2 matrix
%        verbose =   verbose=1: show information, verbose=0: no information,  default=0
%
%OUTPUT: InvMat  =   inverse of matrix, if no inverse exists, return conjugate of matrix
%        Determ  =   determinant of matrix
%
%AUTHOR: Jeroen Goudswaard,  J.C.M.Goudswaard@CTG.TUDelft.NL

%how many variables
if (nargin == 1);
    verbose = 0;
end;
if (verbose ~= 0);
    verbose = 1;
end;

%compute determinant

Determ = Mat(1,1)*Mat(2,2) - Mat(2,1)*Mat(1,2);

%fill Invmat
if (abs(Determ) > 2*eps);
      InvMat(1,1) =   Mat(2,2);
      InvMat(2,2) =   Mat(1,1);
      InvMat(1,2) = - Mat(1,2);
      InvMat(2,1) = - Mat(2,1);
      InvMat      =   InvMat ./ Determ;
else;
%return conjugate
      if(verbose == 1);
         disp(['No Inverse possible: Determ almost ' int2str(Determ)' ==> transpose returned']);
         % sure, I could return 0 instead of int2str(Determ),
         % but I will use this also in the f90 source, and
         % it is not trivial overthere.
      end;
      InvMat = Mat';
end;

STEP 2

Now rewrite the matlab code into f90 subroutine coding. And now we are at the beautiful part of it: it won't cost you too much effort,
as their coding standards are very much alike.
However, how impatient you might be, please read the RED COMMENTS!
It is also advisable only to ALLOCATE arrays only in the gateway subroutine, USE ASSUMED SIZE ARRAYS PREFERABLY!.
Many of my own INEXPLICABLE BUGS have been solved by removing all allocatable arrays from the actual subroutines.

EXAMPLE
 

SUBROUTINE inverse(Mat,verbose,InvMat,Determ)

  !INVERSE compute inverse of 2x2 matrix
  !
  !AUTHOR: Jeroen Goudswaard,  J.C.M.Goudswaard@CTG.TUDelft.NL

  IMPLICIT NONE

  INTEGER :: i,j,verbose
  DOUBLE PRECISION :: Determ,Mat(2,2),InvMat(2,2)
  CHARACTER*80 :: string
 

  ! compute determinant
  Determ = Mat(1,1)*Mat(2,2) - Mat(2,1)*Mat(1,2)

  !fill Invmat
  IF (ABS(Determ) > 2*eps) THEN
        InvMat(1,1) =   Mat(2,2)
        InvMat(2,2) =   Mat(1,1)
        InvMat(1,2) = - Mat(1,2)
        InvMat(2,1) = - Mat(2,1)
        InvMat      =   InvMat / Determ
  ELSE
  !return conjugate-transpose
      IF(verbose == 1)THEN
                        WRITE(string,'(f8.3)')Determ
         CALL mexprintf('No Inverse possible: Determ = '// string //'&
            & conjugate-transpose returned'//CHAR(10))   !CHAR(10) is a <CR>
! first pitfall: you CANNOT use write or print statements, because matlab has control
!     over your display so you HAVE TO USE mexprintf, which is in the f77-library of matlab
      END IF
      InvMat = TRANSPOSE(Mat);
  ENDIF
END SUBROUTINE inverse
 

STEP 3

Now the most important task comes: writing the gateway-subroutine `mexfunction', much of it is standard, but I will document this one very carefully

EXAMPLE

SUBROUTINE mexfunction(nlhs,plhs,nrhs,prhs)
   IMPLICIT NONE

   INTEGER*8 :: plhs(*), prhs(*)
   !pointers to input data: prhs and output: plhs, always take INTEGER*8, to let it work on 64-bit
   !machines (SGI e.g.) 32-bit compilers will correct this to INTEGER*4, so don't worry about the
   !warning(s) on this during compilation. (a #ifdef can also be used in a include file)
   INTEGER :: nlhs, nrhs
   !number of inputs on right and left hand side
   INTEGER*8 :: mxcreatefull, mxgetpr
   !pointers to intrinsic functions of the matlab library
   INTEGER :: mxgetm, mxgetn
   !declare the name of these intrinsics
   !up till here all statements are obliged
   INTEGER*8 :: pr_Mat, pr_verbose, pr_InvMat, pr_Determ
   !define pointer to specific data in input prhs, one pointer per variable/array
   INTEGER :: mMat, nMat
   !integers that will contain the 'size' of Mat
   INTEGER :: nverbose
   !integer to copy the value of verbose to
   DOUBLE PRECISION :: verbose, Determ
   !everything that comes out of matlab is double precision!!
   DOUBLE PRECISION, ALLOCATABLE :: Mat(:,:), InvMat(:,:)
   !let's suppose we do not know the size of Mat in advance
   !so we will use ALLOCATABLE arrays

   IF (nrhs>2) THEN
      CALL mexerrmsgtxt('more than two arguments')
   END IF
   IF (nrhs<1) THEN
      CALL mexerrmsgtxt('no input matrix specified')
   END IF
   !some checks on input, we do not want matlab to give the error, but
   !we check it ourselves
   mMat  = mxgetm(prhs(1))
   !get the first dimension of the first input, i.e. Mat
   nMat  = mxgetn(prhs(1))
   !get the second dimension of the first input, i.e. Mat

   ALLOCATE(Mat(mMat,nMat))
   !allocate memmory for the input matrix

   pr_Mat  = mxgetpr(prhs(1))
   !let the pointer point to the data
   CALL mxcopyptrtoreal8(pr_Mat , Mat , mMat*nMat)
   !copy the values from matlab to the local variable, size is mMat*nMat

   IF (nrhs==2)THEN    !check if there is a second input: verbose
       pr_verbose = mxgetpr(prhs(2))
       CALL mxcopyptrtoreal8(pr_verbose , verbose , 1)
       !verbose is of size 1, naturally
   ELSE
       verbose = 0D0
   END IF

   nverbose = NINT(verbose)
   !nverbose is the integer representation of double precision verbose

   ALLOCATE(InvMat(mMat,nMat))
   !ALLOCATE the outgoing array InvMat

   plhs(1) = mxcreatefull(mMat,nMat, 0)
   !create matrix for output
   plhs(2) = mxcreatefull(   1,   1, 0)
   !create variable for output

   pr_InvMat = mxgetpr(plhs(1))
   pr_Determ = mxgetpr(plhs(2))
   !create the output pointers

   CALL inverse(Mat,nverbose,InvMat,Determ)
   !call the Fortran90 subroutine

   CALL mxcopyreal8toptr(InvMat, pr_InvMat, mMat*nMat)
   CALL mxcopyreal8toptr(Determ, pr_Determ, 1)
   !copy the data back to matlab

   DEALLOCATE(Mat,InvMat)
   !NEVER EVER forget to DEALLOCATE dynamically allocated arrays, without this, either
   !your mex-file or matlab is very likely to CRASH!

END SUBROUTINE mexfunction
 

STEP 4

Now you are ready to compile your mex-file either from matlab or from the commandline.
First put all subroutines together in one file, and let it end on .f
matlab does not compile sourcecode ending on .f90,

In this case copy the routines from STEP2 and STEP3 together in a file inverse.f
then go to your mexopts.sh file,
for your own platform you have to specify some information, like your preferred compiler, linker and some options
in ANY case you need the option that says that it should compile .f -files as f90 source code!
Other necessary options are HIGHLY DEPENDENT on your particular platform and compiler.

。。。。。so if any of you has working mexopts.sh configurations, I would be glad to publish them on this page

And then the last part of it:

mex -O inverse.f
and now you can use inverse.mex just as inverse.m

COMMENTS, SUGGESTIONS, JUST SAYING HELLO: Jeroen@Goudswaard.org <script language="JavaScript" src="http://m1.nedstatbasic.net/basic.js"> </script> <script language="JavaScript"> </script> /

8. Can I use CXML, IMSL or the Intel MKL libraries in FMEX files?

 

Yes, I do it quite frequently. The following is written using Matlab 5.3, Compaq Visual Fortran 6.6, and Windows 2000. The example here uses Fortran95 source, though this can be done in FORTRAN77 just as easily. There are examples of MEX files that use MKL: vvmult_mkl.f90, IMSL: vvmult_imsl.f90, and CXML: vvmult_cxml.f90.

 

First, you will need to modify your mexopts.bat file. On R13 (v6.5) modify the file c:/matlab6p5/bin/win32/mexopts/df66opts.bat. On R11 and R12, you will need to issue the command prefdir(1) to tell you where your mexopts.bat is located. You can see my mexopts.bat file for R11/R12 and R13 at these links.

 

In order to add IMSL, you will need to add the following lines to your mexopts file:

 

set LINKFLAGS=%LINKFLAGS% sstatd.lib sstats.lib

set LINKFLAGS=%LINKFLAGS% smathd.lib smaths.lib sf90mp.lib

set LIB=%LIB%;%DFDir%/IMSL/LIB;

 

For MKL, you will need to add these lines, you must also uncomment one of the last three lines to specify the processor type. Note that MEX files linked against mkl_p4.lib will only run on a P4. MEX files linked against mkl_p3.lib will only run on a P3 or P4.

 

set LINKFLAGS=%LINKFLAGS% mkl_s.lib mkl_lapack.lib

set LIB=%LIB%;C:/Program Files/Intel/MKL/ia32/lib

rem ********************************************************************

rem Choose One (p4=Pentium 4, p3=Pentium 3, def=all others)

rem ********************************************************************

rem set LINKFLAGS=%LINKFLAGS% mkl_def.lib

rem set LINKFLAGS=%LINKFLAGS% mkl_p3.lib

rem set LINKFLAGS=%LINKFLAGS% mkl_p4.lib

 

For CXML, add the following lines:

 

set LINKFLAGS=%LINKFLAGS% cxml.lib

set LIB=%LIB%;%DFDir%/CXML/LIB

set INCLUDE=%INCLUDE%;%DFDir%/CXML/INCLUDE

 

Note that many of the routines in these libraries conflict. Many conflict across all 3 libraries, in fact. Therefore, if you think that you might be using different libs on different projects, you can also leave all of the LINKFLAG variables commented, and specify the libraries on the command line. So if you want to use MKL, you build the file with the command

 

? mex vvmult_mkl.f90 mkl_s.lib mkl_def.lib mkl_lapack.lib

 

Likewise, the syntax for IMSL is:

 

? mex vvmult_imsl.f90 sstatd.lib sstats.lib smathd.lib smaths.lib

 

And finally, CXML:

 

? mex vvmult_cxml.f90 cxml.lib

 

And there you go!

 

10. Can I optimize my MEX files?

 

Sure, I do it, so can you. The optimization happens in the mexopts.bat file. If you have R11/R12, issue the command prefdir(1), which will reveal its location. For R13, I hack on the mexopts file located at c:/matlab6p5/bin/win32/mexopts/df66opts.bat. If you do this, be sure to save a backup of the original! And if you use R13, you will need to run mex -setup after each change.

 

You can see my R11/R12 mexopts file here, and here is my R13 file.

/

r13_df66opts.bat.txt

@echo off

rem DF66OPTS.BAT

rem

rem    Compile and link options used for building MEX-files

rem    using the Compaq Visual Fortran compiler version 6.6

rem

rem    $Revision: 1.1 $  $Date: 2002/06/03 12:17:13 $

rem

rem ********************************************************************

rem General parameters

rem ********************************************************************

set MATLAB=%MATLAB%

set DF_ROOT=%DF_ROOT%

set VCDir=%DF_ROOT%/VC98

set MSDevDir=%DF_ROOT%/Common/msdev98

set DFDir=%DF_ROOT%/DF98

set PATH=%MSDevDir%/bin;%DFDir%/BIN;%VCDir%/BIN;%PATH%

set INCLUDE=%DFDir%/INCLUDE;%INCLUDE%

set LIB=%DFDir%/LIB;%VCDir%/LIB;%LIB%

 

rem ********************************************************************

rem Compiler parameters

rem ********************************************************************

set COMPILER=df

set COMPFLAGS=/fpp:"/m /S%MATLAB%/extern/include" /compile_only  /tune:host /nologo /DMATLAB_MEX_FILE /keep

set OPTIMFLAGS = /threads /optimize:4 /fast /arch:generic /assume:nodummy_aliases /inline:speed

rem set OPTIMFLAGS=/libs:dll /threads /optimize:4 /fast /arch:generic /assume:nodummy_aliases /inline:speed /DNDEBUG

set DEBUGFLAGS=/libs:dll /threads /dbglibs /debug:full /pdbfile

set NAME_OBJECT=/Fo

 

rem ********************************************************************

rem Linker parameters

rem ********************************************************************

set LIBLOC=%MATLAB%/extern/lib/win32/digital/df60

set LINKER=link

set LINKFLAGS=/DLL /EXPORT:_MEXFUNCTION@16 /LIBPATH:"%LIBLOC%" libmx.lib libmex.lib libmat.lib /implib:%LIB_NAME%.lib /NOLOGO

set LINKOPTIMFLAGS=

set LINKDEBUGFLAGS=/debug

set LINK_FILE=

set LINK_LIB=

set NAME_OUTPUT="/out:%OUTDIR%%MEX_NAME%.dll"

set RSP_FILE_INDICATOR=@

 

rem ********************************************************************

rem Resource compiler parameters

rem ********************************************************************

set RC_COMPILER=rc /fo "%OUTDIR%mexversion.res"

set RC_LINKER=

 

rem ********************************************************************

rem There is little conflict between IMSL & CXML or MKL, but

rem CXML & MKL conflict. I typically will leave all of the LIB variables

rem uncommented, and I uncomment the IMSL LINK_FILE, but leave the

rem other LINK_FILE variables commented. Then if I need MKL or CXML I

rem specify the libs on the mex command line.

rem ********************************************************************

 

rem ********************************************************************

rem Uncomment these lines for MKL 5.2 libs

rem One of the architecture-specific lines below must also be uncommented

rem ********************************************************************

set LINKFLAGS=%LINKFLAGS% mkl_s.lib mkl_lapack.lib

set LIB=%LIB%;C:/Program Files/Intel/MKL/ia32/lib

rem ********************************************************************

rem Choose One (p4=Pentium 4, p3=Pentium 3, def=all others)

rem ********************************************************************

rem set LINKFLAGS=%LINKFLAGS% mkl_def.lib

set LINKFLAGS=%LINKFLAGS% mkl_p3.lib

rem set LINKFLAGS=%LINKFLAGS% mkl_p4.lib

 

rem ********************************************************************

rem Uncomment these lines for IMSL libs

rem ********************************************************************

set LINKFLAGS=%LINKFLAGS% sstatd.lib sstats.lib smathd.lib smaths.lib sf90mp.lib

set LIB=%LIB%;%DFDir%/IMSL/LIB;

 

rem ********************************************************************

rem Uncomment these lines for CXML libs

rem ********************************************************************

rem set LINKFLAGS=%LINKFLAGS% cxml.lib

set LIB=%LIB%;%DFDir%/CXML/LIB

set INCLUDE=%INCLUDE%;%DFDir%/CXML/INCLUDE

/

3. Why can't I compile FORTRAN90/95 source code using the MEX command?

 

(This information has only been tested on Win32 using DVF6.x. On other platforms, mileage may vary.)

 

R13 (v6.5)

 

The Mathworks have finally updated mex.bat to handle .f90 files properly, however, they have specified compiler options that cause many Fortran9x files that worked perfectly with R11 and R12 to fail. You will need to alter your df66opts.bat file located at c:/matlab6p5/bin/win32/mexopts/ and then run >>mex -setup again

 

If you examine your mexopts.bat file, you will find variables called OPTIMFLAGS and COMPFLAGS that look like

 

set OPTIMFLAGS=/MD -Ox -DNDEBUG

 

set COMPFLAGS=/fpp:"/m /S%MATLAB%/extern/include" -c -nokeep -G5 -nologo -DMATLAB_MEX_FILE /fixed

 

set DEBUGFLAGS=/MDd -Zi

 

Replace these lines with:

 

set OPTIMFLAGS=/MT -Ox -DNDEBUG

 

set COMPFLAGS=/fpp:"/m /S%MATLAB%/extern/include" -c -nokeep -G5 -nologo -DMATLAB_MEX_FILE

 

set DEBUGFLAGS=/MTd -Zi

 

(delete /fixed and change /MD -> /MT and /MDd -> /MTd) You can see my df66opts.bat for more details.

 

R12 (v6) 

 

(This was submitted to me by Rob Brendel, I haven't tried it yet, as I don't have R12)

 

Make changes, as above, to lines 1216, 1283 and 1618 of c:/matlab/bin/mex.bat .

 

That should work just fine.

 

R11 (version 5.x)

 

Basically the problem is that mex.bat doesn't recognize files with extension .f90 as being source files, it thinks that they are object files and passes them straight to the linker. This is obviously problematic. This change will make your F90 code work fine:

 

(1) open c:/matlab/bin/mex.bat using notepad (or pico or emacs, for your *NIXers)

(2) line 2338 is:

    if ($EXTENSION =~ /(c|f|cc|cxx|cpp|for)$/i ) {

change it to

    if ($EXTENSION =~ /(c|f|cc|cxx|cpp|for|f90)$/i ) {

(notice the "|f90"?)

 

And that is all you will need! You should be able to run just fine.