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 development 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, returnity of matrix% determ = determinant of matrix% Author: Jeroen goudswaard, jcmgoudswaard@ctg.tudelft.nl
% How Many Variables IF (Nargin == 1); Verbose = 0; End; if (Verbose ~ = 0); Verbose = 1;
% 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; input; s; step 2
Now rewrite the matlab code into f90 subroutine coding And now we are at the beautiful part of it:. It will not 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 - TRANSE RETURNED' // Char (10)))! Char (10) IS A
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 :: pl Hs (*), PRHS (*)! Pointers to Input Data: POINTERS TO INPUT DATA: PRHS AND OUTPUT: PLHS, ALWAYS Take Integer * 8, To Let IT Work On 64-bit! Machines (SGI EG) 32-bit compilers Will Correct this to integer * 4, so don't worry About the! warning (s) on this during compiration. (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 :: nMat mMat,! integers that will contain the 'size' of Mat INTEGER :: nverbose! integer to copy the value of verbose to Double Precision :: Verbose, Determ! Everything That Out of Matlab Is Double Precision !! Double Precision, A Llocatable :: Mat (:, :), invmat (:)! Let's suppose we do not know the size of mat in advance! So we will use allocatable arraysif (NRHS> 2) THEN CALL MEXERRMSGTXT ('More Two Arguments ') Endiff (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 OURSLVES MMAT = MxGETM (PRHS (1))! Get the first dimension of the first input, IE mat nmat = mxgetn (PRHS (1))! Get the second dimension of the first Input, IE 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 * NMATIF (NRHS == 2) Then! Check if there is a second input: verbose pra_verbose = mxgetpr (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, EIGER! 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 INE 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 have working mexopts.sh configurations, I would be glad to public be gled to public the
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
