How can I run my first Mex function?

Hi,

I am trying to run a part of my matlab code in c++.

Following this example: https://au.mathworks.com/help/matlab/matlab_external/handling-complex-data-c.html

I am working with Matlab2018a, downloaded and installed mingw-w64 c/c++ compiler from the add on menu. i did set up by

mex -setup
mex -setup cpp

After number of corrections, I have this code:

#include "mex.h"
#include <iostream>
#include <stdlib.h>
#include <complex>

using namespace std;

/* computational subroutine */
void myMexFunction(const mxArray * lambda,const mxArray * ll_n, mxArray * output,
            int N)
{
    mwSignedIndex n,s,r;
    mwSignedIndex M  = N/2-1;
    mwSignedIndex M1 = M+1;
    mwSignedIndex c,d,rmin,rmax;

      mexPrintf("M%d\n", M);

          /* get pointers to the  arrays */
          mxDouble        * lam = mxGetDoubles(lambda);
          mxComplexDouble * lln = mxGetComplexDoubles(ll_n);
          mxComplexDouble * out = mxGetComplexDoubles(output);
          mexPrintf("inside the function\n");

          /* perform the nestedloop op */
            for(n = -M; n < M; ++n){
                mexPrintf("inside the first loop\n");
                c = min(M-n,M);
                mexPrintf("c %d\n", c);
                d = max(-M,-M-n);
                mexPrintf("d %d\n", d);
                for(s = -M; s < M; ++s){
                    mexPrintf("inside the second loop\n");
                    rmin=max(s-M,d);
                    rmax=min(s+M,c);
                    mexPrintf("rmin %d\n", rmin);
                    mexPrintf("rmax %d\n", rmax);
                    for(r= rmin; r<rmax; ++r){
                        mexPrintf("inside the third loop\n");

                         out[1].real=1;
                         out[2].imag=1;
  //                       out[N*(s+M)+(n+M1)].real =
  //                       out[N*(s+M)+(n+M1)].real + (lln[s-r+M1].real * lln[n+r+M1].real - lln[s-r+M1].imag * lln[n+r+M1].imag)*lam[r+M1];
  //                       out[N*(s+M)+(n+M1)].imag =
  //                       out[N*(s+M)+(n+M1)].imag + (lln[s-r+M1].real * lln[n+r+M1].imag + lln[s-r+M1].imag * lln[n+r+M1].real)*lam[r+M1];
                    }
                }
            }

}

/* The gateway routine. */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
    int N = (int) mxGetScalar(prhs[3]);;

      /* coppy array and set the output pointer to it */

      plhs[0] = mxDuplicateArray(prhs[2]);

      /* call the C subroutine */
      myMexFunction(prhs[0], prhs[1], plhs[0], N);

      return;
  }

When I call this function from matlab using:

mex myMexFunction.cpp -R2018a
Lsn_=myMexFunction(lambda,ll_n,Lsn,N);

I can see that the function goes into all loops but if I assign any value to out[] matlab crashes. I cannot assign even a fixed value to out[1] so it is not index related.

What can i do to get it work?

Thanks

Answer 1

[James Tursa]

Some problems:

1)

      /* coppy array and set the output pointer to it */
      plhs[0] = prhs[2];

The above line does not copy the input array to the output array. It only copies a pointer to a pointer. Furthermore, you cannot copy input prhs[] pointers to output plhs[] pointers directly. This violates the rules and can cause a MATLAB crash. Probably what you want instead is a deep copy:

      plhs[0] = mxDuplicateArray(prhs[2]);

2)

out[s+M1,n+M1].real = out[s+M1,n+M1].real + etc.

The above line does not do the indexing like you think it is doing. This is not the way to index into a 2D array in C/C++. What you are actually doing in the above line is evaluating the result of a comma operator and using that as a single index. Since the result of a comma operator is the value of the 2nd expression, the above line is the equivalent of

out[n+M1].real = out[n+M1].real + etc.

You will need to manually calculate the equivalent linear index corresponding to the 2D element you are trying to update, and use that linear index instead of what you are doing.

3)

    int N = prhs[3];

The above line copies an mxArray * type into an integer. This does not get the value into N. To do that you need to do this instead:

    int N = (int) mxGetScalar(prhs[3]);

There are probably other errors too, but I can't read your screen shot above. Please delete the screen shot and post the error messages as text instead.

