Reducing run time of a numerical calculation using a mex file in Matlab

I wrote a Matlab code that involves doing a numeric calculation (relaxation), but it is quite slow. I learned of the possibility of using a mex file to run a C code and integrate it into Matlab, so I was thinking of doing the numerical calculation (which is relatively simple but involves loops and takes time) in C, and the rest (before and after) in Matlab.

The part of my Matlab code where the calculation is done:

% evolution of the potentials %
% note : for the index directions with periodic boundary conditions: index=mod(index-1,L)+1 . for index=index+1 it is mod(index,L)+1 , and for index=index-1 it is mod(index-2,L)+1 %          
for i_t=1:max_relaxation_iterations
    for q=1:length(i_eff_V_bounded) % this is set instead of running i=2:(L-1), j=1:L , k=1:L and ending up going over sites that are 0 in our effective system %
    i=i_eff_V_bounded(q);
    j=j_eff_V_bounded(q);
    k=k_eff_V_bounded(q);
    V0=V(i,j,k);    
    V1=( V(i+1,j,k)+V(i-1,j,k)+V(i,mod(j,L)+1,k)+V(i,mod(j-2,L)+1,k)+V(i,j,mod(k,L)+1)+V(i,j,mod(k-2,L)+1) )/( system(i+1,j,k)+system(i-1,j,k)+system(i,mod(j,L)+1,k)+system(i,mod(j-2,L)+1,k)+system(i,j,mod(k,L)+1)+system(i,j,mod(k-2,L)+1) ); % evolving the potential as the average of its occupied neighbors %
    V(i,j,k)=V0+(V1-V0)*over_relaxation_factor; % evolving the potentials in time with the over relaxation factor %
    delta_V_rms(i_t)=delta_V_rms(i_t)+(V1-V0)^2; % for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg %
    delta_V_abs(i_t)=delta_V_abs(i_t)+abs(V1-V0); % for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg %                
    delta_V_max(i_t)=max(abs(V1-V0),delta_V_max(i_t)); % for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg % 
    end  
end

So in C it should be something like:

#include <stdio.h>
    
                            
int mod(int x,int N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}


double d_abs(double x) /* a function for the absolute value operator */
{
if x<0
    {
    return -x;
    }
else
    {
    return x;
    }
}

                    
double max(double x,double y) /* a function for the max operator */
{
if x>y
    {
    return x;
    }
else
    {
    return y;
    }
}


/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */          
void potentials_evolution(int max_relax_iters,int N_eff_occ_sites,int i_eff_V_bounded[],int j_eff_V_bounded[],int k_eff_V_bounded[],int system[][][],over_relax_fact,double V[][][],double delta_V_rms[],double delta_V_abs[],double delta_V_max[])
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */ 
        }
    }
}

And so in Matlab I will replace the part of my Matlab code shown above with something like:

potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);

How do I implement this? I tried looking around (in places such as here: https://www.mathworks.com/help/matlab/matlab_external/standalone-example.html) for a simple way to do it but I couldn't figure out how to properly do it.

Note 1: This numeric calculation is done not just once but many times for different systems that are generated randomly (there is a `for` loop going over the different systems).

Note 2: I've never used mex files and my C is quite rusty.

EDIT: After hours of trying to figure how to build the mex file, I came up with this:

#include "mex.h"

                            
mwSize mod(mwSize x,mwSize N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}


double d_abs(double x) /* a function for the absolute value operator */
{
if(x<0)
    {
    return -x;
    }
else
    {
    return x;
    }
}

                    
double max(double x,double y) /* a function for the max operator */
{
if(x>y)
    {
    return x;
    }
else
    {
    return y;
    }
}

/* in matlab i will invoke potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max); */

/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */          
void potentials_evolution(mwSize L,mwSize max_relax_iters,mwSize N_eff_occ_sites,mwSize *i_eff_V_bounded,mwSize *j_eff_V_bounded,mwSize *k_eff_V_bounded,mwSize ***system,double over_relax_fact,double ***V,double *delta_V_rms,double *delta_V_abs,double *delta_V_max)
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */ 
        }
    }
}


