I've written a subroutine in Fortran to handle a computationally intensive portion of my code. I want to link it to Matlab using a mex function. Here's the simplest version that gives the relevant error. First, here's the mex function.
#include <fintrf.h>
!======================================================================
#if 0
!
! example.F90
! .F90 file needs to be preprocessed to generate .for equivalent
!
#endif
!======================================================================
! Matlab command
! deltaHat = example(z)
!----------------------------------------------------------------------
! Gateway routine
SUBROUTINE mexFunction(nlhs, plhs, nrhs, prhs)
IMPLICIT NONE
! mexFunction arguments:
MWPOINTER :: plhs(*), prhs(*)
INTEGER :: nlhs, nrhs
! Function declarations:
mwPointer :: mxGetPr
mwPointer :: mxCreateDoubleMatrix
mwSize :: mxGetM, mxGetN
! Pointers to input/output mxArrays:
mwPointer :: deltaHat
mwPointer :: z
! Array information:
mwSize :: N
!------------------------------------------------------------------------
! Get the size of the input array.
N = mxGetN(prhs(1))
! Create Fortran arrays from the input argument.
z = mxGetPr(prhs(1))
! Create matrix for the return argument.
plhs(1) = mxCreateDoubleMatrix(1,24,0)
deltaHat = mxGetPr(plhs(1))
! Call the computational subroutine.
CALL example(%val(deltaHat), %val(z), N)
END SUBROUTINE
!-----------------------------------------------------------------------
! Computational routine for contraction mapping
SUBROUTINE example(deltaHat, z, N)
IMPLICIT NONE
INTEGER, PARAMETER :: dble = selected_real_kind(15)
INTEGER, PARAMETER :: J = 24
! Type delcarations for inputs
mwSize :: N
REAL(KIND=dble), DIMENSION(1,J-1), INTENT(OUT) :: deltaHat
REAL(KIND=dble), DIMENSION(N,1), INTENT(IN) :: z
! Declare local variables
REAL(KIND=dble), DIMENSION(N,J) :: num
num(:,1) = z(:,1)
DO i = 1,J-1
deltaHat(1,i) = num(i,1)
END DO
END SUBROUTINE
I then compile the mex function in Matlab using "mex example.F90." No problems there. Now here's the Matlab code and output.
clear all;
N = 1000;
J = 25;
% Generate covariates-----------------------------------------------------
% Reset random number generators
rng('default');
% Create indicators for parameter heterogeneity
z = randi(2,N,1);
z = z-ones(N,1);
deltaHat = example(z);
>> sum(deltaHat' - z(1:J-1))
ans =
-15
>> deltaHat
deltaHat =
Columns 1 through 11
1.0000 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Columns 12 through 22
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Columns 23 through 24
0.0000 0.0000
>> z(1:J-1)
ans =
1
1
0
1
1
0
0
1
1
1
0
1
1
0
1
0
0
1
1
1
1
0
1
1
The term sum(deltaHat' - z(1:J-1)) is supposed to be zero because deltaHat and z(1:J-1) are supposed to be the same vector. I also printed these vectors above in case there's something going on with variable type definitions that I'm missing (though I have checked that both are double precision). Clearly the data is not being passed to Fortran and back properly. However, if I make the following adjustment in the mex file
...
DO i = 1,J-1
deltaHat(1,i) = z(i,1)
END DO
...
then I get the right answer. Here's the corresponding Matlab output.
...
>> sum(deltaHat' - z(1:J-1))
ans =
0
...
So in this case, the data is being passed from Matlab to the Fortran subroutine and back properly. The problem seems to be in assigning z to the local variable num in the mex file. Skipping this step gives the right answer.
N
passed to SUBROUTINE example
is 1. Hence you are touching out-of-bounds of the array z
when you do assign deltaHat(1,2) = z(2,1)
change line n = mxGetN(prhs(1))
to n = mxGetM(prhs(1))
You should learn how to debug your MEX functions, so you can easily fix this kind of trivial errors.