我写了一个涉及进行数值计算（松弛）的 Matlab 代码，但它很慢。我了解到使用 mex 文件运行 C 代码并将其集成到 Matlab 中的可能性，所以我正在考虑在 C 中进行数值计算（这相对简单但涉及循环并且需要时间），其余的（之前及之后）在 Matlab 中。
我的 Matlab 代码中完成计算的部分：
% 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
所以在 C 中它应该是这样的：
#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 */
}
}
}
因此，在 Matlab 中，我将上面显示的部分 Matlab 代码替换为：
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);
我该如何实施？我试图寻找一种简单的方法来做到这一点，但我无法弄清楚如何正确地做到这一点。
注意 1：对于随机生成的不同系统（存在一个for循环遍历不同系统），此数值计算不仅执行一次，而且执行多次。
注 2：我的 C 很生锈。