/* in matlab: [V,delta_V_rms,delta_V_abs,delta_V_max] = potentials_evolution(i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V_0,delta_V_rms_0,delta_V_abs_0,delta_V_max_0) */
void mexFunction(int nlhs,mxArray *plhs,int nrhs,const mxArray *prhs[])
{
size_t L, *dim_V, max_relax_iters, N_eff_occ_sites;

size_t *i_eff_V_bounded, *j_eff_V_bounded, *k_eff_V_bounded;   
size_t ***system;       
double over_relax_fact;        
double ***V, ***V_0;        
double *delta_V_rms, *delta_V_abs, *delta_V_max, *delta_V_rms_0, *delta_V_abs_0, *delta_V_max_0;        


/* additional input for the calculation */
dim_V = mxGetDimensions(prhs[5]);
L = dim_V[0];
max_relax_iters = mxGetN(prhs[6]);
N_eff_occ_sites = mxGetN(prhs[0]);

/* input arguments */
i_eff_V_bounded = mxGetInt32s(prhs[0]);
j_eff_V_bounded = mxGetInt32s(prhs[1]);
k_eff_V_bounded = mxGetInt32s(prhs[2]);
system = mxGetInt32s(prhs[3]);
over_relax_fact = mxGetScalar(prhs[4]);
V_0 = mxGetDoubles(prhs[5]);
delta_V_rms_0 = mxGetDoubles(prhs[6]);
delta_V_abs_0 = mxGetDoubles(prhs[7]);
delta_V_max_0 = mxGetDoubles(prhs[8]);        

/* output arguments */
plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[5]),mxGetDimensions(prhs[5]),mxDOUBLE_CLASS,mxREAL);
V = mxGetDoubles(plhs[0]);
plhs[1] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_rms = mxGetDoubles(plhs[1]);
plhs[2] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_abs = mxGetDoubles(plhs[2]);  
plhs[3] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_max = mxGetDoubles(plhs[3]);  
        
/* performing the calculation */
potentials_evolution((mwSize)L,(mwSize)max_relax_iters,(mwSize)N_eff_occ_sites,(mwSize)i_eff_V_bounded,(mwSize)j_eff_V_bounded,(mwSize)k_eff_V_bounded,(mwSize)system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);
}

Is this correct syntax/form-wise and does this match my Matlab code (shown above)? Do I need something else for it to run?

Answer 1

[Jorik]

If you have access to MATLAB Coder, it can let you generate a C-MEX implementation for your MATLAB code. If you use the Coder App, it will guide you through the steps to generate a MEX-function and validate it runs correctly.

If you'd like to write it manually, you'll need to learn the basics of the MEX API. I would start by copying and modifying a simple function like timestwo.c:

edit(fullfile(matlabroot, 'extern', 'examples', 'refbook', 'timestwo.c'))

You can lookup the API functions in the documentation to better understand the MEX and MATRIX API, e.g.:

https://www.mathworks.com/help/releases/R2019a/matlab/apiref/mexfunction.html

and

https://www.mathworks.com/help/releases/R2019a/matlab/apiref/mxgetm.html

[tensorisation]

I have Matlab coder (although I've never used it), but that is for turning Matlab code into C code no? I've already wrote what my C code will be based on the part of my Matlab code that I want to run in C (that I also provided here).

I tried looking at this timestwo.c example but I cannot understand from it how to properly implement this for my case.

[tensorisation]

After hours of trying to figure how to build the mex file, I came up with this:

#include "mex.h"

                            
mwSize mod(mwSize x,mwSize N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}


double d_abs(double x) /* a function for the absolute value operator */
{
if x<0
    {
    return -x;
    }
else
    {
    return x;
    }
}

                    
double max(double x,double y) /* a function for the max operator */
{
if x>y
    {
    return x;
    }
else
    {
    return y;
    }
}

/* in matlab i will invoke potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max); */

/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */          
void potentials_evolution(mwSize max_relax_iters,mwSize N_eff_occ_sites,mwSize *i_eff_V_bounded,mwSize *j_eff_V_bounded,mwSize *k_eff_V_bounded,mwSize ***system,double over_relax_fact,double ***V,double *delta_V_rms,double *delta_V_abs,double *delta_V_max)
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */ 
        }
    }
}


