我正在研究具有周期边界条件的厌世过程的扩散问题。在这个过程中，我考虑了一个具有M个位置和N个排斥相互作用粒子的一维离散晶格;每个位置可以在任何时间被任意数量的颗粒占据。
我已经有工作脚本来计算这个在我的个人工作站（4个并行工作者，没有GPU）的方式太慢，为我的目的。我可以访问以下共享硬件：10核英特尔酷睿i9 - 7900X处理器[20线程超线程]和128GB的RAM与Nvidia泰坦V GPU与12GB的GPU RAM。
我不是GPUMaven，但我想以某种方式利用可用的Nvidia Titan。我的代码是用Matlab编写的，共享的机器有Matlab 2018a。模拟记录的合奏的粒子和环的大小的合奏，以及在一个相互作用强度参数的时间演化。
这是嵌套的，如下所示：

loop over interaction strength parameter (for loop beta)
    loop over N                          (parallelize with parfor)
        loop over trajectories           (for loop n_traj)
            loop over time               (while loop t_f)

两个最大的循环是轨迹循环和时间循环。将我的任何变量（可能是n_traj）作为Matlab上的GPUarray是否会加快我的代码？
下面是脚本（用于快速测试：n_traj = 100，t_f = 30）

Repulsion=[6 5 2 1 0.5];
%--- lattice size
L=2; M=10; x=1:L*M;
%--- declares hopping rates  
k0=0.1; X=2; k_r=exp(X)*k0;   k_l=exp(-X)*k0; 
%--- non-interacting single particle drift and diffusion
D=(k_l+k_r)/2; v=(k_r-k_l);
%--- pre-allocates interacting particles drift and diffusion  
driftMEx=zeros(L*M,length(Repulsion)); 
diffuMEx=zeros(L*M,length(Repulsion));
%--- gets ensemble size and time duration
prompt=sprintf('Give the value No of trajectories:\n n_traj = ');
n_traj = input(prompt);
prompt=sprintf('\nGive the value time lenght:\n t_f = ');
t_f = input(prompt);

j=1;
for beta=Repulsion
    parfor i=1:L*M 
            [driftMEx(i,j),diffuMEx(i,j)]=simulation(i,M,X,k0,beta,n_traj,t_f);
    end 
    j=j+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                        Functions used in script                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%----       Montecarlo Simulations for Misanthrope Process      ----%%%%
function [drift,diffusion] = simulation(N,M,X,k_o,beta,n_traj,t_f)
%%%%%%%%%%%%%%%%%%%%%%% parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k_r=exp(X)*k_o; k_l=exp(-X)*k_o; 
s_min=1; % lowest state in boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%  Pre-allocation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_final=zeros(n_traj,1);
P=zeros(2*N+1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for jj=1:n_traj
    tc=0; % current time starts at 0
    state0=1:N;
    state0=mod(state0,M);
    state0(state0==0)=M;
    state=histcounts(state0,0.5+(0:M));
    bc_counter=0;
    hop=1;
    while tc<t_f  % loop over time
        % findMisanthrope finds index of where motors are this state
        N_pos=findMisanthrope(state,N);  
        T=zeros(N,M); % reset transition matrix
        % temporarily assign state of hopping fowards and backwards
        N_plus=N_pos+1;
        N_minus=N_pos-1;
        nn=histcounts(N_pos,0.5+(0:M));
        % finds possible moves
        for qq=1:N     
            % findMisanthrope finds index of where motors are this state
            N_pos=findMisanthrope(state,N);              
            if N_pos(qq)==1
                N_minus(qq)=M;
            end            
            if N_pos(qq)==M
                N_plus(qq)=1;
            end
            T(qq, N_minus(qq))=k_l; % rate of hopping left
            T(qq,N_plus(qq))=k_r;   % rate of hopping right
            
            alpha=0;
            bb=exp(-beta.*nn);
            aa=exp(alpha.*(nn-1));

            temp_qq=N_pos(qq);
            T(qq,setdiff(1:end,temp_qq))=T(qq,setdiff(1:end,temp_qq)).*bb(setdiff(1:end,temp_qq)).*aa(temp_qq);           
        end

        for mm=2:2:2*N
            P(mm)=T(mm/2,N_minus(mm/2))/sum(T(:)) + P(mm-1);
            P(mm+1)=T(mm/2,N_plus(mm/2))/sum(T(:)) + P(mm);
        end
        tau=(-log(rand))/(sum(T(:)));  % find time until next hop
        tc=tc+tau; % update current time
        R=rand; % pick another random number from 0 to 1 (r_2)
        iii=find(P<=R,1,'last'); % find which entry in vector correspons to r_2
        N_hopp=round(iii/2);  % which motor changes position
            state(N_pos(N_hopp))=state(N_pos(N_hopp))-1; % annihilate motor to hop
            N_plus=N_pos(N_hopp)+1;
            N_minus=N_pos(N_hopp)-1;
            
            if mod(iii/2,1)~=0   % if odd then motor hops left
                if N_pos(N_hopp)==s_min
                    N_minus=M;
                    bc_counter=bc_counter-1; % take note of how many times a motor crosses the boudary
                end
                state(N_minus)=1+state(N_minus); % create motor where it hops to
            else   % even then hops right
                if N_plus>M
                    N_plus=s_min;
                    bc_counter=bc_counter+1;
                end               
                state(N_plus)=1+state(N_plus);       
            end
        hop=hop+1; % advance hop index
    end
    N_f=findMisanthrope(state,N); % final position of motors
    delta=N_f-state0; % find change in state for drift calculation
    x_final(jj)=(sum(delta)+M*bc_counter); % add boundary condition info on. To find total distance travelled
end
% calculates diffusion
ave=mean(x_final);
diff_n(1:n_traj)=(x_final(1:n_traj)-ave).^2;
diffusion=sum(diff_n)/(2*n_traj*t_f);
% calculates drift
drift=sum(x_final(1:n_traj)/t_f)/n_traj;
end

%%%%----     Find indices of particles for Misanthrope Process     ----%%%%
function [rr]=findMisanthrope(tt,m)

    if 1*isempty(tt(tt>=2))==0    
        ttt=tt(tt~=0);
        rrr=find(tt);
        rr=zeros(1,m);
        rrt=cell(1,length(rrr));    
        for i=1:length(rrr)
            rrt{i}=rrr(i)*ones(1,ttt(i));
        end
        rr(1:m)=cell2mat(rrt);    
    else
        rr=find(tt);    
    end

end