How to Loop in mexfunction?

Hello,

I am working on 2d array in mexfunction. And I have a couple of questions about it.

This is the example that I modify a code that I found on the website.

/*==========================================================
 * arrayMultiplier.c 
 *
 * Multiplies a 1xN matrix with a N X N matrix
 * and outputs a 1xN matrix 
 *
 * The calling syntax is:
 *    outMatrix = arrayProduct(vector, array)
 *========================================================*/
    #include "mex.h"
    /* The computational routine */
    void arrayMultiplier( double *vec, double *arr, int n, int temp ,double *z){
        temp = temp + 1;
        for (int i=0; i<n; i++) {
            for(int j=0; j<n; j++) {
            z[i] += vec[j] * arr[j+i*n]; 
            }
        }   
    }
    /* The gateway function */
    void mexFunction( int nlhs, mxArray *plhs[],
                      int nrhs, const mxArray *prhs[]) {
        double *vec;              
        double *arr;               
        int nr;
        double *outMatrix;
        int tempr;
        vec = mxGetPr(prhs[0]);
        arr = mxGetPr(prhs[1]);       /* This is where the warning appears */
        nr = mxGetScalar(prhs[2]);   /*The "n" parameter of the 1 x n output matrix*/
        tempr = mxGetScalar(prhs[3]);
        /* creating the output matrix */
        plhs[0] = mxCreateDoubleMatrix(1,nr,mxREAL);
        /* get a pointer to the real data in the output matrix */
        outMatrix = mxGetPr(plhs[0]);
        /* call the computational routine */
        arrayMultiplier(vec,arr,nr,tempr,outMatrix);
  }

After this I made .m file

clear all;
clc;

tic
cd 'C:\Users\chang\Desktop\New Folder'
mex testing.c
% ================ 1. Parameters and Constants ============================
numb = 1
A = [1,2,3,4];
B = [1,2,3,4;5,6,7,8;1,2,3,4;5,6,7,8];
AAA = testing(A, B, 4, numb);

I thought by running this code, I should get numb = 2 becasue I added temp = temp +1 in the .c file

Then I got right answer for A*B but numb is still 1 which is quite different from what I expected.

This is quite a simple version of what I am trying to do.

I am trying to use 'while' in the code but everytime I run the code it just run one time and never repeat.

