# 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 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: My C is quite rusty.

• I think it would be easier if you stated your mathematical problem instead of writing the code... In this form it is not clear if your computation can be written in a vectorized way, thus improving the speed. – Beni Bogosel Aug 9 '19 at 11:20
• @BeniBogosel That would be a different question. Taking the calculation as is (as presented here), my aim is to run it in C (so it will be much faster) using a mex file. Looking around, I couldn't figure out how to properly implement it. – TensoR Aug 9 '19 at 11:45
• It would be a different question, but people would tell you their opinions on it. Off the top of my head, I wouldn't use a three dimensional array, but a one dimensional one for V. In this way, the computation of V1 could be implemented as a matrix-vector product, eliminating the need of the inner loop. – Beni Bogosel Aug 9 '19 at 19:59
• @BeniBogosel Turning V into a one dimensional array doesn't make much sense to me, but even if you can do it properly in a manner that makes sense (not sure how at the moment), it would still be much slower than running the code in C, no? – TensoR Aug 10 '19 at 15:13
• If you could simplify the inner loop using some sparse matrix *vector multiplications, those are quite efficient in Matlab. In some of my codes, using array multiplications instead of loops increased the speed 100 fold. Of course, it all depends if this is possible in your case or not. – Beni Bogosel Aug 11 '19 at 16:14