[ttopal]

Thank you James, errors 1 and 3 are very clear, for error 2 i will do some reading do you think i cant use indexing like out[][] either? Also the main problem looks like compiler cant find mex functions such as mxGetDoubles and mxComplexDouble am I wrong? I have updated the code and error messages.

[James Tursa]

4)

It can't find the mxGetDoubles and mxGetComplexDoubles because you are not linking with the correct API library. The default library is -R2017b, which doesn't have those functions. You need to explicitly use the -R2018a compile option to link with the new functions. E.g.,

mex myMexFunction.cpp -R2018a

5)

Apply const to the prhs[] arguments of your routine to get rid of the conversion errors. E.g.,

void myMexFunction(const mxArray * lambda, const mxArray * ll_n, mxArray * output, int N)

2a)

You still don't have the indexing correct. Your current code is this:

out[s+M1][n+M1].real = out[s+M1][n+M1].real + etc;

The variable "out" is a pointer to mxComplexDouble, it is not a pointer to pointer to mxComplexDouble. As such, "out" can only be dereferenced once, not twice like you are doing.

I am not totally sure what the dimensions of your variables are, but it looks like maybe the matrices are NxN. Assuming this is the case, your dereferencing should probably look something like this:

out[(s+M1)+N*(n+M1-1)].real = out[(s+M1)+N*(n+M1-1)].real + etc;

or perhaps this

out[N*(s+M1-1)+(n+M1)].real = out[N*(s+M1-1)+(n+M1)].real + etc;

Or maybe some similar variation of these. The correct one will depend on the order you want to access the data, which again is not entirely clear to me from what you have provided. Basically, you have to manually calculate the equivalent linear index into the 2D matrix.

That being said, your indexing looks suspect to me. Ordinarily I would expect the indexing to notionally start with the [0] element. However, in your code the first element accessed is when s=-M and n=-M, which gives you this:

out[N*(-M+M1-1)+(-M+M1)]
out[N*(-M+M+1-1)+(-M+M+1)]
out[N*(0)+(1)]
out[1]

So your algorithm starts with the 2nd element of the variable, and not the 1st element of the variable. You might double check your indexing formulas to see if they traverse the matrix correctly.

[ttopal]

Thank you very much for your help and time James, i got it to compile, but only to crash. I am updating the post and error again.

[ttopal]

I figured the crash as indices needed to be transformed into mwSize and size_t. But this time i cant get it to go into the loops. i am updating the question and error once again.

[James Tursa]

mwSize is probably a size_t in the background, which is unsigned. Setting an unsigned variable to a negative value will automatically wrap the value around to a positive value in a modulo sense. For the -M value you are starting with, this means your for-loop indexing value starts with a very large positive number and never satisfies n < M requirement. Hence the for-loop never runs. If you want to use negative values with your for-loop index, you will need to use signed integers for your indexing variable.

[ttopal]

James, I really appreciate your help, I am learning as I am correcting this errors one by one. For the signed integers,mwSignedIndex works. Now all the loops are working. But something wrong with out matrix. Even assigning out[1].real=1: crashes matlab. As for the indexes you are right I am starting from 1, the outcome of those rmin and rmax looks like this:

M = 1  
rmin =
       0    -1    -1
       0    -1    -1
       0     0     0
rmax =  
       0     0     0
       1     1     0
       1     1     0
M=2
rmin =
     0    -1    -2    -2    -2
     0    -1    -2    -2    -2
     0    -1    -2    -2    -2
     0    -1    -1    -1    -1
     0     0     0     0     0
rmax =
     0     0     0     0     0
     1     1     1     1     0
     2     2     2     1     0
     2     2     2     1     0
     2     2     2     1     0

any idea why i can not assign a value to out[] matrix? I have updated the code again.

[ttopal]

I have finally figured it. The initial array that i duplicated (zeros matrix from matlab) was not complex type, hence the crash! I was doing it so that I didnt have to write to loops to make a zeros matrix in c++. One final question is how i can make out[][] indexing work?

[James Tursa]

"how i can make out[][] indexing work?"

You can make the [][] indexing syntax work, sort of, but IMO it usually isn't worth the trouble and I prefer simple linear indexing instead. MATLAB data storage ordering in memory is column-based, but native C/C++ arrays are stored in memory row-based. So if you use the [][] syntax on a native C/C++ 2D array it will naturally be interpreted as [row][column]. But if you attempt to use [][] on a MATLAB mxArray 2D matrix memory block it will be interpreted as [column][row] ... i.e. basically the transpose of the C/C++ interpretation. So it can get a bit confusing to the reader.