while (it <= itmax) {
        it = it ++;

Then I found out that 'temp = temp + 1' doesn't work as what I expected.

Any solutions?

Thanks in advance.

P.S

Is there a way to have matrix multiplication inside mex function?

So that I can just just have A = matrixMultiply(B,C) instead of using for loop to get A.

Answer 1

[James Tursa]

You have a fundamental misunderstanding of how the C language works with pass-by-value scalar arguments. In this code:

    void arrayMultiplier( double *vec, double *arr, int n, int temp ,double *z){  /* (1) */
        temp = temp + 1;  /* (2) */
        :
    void mexFunction( int nlhs, mxArray *plhs[],
                      int nrhs, const mxArray *prhs[]) {
        int tempr;
        tempr = mxGetScalar(prhs[3]);  /* (3) */
        arrayMultiplier(vec,arr,nr,tempr,outMatrix);  /* (4) */

You pass in numb from MATLAB. In (3) above you get a copy of numb and place it in the tempr variable. From that point on you are working with the copy, not the original. So there is no way to alter the original numb by modifying tempr. Additionally, in (4) you pass tempr by value to the arrayMultiplier function. That means in (1) the temp variable is a copy of the tempr variable. You are not working with tempr directly anymore. So in (2) when you add 1 to temp, you are adding 1 to the copy and this has no effect on the original tempr variable.

There are ways to alter the MATLAB numb variable inplace in a mex routine using data pointers, but this violates the const attribute of prhs[ ] and is not advised because it can have nasty side effects.

You also asked about how to multiply matrices in mex functions without using loops. The short answer is to simply call the same BLAS library function that MATLAB uses in the background for matrix multiplies, which is dgemm for full double precision variables. You can see the interface here:

https://www.netlib.org/blas/

Some instructions for doing this in MATLAB can be found here:

https://www.mathworks.com/help/matlab/matlab_external/calling-lapack-and-blas-functions-from-mex-files.html?searchHighlight=BLAS&s_tid=srchtitle

[Chang seok Ma]

Hello James, thanks for the answer. I understand the problem. Can I ask you one more?

This is the new one I just changed it to make it simle. I think I got your point somewhat but I think I am still lost.

So basically, I want this arrayProduct to run 3 times.

#include "mex.h"

/* The computational routine */
void arrayProduct(double x, double *y, int temp, double *z, mwSize n)
{   
    mwSize i;
    while(temp<3){
    temp = temp + 1;
    /* multiply each element y by x */
    for (i=0; i<n; i++) {
        z[i] = x * y[i];
    }
    }
}

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    double multiplier;              /* input scalar */
    double *inMatrix;               /* 1xN input matrix */
    size_t ncols;                   /* size of matrix */
    double *outMatrix;              /* output matrix */
    int tempr;
    
    multiplier = mxGetScalar(prhs[0]);
    tempr = mxGetScalar(prhs[2]);
    inMatrix = mxGetPr(prhs[1]);
    ncols = mxGetN(prhs[1]);
    plhs[0] = mxCreateDoubleMatrix(1,(mwSize)ncols,mxREAL);
    outMatrix = mxGetPr(plhs[0]);
    arrayProduct(multiplier,inMatrix,tempr,outMatrix,(mwSize)ncols);
}

clc
clear
cd 'C:\Users\chang\Desktop\Test'
mex arrayProduct.c

s = 5; 
A = [1.5, 2, 9];
numb = 0;
B = arrayProduct(s,A,numb)

But if I run this I only get [7.5 10 45]. Even if 'temp' variable that I am using is the copy not the original, I think the loop still has to work because I am working with the copy. I think I am making same mistake you just pointed out.

Is there any way I could make it run 3 times not once?

[James Tursa]

You are just doing this three times:

z[i] = x * y[i];

That is, you are calculating z from y each time, but y doesn't change so you get the same answer three times.

What you need to do is have z on the right hand size so that the multiplication effect compounds each time.

So something like this:

        z[i] = x * z[i]; /* z gets calculated from the current z each time */
            :
    plhs[0] = mxDuplicateArray(prhs[0]); /* create z as a copy of the input */

[Chang seok Ma]

Thanks. Now I fully understand your comments.

[Chang seok Ma]

Sorry to keep asking but I actually have one more regarding with this problem.

#include "mex.h"
/* The computational routine */
void testing(double x, double *y, int tempr, double *z, mwSize n)
{   
    mwSize i;
    while(tempr<3){
    tempr = tempr + 1;
    /* multiply each element y by x */
    for (i=0; i<n; i++) {
        z[i] = x * y[i];
        y[i] = z[i];
    }
    }
}    
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    double multiplier;              /* input scalar */
    double *inMatrix;               /* 1xN input matrix */
    size_t ncols;                   /* size of matrix */
    double *outMatrix;              /* output matrix */
    int tempr;
    
    multiplier = mxGetScalar(prhs[0]);
    tempr = mxGetScalar(prhs[2]);
    inMatrix = mxGetPr(prhs[1]);
    ncols = mxGetN(prhs[1]);
    plhs[0] = mxCreateDoubleMatrix(1,(mwSize)ncols,mxREAL);
    outMatrix = mxGetPr(plhs[0]);
    testing(multiplier,inMatrix,tempr,outMatrix,(mwSize)ncols);
}

clc
clear all
clear mex
cd 'C:\Users\chang\Desktop\Test'
mex testing.c
s = 5;
a = [1,2,3];
tempr = 0;
b = testing(s,a,tempr)

If I run this I get b = [125,250, 375] as I expected and now I understand 'tempr' should be still 0 because the one that I am using in .c is copy version of the original variable 'tempr'

Copy of 'a' is 'inMatrix' and copy of 'inMatrix' is 'y' right?

However, now I get a = [125,250,375].

If y (in the .c file) is the copy version of [1,2,3], then shouldn't this value be still [1,2,3] because the only result I am returning is z?

[James Tursa]

You get a different 'a' because you have violated the rules and changed an input argument inplace. That is, 'a' is 'inMatrix' and 'inMatrix' is 'y', but these are not copies ... they are pointers into the original array! So changing something that the pointer points to via this line:

y[i] = z[i];

will actually change the original input variable 'a'.

[Chang seok Ma]

Just to make sure I got this concept right...

Again, this is the code that I wrote just for practice.

I want to use 'mydef' inside 'secondfn'

So this is what I intended...

1. In gateway function, Am Pm Bm are pointers saving addresses of A, B (from Matlab) respectively
2. In 'secondfn', As Bs are pointers saving addresses of Am, Bm.
3. In `mydef', this returns Cf(pointers?) using Af which points As array from 'secondfn'
4. All the m,n values I am using is just a copy of original m, n values

Then what I expected first was I would get 'Result' in matlab which shoud be equal to A.

However I got zeros matrix of 2 by 3.

So I was wondering if there is something wrong with return function in 'mydef' or it is another pointer issue I am still confusing with.

#include "mex.h"

double *mydef(double *Af, double *Bf, size_t m, size_t n, size_t p){    
    double *Cf;
    for(int j = 0; j < m*n; j++){
        Cf[j] = Af[j];
    }
return Cf;    
}

void secondfn(double *As, double *Bs, size_t m, size_t n, size_t p, double *Ps){
     double *Cs;
     Ps = mydef(As, Bs, m, n, p);
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *Am, *Bm, *Cm, *Pm; /* pointers to input & output matrices*/
    size_t m,n,p;      /* matrix dimensions */
    Am = mxGetPr(prhs[0]); /* first input matrix */
    Bm = mxGetPr(prhs[1]); /* second input matrix */
    m = mxGetM(prhs[0]);  /* dimensions of input matrices */
    p = mxGetN(prhs[0]);
    n = mxGetN(prhs[1]);
    plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
    Pm = mxGetPr(plhs[0]);
    secondfn(Am, Bm, m, n, p, Pm);
}

clear all
clear mex
cd 'C:\Users\chang\Desktop\New Folder'

mex additionalcodemine.c 
A = [1 2 3;1 2 3];
B = [3 4 5 ; 6 7 8 ; 9 10 11];
Result = additionalcodemine(A,B,2,3,3);

[James Tursa]

In this code:

    double *Cf; // this as an uninitialized pointer
    for(int j = 0; j < m*n; j++){
        Cf[j] = Af[j]; // this writes into an uninitialized address
    }

Since you are writing into an uninitialized address, anything can happen such as bad results or a MATLAB crash.

Maybe this is what you are trying to do:

void mydef(double *Af, double *Bf, size_t m, size_t n, double *Pf){    
    for(size_t j = 0; j < m*n; j++){
        Pf[j] = some appropriate calculation involving Af and Bf;
    }
}
void secondfn(double *As, double *Bs, size_t m, size_t n, size_t p, double *Ps){
     mydef(As, Bs, m, n, Ps);
}

[Chang seok Ma]

Exactly... that is what I a trying to do in the end. And as you said, uninitalized address was the problem.. your code works perfect.

If I want to also have some calculation other than `mydef' inside 'secondfn' as below,

Now I don't have any uninitialized address but I add one more matrix 'C' in the function.

I add `C' because I just need to store the result from 'mydef' and use them inside 'secondfn'

So I didn't assign anything to 'C' in 'mexFunction'

Can this be also another issue? (this also keeps crashing...)

#include "mex.h"

void mydef(double *Af, size_t m, size_t n, double *Cf){    
    for(size_t j = 0; j < m*n; j++){
        Cf[j] = Af[j];
    }
}
void secondfn(double *As, size_t m, size_t n, double *Cs, double *Ps){
     mydef(As, m, n, Cs);
     for(size_t j = 0; j < m*n; j++){
        Ps[j] = 2*Cs[j];
    }
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *A, *C, *P; /* pointers to input & output matrices*/
    size_t m,n;      /* matrix dimensions */
    A = mxGetPr(prhs[0]); /* first input matrix */
    m = mxGetM(prhs[0]);  /* dimensions of input matrices */
    n = mxGetN(prhs[0]);
    plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
    P = mxGetPr(plhs[0]);
    secondfn(A, m, n, C, P);
}

mex sample.c 
A = [1 2 3;1 2 3];
Result = sample(A);

[James Tursa]

Yes, same problem. C is an uninitialized pointer in mexFunction. So you pass it down to your functions and write into it and bad things happen. You need to allocate C first. E.g., you could put this in your mexFunction:

C = (double *) mxMalloc(m*n*sizeof(*C));
secondfn(A, m, n, C, P);
mxFree(C):

[Chang seok Ma]

Oh I get it.

I followed your comments and finally build my own code.

However I think this only works in Case1. which only have small size matrices.

If I use somewhat bigger size matrices as it is in Case2, I get different results each time I run the code.

Do you think this is just a code error or is there any other thing that I am missing?

Thank you so much for all the comments. I really appreciate it.

#include "mex.h"
void firstfn(double *Af, double *Bf, size_t m, size_t n, size_t p, double *Df){
    for(size_t k = 1; k <= n; k++){
        for(size_t i = m*(k-1); i < m*k; i++){
                size_t temp = i - m*(k-1);
                for(size_t j = 0; j < p; j++){
                    Df[i] += Af[temp+j*m] * Bf[j+p*(k-1)] ;
            }
        }
    }
}
void secondfn(double *As, double *Bs, size_t m, size_t n, size_t p, double *Ds, double *Ps){
     firstfn(As, Bs, m, n, p, Ds);
     for (size_t i = 0; i < m*n; i++){
     Ps[i] = Ds[i];
     }
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *A, *B, *D, *P; /* pointers to input & output matrices*/
    size_t m,n,p;      /* matrix dimensions */
    A = mxGetPr(prhs[0]); /* first input matrix */
    B = mxGetPr(prhs[1]); /* second input matrix */
    m = mxGetM(prhs[0]);  /* dimensions of input matrices */
    p = mxGetN(prhs[0]);
    n = mxGetN(prhs[1]);
    plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
    P = mxGetPr(plhs[0]);
    D = (double *) mxMalloc(m*n*sizeof(*D));
    //secondfn(A, B, m, n, p, D, P);
    secondfn(A, B, m, n, p, D, P);
    mxFree(D);
}

Case 1.

mex additionalcodemine2.c 
A = [1 2 3 4 ;1 2 3 4 ];
B = [3 4 5 6 ; 7 8 9 10 ; 11 12 13 14 ; 7 8 9 10];
Result = additionalcodemine2(A,B);

Case 2.

mex additionalcodemine2.c 
A= 2* ones(320,3);
B = [1/3 1/3 1/3]';
Result = additionalcodemine2(A,B);