/* in matlab: [V,delta_V_rms,delta_V_abs,delta_V_max] = potentials_evolution(i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V_0,delta_V_rms_0,delta_V_abs_0,delta_V_max_0) */
void mexFunction(int nlhs,mxArray *plhs,int nrhs,const mxArray *prhs[])
{
size_t max_relax_iters, N_eff_occ_sites;

size_t *i_eff_V_bounded, *j_eff_V_bounded, *k_eff_V_bounded;   
size_t ***system;       
double over_relax_fact;        
double ***V, ***V_0;        
double *delta_V_rms, *delta_V_abs, *delta_V_max, *delta_V_rms_0, *delta_V_abs_0, *delta_V_max_0;        


/* additional input for the calculation */
max_relax_iters = mxGetN(prhs[6]);
N_eff_occ_sites = mxGetN(prhs[0]);

/* input arguments */
i_eff_V_bounded = mxGetInt32s(prhs[0]);
j_eff_V_bounded = mxGetInt32s(prhs[1]);
k_eff_V_bounded = mxGetInt32s(prhs[2]);
system = mxGetInt32s(prhs[3]);
over_relax_fact = mxGetScalar(prhs[4]);
V_0 = mxGetDoubles(prhs[5]);
delta_V_rms_0 = mxGetDoubles(prhs[6]);
delta_V_abs_0 = mxGetDoubles(prhs[7]);
delta_V_max_0 = mxGetDoubles(prhs[8]);        

/* output arguments */
plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[5]),mxGetDimensions(prhs[5]),mxDOUBLE_CLASS,mxREAL);
V = mxGetDoubles(plhs[0]);
plhs[1] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_rms = mxGetDoubles(plhs[1]);
plhs[2] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_abs = mxGetDoubles(plhs[2]);  
plhs[3] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_max = mxGetDoubles(plhs[3]);  
        
/* performing the calculation */
potentials_evolution((mwSize)max_relax_iters,(mwSize)N_eff_occ_sites,(mwSize)i_eff_V_bounded,(mwSize)j_eff_V_bounded,(mwSize)k_eff_V_bounded,(mwSize)system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);
}

Is this correct syntax/form-wise and does this match my Matlab code (shown above)? Do I need something else for it to run?

[Jorik]

The nice thing about MATLAB Coder, is that besides generating standalone C code, it allows you to create a MEX-Function with the generated C code to verify in MATLAB with test stimuli if it indeed gives the same results as the MATLAB code.

Did you try compiling this in MATLAB? I think I see a couple of violations of ANSI C syntax. For instance:

if x>y

in C this should be

if (x>y)

After fixing this, the next compiler error shows me:

error: use of undeclared identifier 'L'

BTW, I'm surprised to see a custom implementation for d_abs as I would expect fabs() to be good enough?

Anyway, if you're trying to get a faster implementation of your function by leveraging C code, I would personally always start with MATLAB Coder (if available) and then maybe tweak it a bit.

But of course if you'd like to know how to write C MEX-Functions, this could be a good candidate.

I would recommend (as this is your first MEX-Function), to attach a debugger when you run it the first time when building it succeeds and step through the code. It's quite easy to make MATLAB crash with MEX-Functions that poke at wrong memory locations ;-). In general it's easy to make small mistakes with "big" effects in C.

[tensorisation]

1.) I'm quite rusty in C - I fixed the if. I also forgot to add an input of L in the calculation of void potentials_evolution - added it now.

I now have:

#include "mex.h"

                            
mwSize mod(mwSize x,mwSize N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}


double d_abs(double x) /* a function for the absolute value operator */
{
if(x<0)
    {
    return -x;
    }
else
    {
    return x;
    }
}

                    
double max(double x,double y) /* a function for the max operator */
{
if(x>y)
    {
    return x;
    }
else
    {
    return y;
    }
}

/* in matlab i will invoke potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max); */

/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */          
void potentials_evolution(mwSize L,mwSize max_relax_iters,mwSize N_eff_occ_sites,mwSize *i_eff_V_bounded,mwSize *j_eff_V_bounded,mwSize *k_eff_V_bounded,mwSize ***system,double over_relax_fact,double ***V,double *delta_V_rms,double *delta_V_abs,double *delta_V_max)
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */ 
        }
    }
}