I wrote a detailed mex example of all of this on Answers some time ago, but I can't find the post at the moment. In lieu of that, here is a short example. The code simply prints out a full double complex 2D matrix to the screen:

#include "mex.h"
#define PR(i,j) pr[(i)+M*(j)]

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mxComplexDouble *pr;
    mxComplexDouble **ptr;
    size_t M, N, i, j;

      if( nrhs != 1 || !mxIsDouble(prhs[0]) || !mxIsComplex(prhs[0]) ||
          mxIsSparse(prhs[0]) || mxGetNumberOfDimensions(prhs[0]) != 2 ) {
          mexErrMsgTxt("Need one full double complex 2D input");
      }
      if( nlhs ) {
          mexErrMsgTxt("Too many outputs");
      }
      pr = mxGetComplexDoubles(prhs[0]);
      M = mxGetM(prhs[0]);
      N = mxGetN(prhs[0]);

      /* Linear Indexing */
      mexPrintf("Linear Indexing\n");
      for( i=0; i<M; i++ ) {
          mexPrintf("Row %2d = ",i+1);
          for( j=0; j<N; j++ ) {
              mexPrintf("%g + %gi ,",pr[i+M*j].real,pr[i+M*j].imag);
          }
          mexPrintf("\n");
      }
      mexPrintf("\n");

      /* Macro Indexing */
      mexPrintf("Macro Indexing\n");
      for( i=0; i<M; i++ ) {
          mexPrintf("Row %2d = ",i+1);
          for( j=0; j<N; j++ ) {
              mexPrintf("%g + %gi ,",PR(i,j).real,PR(i,j).imag);
          }
          mexPrintf("\n");
      }
      mexPrintf("\n");

      /* [][] Indexing */
      mexPrintf("[][] Indexing\n");
      ptr = (mxComplexDouble **) mxMalloc( N * sizeof(*ptr) ); /* Allocate the column pointer array */
      /* Fill in the column pointer array */
      ptr[0] = pr;
      for( j=1; j<N; j++ ) {
          ptr[j] = ptr[j-1] + M;
      }
      for( i=0; i<M; i++ ) {
          mexPrintf("Row %2d = ",i+1);
          for( j=0; j<N; j++ ) {
              mexPrintf("%g + %gi ,",ptr[j][i].real,ptr[j][i].imag); /* Note the order of [j][i] */
          }
          mexPrintf("\n");
      }
      mexPrintf("\n");
      mxFree(ptr); /* Free the column pointer array */
  }

And a sample run:

>> mex indexing2D.c -R2018a
Building with 'Microsoft Visual C++ 2013 (C)'.
MEX completed successfully.
>> x = reshape(1:12,3,4) - reshape(1:12,3,4)*i
x =
   1.0000 - 1.0000i   4.0000 - 4.0000i   7.0000 - 7.0000i  10.0000 -10.0000i
   2.0000 - 2.0000i   5.0000 - 5.0000i   8.0000 - 8.0000i  11.0000 -11.0000i
   3.0000 - 3.0000i   6.0000 - 6.0000i   9.0000 - 9.0000i  12.0000 -12.0000i
>> indexing2D(x)
Linear Indexing
Row  1 = 1 + -1i ,4 + -4i ,7 + -7i ,10 + -10i ,
Row  2 = 2 + -2i ,5 + -5i ,8 + -8i ,11 + -11i ,
Row  3 = 3 + -3i ,6 + -6i ,9 + -9i ,12 + -12i ,

Macro Indexing
Row  1 = 1 + -1i ,4 + -4i ,7 + -7i ,10 + -10i ,
Row  2 = 2 + -2i ,5 + -5i ,8 + -8i ,11 + -11i ,
Row  3 = 3 + -3i ,6 + -6i ,9 + -9i ,12 + -12i ,

[][] Indexing
Row  1 = 1 + -1i ,4 + -4i ,7 + -7i ,10 + -10i ,
Row  2 = 2 + -2i ,5 + -5i ,8 + -8i ,11 + -11i ,
Row  3 = 3 + -3i ,6 + -6i ,9 + -9i ,12 + -12i ,