How to adapt my for loop to get it working with parfor? (ODE solver)

Hi guys,

I finally got as far to changing my for loop code in the ODE solver into parfor code. However, something is off which keeps me from getting my old for-loop results. Im really stuck. Could someone please be so kind to look at it? I dont know what else to do besides put the skimmed down code here. Please ask for more/less details!

What youre looking at is there are 66+48+48 ODEs per SACn, for which I reserve blocks of 200 ODEs. My thought was the following: per time step each SACn block of 200 ODEs is uncoupled. Only in the NEXT time step they are coupled at the line gGABA /gtemp = ...Codraai(...). 'Codraai' is a lookup matrix to check to which other ODE of another SACn block each ODE is coupled.

For version:

function [dy t] = encisoODE(t,y,Sgt,Syngg,SACnum,Codraai)
ODEs=SACnum*200;
dy = zeros(ODEs,size(y,2));                                                          
gGABA=zeros(ODEs,size(y,2));

...........some constants..........

gGABA=gGABA+gGABA_rest;

Syg = interp1q(Sgt,Syngg,t);                             % INTERP1 = interpolate for noninteger values of t
Sg = repmat(Syg,size(y,2),1);                            % template voor alle mogelijkheden y voor dezelfde t

for SACn=200*(0:SACnum-1)
    i=19+SACn:1:66+SACn;
    dy(i-18+100,:) = alpha.*(1-y(i-18+100,:)).*(1./(1+exp(-((y(i,:)  -th1)./k1))))-beta.*y(i-18+100,:);
    dy(i-18+150,:) = alpha.*(1-y(i-18+150,:)).*(1./(1+exp(-((y(i-18+100,:)-th2)./k2))))-beta.*y(i-18+150,:);

      i=1+SACn:1:66+SACn;  
      gGABA(i,:)=gGABA_rest+(gGABA_bound-gGABA_rest).*...
                 (y(Codraai(i,1)-18+150,:)+y(Codraai(i,2)-18+150,:)+y(Codraai(i,3)-18+150,:)+y(Codraai(i,4)-18+150,:)+...
                  y(Codraai(i,5)-18+150,:)+y(Codraai(i,6)-18+150,:)+y(Codraai(i,7)-18+150,:)+y(Codraai(i,8)-18+150,:));       

      i=1+SACn:3:18+SACn;   
      dy(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i)'.*(Egluta-y(i,:))+gGABA(i,:).*(EGABAp-y(i,:))+gK.*(EK-y(i,:))+delta.*((y(i+1,:)-y(i,:))));
      i=2+SACn:3:18+SACn;  
      dy(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i)'.*(Egluta-y(i,:))+gGABA(i,:).*(EGABAp-y(i,:))+gK.*(EK-y(i,:))+delta.*((y(i-1,:)-y(i,:))+(y(i+1,:)-y(i,:))));
      i=3+SACn:3:18+SACn;  
      dy(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i)'.*(Egluta-y(i,:))+gGABA(i,:).*(EGABAp-y(i,:))+gK.*(EK-y(i,:))+delta.*((y(i-1,:)-y(i,:))+(y(i+17,:)-y(i,:))));
      i=20+SACn:3:36+SACn;  
      dy(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i)'.*(Egluta-y(i,:))+gGABA(i,:).*(EGABAd-y(i,:))+gK.*(EK-y(i,:))+delta.*((y(i-17,:)-y(i,:))+(y(i+19+(i-(SACn+20))/3*2,:)-y(i,:))));
      i=39+SACn:5:66+SACn; 
      dy(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i)'.*(Egluta-y(i,:))+gGABA(i,:).*(EGABAd-y(i,:))+gK.*(EK-y(i,:))+delta.*((y(i-19-(i-(SACn+39))/5*2,:)-y(i,:))));
  end

end

parfor version:

function [dy t ggabavt] = encisoODE(t,y,Sgt,Syngg,SACnum,L,Co,Codraai)
ODEs=SACnum*200;
dy = zeros(ODEs,size(y,2));                                                                                                                              
gGABA=zeros(ODEs,size(y,2));

..........some constants...............

    gGABA=gGABA+gGABA_rest;

Syg = interp1q(Sgt,Syngg,t);                             % INTERP1 = interpolate for noninteger values of t
Sg = repmat(Syg,size(y,2),1);                            % template voor alle mogelijkheden y voor dezelfde t

%parfor temporaries
dsgt1=zeros(SACnum,48,size(y,2));
dsgt2=zeros(SACnum,48,size(y,2));
sgt1=zeros(SACnum,48,size(y,2));
sgt2=zeros(SACnum,48,size(y,2));
gt=zeros(SACnum,66,size(y,2));
dyt3=zeros(SACnum,66,size(y,2));
yt3=zeros(SACnum,66,size(y,2));