/* in matlab: [V,delta_V_rms,delta_V_abs,delta_V_max] = potentials_evolution(i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V_0,delta_V_rms_0,delta_V_abs_0,delta_V_max_0) */
void mexFunction(int nlhs,mxArray *plhs,int nrhs,const mxArray *prhs[])
{
size_t L, *dim_V, max_relax_iters, N_eff_occ_sites;

size_t *i_eff_V_bounded, *j_eff_V_bounded, *k_eff_V_bounded;   
size_t ***system;       
double over_relax_fact;        
double ***V, ***V_0;        
double *delta_V_rms, *delta_V_abs, *delta_V_max, *delta_V_rms_0, *delta_V_abs_0, *delta_V_max_0;        


/* additional input for the calculation */
dim_V = mxGetDimensions(prhs[5]);
L = dim_V[0];
max_relax_iters = mxGetN(prhs[6]);
N_eff_occ_sites = mxGetN(prhs[0]);

/* input arguments */
i_eff_V_bounded = mxGetInt32s(prhs[0]);
j_eff_V_bounded = mxGetInt32s(prhs[1]);
k_eff_V_bounded = mxGetInt32s(prhs[2]);
system = mxGetInt32s(prhs[3]);
over_relax_fact = mxGetScalar(prhs[4]);
V_0 = mxGetDoubles(prhs[5]);
delta_V_rms_0 = mxGetDoubles(prhs[6]);
delta_V_abs_0 = mxGetDoubles(prhs[7]);
delta_V_max_0 = mxGetDoubles(prhs[8]);        

/* output arguments */
plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[5]),mxGetDimensions(prhs[5]),mxDOUBLE_CLASS,mxREAL);
V = mxGetDoubles(plhs[0]);
plhs[1] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_rms = mxGetDoubles(plhs[1]);
plhs[2] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_abs = mxGetDoubles(plhs[2]);  
plhs[3] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_max = mxGetDoubles(plhs[3]);  
        
/* performing the calculation */
potentials_evolution((mwSize)L,(mwSize)max_relax_iters,(mwSize)N_eff_occ_sites,(mwSize)i_eff_V_bounded,(mwSize)j_eff_V_bounded,(mwSize)k_eff_V_bounded,(mwSize)system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);
}

2.) In regards to fabs, if I want to use it, I need to put:

#include <math.h>

in the file as well?

Won't it run faster if i use my own d_abs instead of fabs?

3.) in regards to the lines I just added (for the sake of L):

dim_V = mxGetDimensions(prhs[5]);
L = dim_V[0];

Is it possible to get L straight without needing to define and use dim_V by doing something like:

L = mxGetDimensions(prhs[5])[0];

?

4.) In my Matlab code, the 3 dimensional array (L-L-L array) system is an array of 0's and 1's that at first is generated randomly and later on is turned into a logical array (all this is still before doing the calculation with potentials_evolution). In Matlab, summing elements of system when it's a logical array is not a problem (for example if system(1,1,1) is 1 (true) and system (1,1,2) is 1 (true), then system(1,1,1)+system(1,1,2) will give 2), how do I handle this in C to work the same way?

Is there even such a thing as a logical array in C? Assuming system is a logical Array from Matlab, what class do I make it in C and how do I handle it properly? If you look at the code above you can see that I used mxGetInt32s for it, but it's probally incorrent and inefficient.

Should I use mxGetLogicals for it instead? But then again - what will happen when I add elements of system together? I want it to function like in Matlab - adding them like numeric 1's and 0's.

[Jorik]

2)

Yes you need to include math.h to use fabs(). I would expect the fabs() implementation to be at least as fast, e.g. it may very well be implemented as a macro, so that avoids the overhead of a function call.

3)

Yes I think that should be possible.

4)

There is no bool/logical type in ANSI C (89/90), so you'll have to use an integer. You can then decide, short int, int, long int, etc. balancing between maximum value ranges and possibly memory usage.

BTW, did you already try MEX'ing your source code? I still can't get it to compile, because of more errors (and warnings), e.g. the mexFunction function declaration is different in mex.h than in your code.

[tensorisation]

1.) so just to be clear, this is a correct syntax in C?

L = mxGetDimensions(prhs[5])[0];

2.) what about mxGetLogicals and mxLOGICAL_CLASS then? And if in Matlab my array is logical, will it work fine if I use mxGetUint8s straight on it as is or do I have to turn the array to be class uint8 in matlab beforehand for it to work?

3.) What exactly do you mean by " the mexFunction function declaration is different in mex.h than in your code"?

[Jorik]

1) I don't think I'd be helping you a lot if I would fully answer this question. Please look up the reference page, it has links to 6 examples that use this API to get the dimensions. Maybe you want to make a copy of one of the examples, change the code, compile it and inspect how it works. Just as a hint, this line is in the first example from the reference page:

if (temp[x]> ((mxGetDimensions(prhs[0]))[x]) ){

That looks quite similar, so in any case you definitely don't need to use a temporary variable.

2) Did you look at the mxGetLogicals function reference page? It links to an example that shows how you can use it. No, you should not use mxGetUint8s on an mxArray that is of class mxLOGICAL_CLASS. The function reference page for mxGetUint8s mentions you need to provide a pointer to an mxArray of class mxUINT8_CLASS.

