Multidimensional matrix multiplication with mex

Hi

I am trying to do matrix multiplication via Mex.

So fare i have, using this as guide :

#include "mex.h"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

      float *x;
      float *y;
      double *z;

      mwSize rFirst, cFirst, cSecond;

      x = mxGetData(prhs[0]);
      y = mxGetData(prhs[1]);

      rFirst = mxGetM(prhs[0]);
      cFirst = mxGetN(prhs[0]);
      cSecond = mxGetN(prhs[1]);

      plhs[0] = mxCreateDoubleMatrix(rFirst, cFirst, mxREAL);

      z = mxGetPr(plhs[0]);

      int i;
      int j;
      int k;

      for(k = 0; k <= cSecond - 1 ; k++) {
          for(i = 0; i <= rFirst - 1 ; i++) {

              z[i][k] = 0;

              for(j = 0; j <= cFirst - 1; j++) {

              z[i][k] = z[i][k] + x[i][j]*y[j][k];

              }
          }
      } 
  }

When i try and implement this i get the error:

    > mex matrixMultiplication.c
  error: subscripted value is neither array nor pointer

I then tried changing the script to:

    float **x;
    float **y;
    double **z;

But that results in:

    > mex matrixMultiplication.c
  warning: assignment from incompatible pointer type

Does anyone know how to fix this?

Here is an extra question: I tried importing a mex file from File Exchange: mTimesx. I unzipped the folder and placed the content in my working directory. I then tried initializing it (Or what you call it).

    > mex mTimesx.c
  undefined reference to 'foobar'

I got this error many many times. How do you import a C-mex file correctly?

Hope someone knows the answers to these things :) Thanks in advance!!

Answer 1

[James Tursa]

Your basic problem is a misunderstanding of how pointers work in C. Take your z variable for instance:

    double *z;
       :
    z = mxGetPr(plhs[0]);
               :
            z[i][k] = 0;

The variable z is a "pointer to double". You need to think of it as pointing into a linear array of double values. Since it points to a double, the result of dereferencing it once is a double. E.g., the expression z[i] by itself is a double ... it is not another pointer. So now when you try to apply another dereference with the [k] part, the compiler complains because z[i] is not a pointer so it can't be dereferenced. You can only dereference z once. In other words, you can't use two indexes with z the way you have z defined. You would need to manually calculate the location in memory and use that as a linear subscript. E.g., this would work:

            z[i+k*rFirst] = 0;

Same comment applies to your x and y variables. You would need to manually calculate an equivalent linear index to use as a single subscript. E.g.

            z[i+k*rFirst] = z[i+k*rFirst] + x[i+j*rFirst]*y[j+k*cFirst];

There are ways to use [i][k] subscripting, but you have to define your pointers differently and do some prep work up front to get things to work properly. I will not include that code here but will refer you to this link:

https://www.mathworks.com/matlabcentral/answers/309434-matlab-crashing-when-evaluating-mex

Additionally, you don't have your output dimensions correct. The line

    plhs[0] = mxCreateDoubleMatrix(rFirst, cFirst, mxREAL);

should be

    plhs[0] = mxCreateDoubleMatrix(rFirst, cSecond, mxREAL);

Also, I will point out that your code will likely crash MATLAB if the inputs are not exactly as expected. I.e., check that there are two inputs and they are each single class matrices and that their inner dimensions agree and that they are not complex. If you were to allow double class inputs, you would need to check that they were not sparse as well. You should put in code to check for those things.

Your second problem is that the author of the MTIMESX submission on the FEX has neglected for some time to get it updated for the newer MATLAB versions and for 64-bit machines. He needs to somehow get access to a newer MATLAB version with a C compiler and find some time to work on that. What MATLAB version are you using, and are you running 32-bit or 64-bit?

[Michael Madelaire]

Thanks for the answer. You are indeed correct. I have only had a minor introduction to C.