parfor SACn=0:SACnum-1  
    dsgtemp1=zeros(48,size(y,2));
    dsgtemp2=zeros(48,size(y,2));
    sgtemp1=zeros(48,size(y,2));
    sgtemp2=zeros(48,size(y,2));
    gtemp=zeros(66,size(y,2));
    dytemp3=zeros(66,size(y,2));
    ytemp3=zeros(66,size(y,2));

      i=19:1:66;              
      dsgtemp1(i-18,:) = alpha.*(1-sgtemp1(i-18,:)).*(1./(1+exp(-((ytemp3(i,:)  -th1)./k1))))-beta.*sgtemp1(i-18,:);      
      dsgtemp2(i-18,:) = alpha.*(1-sgtemp2(i-18,:)).*(1./(1+exp(-((sgtemp1(i-18,:)-th2)./k2))))-beta.*sgtemp2(i-18,:); 

      dsgt1(SACn+1,:,:)=dsgtemp1;
      dsgt2(SACn+1,:,:)=dsgtemp2;
      sgt1(SACn+1,:,:)=sgtemp1;
      sgt2(SACn+1,:,:)=sgtemp2;

      i=1:1:66; 
      gtemp(i,:)=gGABA_rest+(gGABA_bound-gGABA_rest).*...
                 (y(Codraai(i+200*SACn,1)-18+150,:)+y(Codraai(i+200*SACn,2)-18+150,:)+y(Codraai(i+200*SACn,3)-18+150,:)+y(Codraai(i+200*SACn,4)-18+150,:)+...
                  y(Codraai(i+200*SACn,5)-18+150,:)+y(Codraai(i+200*SACn,6)-18+150,:)+y(Codraai(i+200*SACn,7)-18+150,:)+y(Codraai(i+200*SACn,8)-18+150,:));

      gt(SACn+1,:,:)=gtemp;

      i=1:3:18;  
      dytemp3(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i+200*SACn)'.*(Egluta-ytemp3(i,:))+gtemp(i,:).*(EGABAp-ytemp3(i,:))+gK.*(EK-ytemp3(i,:))+delta.*((ytemp3(i+1,:)-ytemp3(i,:))));
      i=2:3:18;  
      dytemp3(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i+200*SACn)'.*(Egluta-ytemp3(i,:))+gtemp(i,:).*(EGABAp-ytemp3(i,:))+gK.*(EK-ytemp3(i,:))+delta.*((ytemp3(i-1,:)-ytemp3(i,:))+(ytemp3(i+1,:)-ytemp3(i,:))));
      i=3:3:18;  
      dytemp3(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i+200*SACn)'.*(Egluta-ytemp3(i,:))+gtemp(i,:).*(EGABAp-ytemp3(i,:))+gK.*(EK-ytemp3(i,:))+delta.*((ytemp3(i-1,:)-ytemp3(i,:))+(ytemp3(i+17,:)-ytemp3(i,:))));
      i=20:3:36;  
      dytemp3(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i+200*SACn)'.*(Egluta-ytemp3(i,:))+gtemp(i,:).*(EGABAd-ytemp3(i,:))+gK.*(EK-ytemp3(i,:))+delta.*((ytemp3(i-17,:)-ytemp3(i,:))+(ytemp3(i+19+(i-(+20))/3*2,:)-ytemp3(i,:))));
      i=39:5:66; 
      dytemp3(i,:) = 1./tau.*(ggluta_bound.*Sg(:,i+200*SACn)'.*(Egluta-ytemp3(i,:))+gtemp(i,:).*(EGABAd-ytemp3(i,:))+gK.*(EK-ytemp3(i,:))+delta.*((ytemp3(i-19-(i-(+39))/5*2,:)-ytemp3(i,:))));

      dyt3(SACn+1,:,:)=dytemp3;
      yt3(SACn+1,:,:)=ytemp3;
  end

%unpicking all the temps to fill dy
for SACn=0:SACnum-1
dy(SACn*200+101:SACn*200+148,:)=dsgt1(SACn+1,:,:);
dy(SACn*200+151:SACn*200+198,:)=dsgt2(SACn+1,:,:);
y(SACn*200+101:SACn*200+148,:)=sgt1(SACn+1,:,:);
y(SACn*200+151:SACn*200+198,:)=sgt2(SACn+1,:,:);

gGABA(SACn*200+1:SACn*200+66,:)=gt(SACn+1,:,:);

dy(SACn*200+1:SACn*200+66,:)=dyt3(SACn+1,:,:);
y(SACn*200+1:SACn*200+66,:)=yt3(SACn+1,:,:);
end

end

[Sanne]

Ok apparantly this is very difficult to do since the parallelization must concentrate on one time step.. does anybody know if this is possible?