3) Exactly what I'm saying, they're not the same, compare the plhs input argument below.

I took your mexFunction line and the corresponding lines from mex.h. (BTW, I only found it because my compiler gave an error, I would strongly recommend to try to compile your MEX source code to find errors/warnings for yourself, if your code doesn't compile there's definitely something wrong)

void mexFunction(int nlhs,mxArray *plhs,int nrhs,const mxArray *prhs[])

void mexFunction(
    int           nlhs,           /* number of expected outputs */
    mxArray       *plhs[],        /* array of pointers to output arguments */
    int           nrhs,           /* number of inputs */
    const mxArray *prhs[]         /* array of pointers to input arguments */
    );

[tensorisation]

I've been away for a few days.

(1) in regards to 2.) : I looked in the example, so you can only use logical arrays for the output?

I should just turn my logical array to an integer class in matlab beforehand then for it to work properly, right? But for other arrays that are just default in Matlab I don't have to convert into a specific class beforehand? (for example: a=[1,5,7] in Matlab, needs to be converted to an integer beforehand?)

(2) In regards to 3.) : Oh, you mean it's *plhs instead of *plhs[]? I fixed it now.

(3) I went by the examples but I'm still not entirely sure why it uses mwSize and size_t instead of just int, can you explain this?

I now have:

#include "mex.h"


mwSize mod(mwSize x,mwSize N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}


double d_abs(double x) /* a function for the absolute value operator */
{
if(x<0)
    {
    return -x;
    }
else
    {
    return x;
    }
}

                    
double max(double x,double y) /* a function for the max operator */
{
if(x>y)
    {
    return x;
    }
else
    {
    return y;
    }
}

/* in matlab i will invoke potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max); */

/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */          
void potentials_evolution(mwSize L,mwSize max_relax_iters,mwSize N_eff_occ_sites,mwSize *i_eff_V_bounded,mwSize *j_eff_V_bounded,mwSize *k_eff_V_bounded,mwSize ***system,double over_relax_fact,double ***V,double *delta_V_rms,double *delta_V_abs,double *delta_V_max)
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */ 
        }
    }
}


/* in matlab: [V,delta_V_rms,delta_V_abs,delta_V_max] = potentials_evolution(i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V_0,delta_V_rms_0,delta_V_abs_0,delta_V_max_0) */
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
size_t L, max_relax_iters, N_eff_occ_sites;

size_t *i_eff_V_bounded, *j_eff_V_bounded, *k_eff_V_bounded;   
size_t ***system;       
double over_relax_fact;        
double ***V, ***V_0;        
double *delta_V_rms, *delta_V_abs, *delta_V_max, *delta_V_rms_0, *delta_V_abs_0, *delta_V_max_0;        


/* additional input for the calculation */
L = mxGetDimensions(prhs[5])[0];
max_relax_iters = mxGetN(prhs[6]);
N_eff_occ_sites = mxGetN(prhs[0]);

/* input arguments */
i_eff_V_bounded = mxGetInt32s(prhs[0]);
j_eff_V_bounded = mxGetInt32s(prhs[1]);
k_eff_V_bounded = mxGetInt32s(prhs[2]);
system = mxGetInt32s(prhs[3]);
over_relax_fact = mxGetScalar(prhs[4]);
V_0 = mxGetDoubles(prhs[5]);
delta_V_rms_0 = mxGetDoubles(prhs[6]);
delta_V_abs_0 = mxGetDoubles(prhs[7]);
delta_V_max_0 = mxGetDoubles(prhs[8]);        

/* output arguments */
plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[5]),mxGetDimensions(prhs[5]),mxDOUBLE_CLASS,mxREAL);
V = mxGetDoubles(plhs[0]);
plhs[1] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_rms = mxGetDoubles(plhs[1]);
plhs[2] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_abs = mxGetDoubles(plhs[2]);  
plhs[3] = mxCreateNumericMatrix(max_relax_iters,1,mxDOUBLE_CLASS,mxREAL);
delta_V_max = mxGetDoubles(plhs[3]);  
        
/* performing the calculation */
potentials_evolution((mwSize)L,(mwSize)max_relax_iters,(mwSize)N_eff_occ_sites,(mwSize)i_eff_V_bounded,(mwSize)j_eff_V_bounded,(mwSize)k_eff_V_bounded,(mwSize)system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);
}

[tensorisation]

Sadly, I'm still trying to make this work.