In regard to the first answer: I made the changes you pointed out, thanks for catching the output mistake :) Before implementing the input check, I thought it necessary to check the performance of the C-mex. (I will post the c-mex script i used at the end.) I came to the conclusion that my matrixMultiplication.c is VERY much slower than just using *. I may be naive, but shouldn't C-mex speed up the process? Or maybe it does but my programming skills lag abit ;)

 a = randi(5, 10000, 10000);
 b = randi(5, 10000, 10);

 tic; c = a*b; toc
Elapsed time is 0.258634 seconds.

 tic; matrixMultiplication(a, b); toc
Elapsed time is 4.673416 seconds.

Is the computational time different because the c-Mex is written inefficiently?

In regard to my second question. I am running Matlab R2016a 64bit on a linux computer.

if true (same trouble as last time, i think it is the pound)
  #include "mex.h"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

    double *x;
    double *y;
    double *z;

    mwSize rFirst, cFirst, cSecond;

    x = mxGetData(prhs[0]);
    y = mxGetData(prhs[1]);

    rFirst = mxGetM(prhs[0]);
    cFirst = mxGetN(prhs[0]);
    cSecond = mxGetN(prhs[1]);

    plhs[0] = mxCreateDoubleMatrix(rFirst, cSecond, mxREAL);

    z = mxGetPr(plhs[0]);

    int i;
    int j;
    int k;

    for(k = 0; k <= cSecond - 1 ; k++) {
        for(i = 0; i <= rFirst - 1 ; i++) {

            z[i+k*rFirst] = 0;

            for(j = 0; j <= cFirst - 1; j++) {

            z[i+k*rFirst] = z[i+k*rFirst] + x[i+j*rFirst]*y[j+k*cFirst];

            }
        }
    } 
}

[James Tursa]

"... I came to the conclusion that my matrixMultiplication.c is VERY much slower than just using *. I may be naive, but shouldn't C-mex speed up the process? ..."

You came to the correct conclusion. MATLAB uses a very highly optimized 3rd party library to do linear algebra calculations. In this case, the library is called BLAS (Basic Linear Algebra Subroutines). There is also a library call LAPACK that gets used for more complicated stuff like matrix decomposition and the like. Bottom line is your handwritten code is competing against code written by a team of professional experts that is optimized for the architecture of your machine, some parts of which may be written in assembly, and in many cases is multi-threaded. There are things you could do to improve your code (e.g., to speed up your indexing), but even with such improvements you have no chance to beat the BLAS for this! :)

[Michael Madelaire]

But the use of mtimex would be faster, correct?

[James Tursa]

MTIMESX calls the same BLAS library directly that MATLAB uses (it links directly to it), so it matches MATLAB exactly for speed and results in simple matrix multiply cases. Where MTIMESX can sometimes beat MATLAB for speed is with multi-dimensional cases where the first two dimensions are regarded as 2D matrix planes, and with some large complex and sparse cases (although MATLAB has gotten better with these over the years). MTIMESX can do these cases without the intermediate data copying that is necessary to do the same calculations at the m-file level. So, in a sense, MTIMESX isn't doing any of the matrix multiplying any faster ... it is just avoiding the time associated with unnecessary data copying. Results are highly dependent on what exactly you are doing, what compiler you are using (and whether it supports OpenMP or not), and the sizes involved. E.g., just using Microsoft Visual Studio instead of a generic compiler like lcc can make a huge difference in speed.

[Michael Madelaire]

I am having a hard time getting mTimesX to work. I have tried the manual compile. But I get an error i cannot identify.

    Unknown file extension ' '

I therefore cannot test it myself. But I am working with some very big 2D matrices (170 GB). As I understand mTimesX could speed up multiplication of matrices of this size? Or is the implementation unnecessary?

[James Tursa]

If you are doing just the straightforward calculation (large 2D matrix) * (large 2D matrix), MTIMESX will not be able to help you.

[Michael Madelaire]

Thank you very much for the help! Although I didn't find some black magic to speed up my simulations :)