Fortran Mex File - Incorrect Values Passed Back to Matlab

Question

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.

Answer 1

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.