problem with unclassified variable in parfor loop

I have a parfor loop in a section of code that is in error as the variable msdx is unclassfied. I understand this means Matlab cannot specifically define the variable. Does anyone know of a workaround? xu(i,j) is defined and works before it gets to this parfor loop. If I don't have the parfor loop (substitute in a for loop), it works fine.

Here is the code:

for i=1:ntimes %number value sets after inytial
    parfor j=1:nH
        % xzdiff(i,j)=xu(i,j)-zu(i,j);
        xuo=xu(1,j);
        yuo=yu(1,j);
        zuo=zu(1,j);
        msdx(i,j)=(xu(i,j)-xuo)^2;
        msdy(i,j)=(yu(i,j)-yuo)^2;
        msdz(i,j)=(zu(i,j)-zuo)^2;
    end
end

[Edric Ellis]

I can't see a problem here - the form of indexing for msdx looks OK for parfor. Which release of MATLAB are you using? Can you post a minimal reproduction that is (or should be) executable?

[richg]

% Here is some simplified code that executes and reproudces the error. Sorry abou the extra lines-inserted by forum edior. I dont find how to correct this.

ntimes=10;
nh=10;
ntemp=2;
tp=zeros(ntemp);
time_dt=0.001; %this time interval in picoseconds from LAMMPS
for i=0:ntimes
     tp=i*time_dt; % should be 0,5000,10000 .. 95000
end

for ik=1:2
    
msdx=zeros(ntimes,nh);
msdy=zeros(ntimes,nh);
msdz=zeros(ntimes,nh);
msdsx=zeros(ntimes);
msdsy=zeros(ntimes);
msdsz=zeros(ntimes);
   
  switch(ik) % intialize data set 1 and datset 2 for test
      case 1
xu1= [1.344 1.794 2.027 1.764 1.583 1.965 1.836 1.858 1.848 2.109];
yu1 = [66.222 37.101 60.380 48.642 21.410 29.483 5.567 2.034 19.569 21.330];
zu1 = [8.087  9.5287 40.727 1.599  68.364 56.351 9.78  1.899 46.770 36.80];
xu2 = [1.794  1.986  1.712  1.478  1.314  2.022  2.253 2.302 1.825 1.847];
yu2 = [66.404 37.063 60.583 48.333 21.54  29.446 5.785 1.861 19.182 21.365];
zu2 = [7.349  9.521  40.82  2.192  68.43  56.460 9.919 1.424 46.80  37.048];
xu3=[2.192  2.044  1.955  1.975  2.165  1.617  1.979 2.085 1.985  1.980]; 
yu3=[65.822 36.843 60.282 48.532 21.673 29.551 5.768 1.553 19.732 21.339];
zu3=[8.077 9.501 40.708 2.205 68.451 56.251 9.368 2.174 46.49 37.135];
xu4=[2.092  1.853  1.595  2.065  1.915  1.855  1.769 1.599 2.202 2.135];
yu4=[66.689 37.187 60.394 48.932 21.089 29.356 5.565 1.808 19.25 21.630]; 
zu4=[7.544  9.655  40.930 1.831  68.155 56.654 9.737 1.872 46.47 36.18];
xu5=[1.911 1.914 1.783 1.494 1.771 1.348 2.186 1.653 2.136 1.882];
yu5=[66.162 36.799 60.356 48.547 21.494 29.258 5.617 1.685 19.556 21.402];
zu5=[7.889 9.420 41.042 1.440 68.170 56.539 9.381 2.011 46.684 36.820];
xu6=[1.582 1.759 1.409 1.840 1.914 1.977 2.315 1.797 2.079 1.842];
yu6=[66.294 37.189 60.424 48.921 21.656 29.191 5.384 1.637 19.291 21.376];
zu6=[7.890 9.452 40.833 1.788 68.256 56.597 9.616 1.730 46.90 37.182];
xu7=[1.793 1.778 1.898 1.724 2.167 1.504 1.396 1.923 2.136 1.721];
yu7=[66.378 36.750 60.517 48.925 21.223 29.80 5.442 1.81 19.071 21.230];
zu7=[7.704 9.344 40.928 1.899 68.459 56.386 9.569 1.478 46.568 37.055];
xu8=[1.626 1.782 1.964 1.851 1.922 1.498 2.260 1.589 3.325 2.289];
yu8=[66.422 36.560 60.424 48.748 21.155 29.169 5.620 1.634 19.904 21.711];
zu8=[7.778 9.310 40.652 1.892 69.405 56.745 9.696 2.066 45.049 37.014];
xu9=[1.636 1.792 1.974 1.861 1.932 1.508 2.270 1.599 3.335 2.299];
yu9=[66.432 36.570 60.434 48.758 21.165 29.179 5.630 1.644 19.914 21.721];
zu9=[7.88 9.320 40.662 1.902 69.415 5.6755 9.706 2.076 45.059 37.029]; 
xu10=[1.856 1.810 1.980 1.624 2.082 -0.442 2.039 1.880 3.783 1.804]; 
yu10=[66.340 36.950 60.720 48.584 21.443 29.197 6.085 1.610 19.646 21.418];
zu10=[7.804 9.664 40.93 2.757 60.291 54.544 9.696 2.799 44.495 37.426];
xu=[xu1;xu2;xu3;xu4;xu5;xu6;xu7;xu8;xu9;xu10]
yu=[yu1;yu2;yu3;yu4;yu5;yu6;yu7;yu9;yu9;yu10]
zu=[zu1;zu2;zu3;zu4;zu5;zu6;zu7;zu8;zu9;zu10]
        otherwise
xu1=[-1.848 1.822 -6.367 0.050 -4.131 -0.041 1.583 9,798 -2.056 11.801];
yu1=[66.319 31.404 60.609 52.433 17.033 39.316 9.783 2.154 1.889 19.268];
zu1=[19.434 0.0208 37.068 0.068 60.793 5.661 17.601 5.488 52.559 36.745];
xu2= [-1.990 1.72113 -5.882 7.828 -4.491 -0.012 1.385 9.566 0.0174 11.760];
yu2=[65.859 31.090 60.218 56.629 17.331 38.516 9.503 1.961 3.905 19.306];
zu2=[19.111 0.068 37.102 -0.355 69.818 52.321 17.314 5.642 52.8-7 37.314];
xu3=[-2.229 1.490 -6.213 7.823 -3.749-2.986 1.800 9.683 -.430 13.591];
yu3=[66.126 31.115 60.265 52.399 13.546 35.814 9.722 1.633 3.909 19.226]; 
zu3=[19.301 -0.214 36.958 -0.081 66.214 53.645 17.292 5.719 52.124 35.105];
xu4=[-2.454 1.583 -5.942 7.672 -5.775 -2.105 1.470 8.110 -0.0454 11.214];
yu4=[64.070 32.435 60.623 52.671 13.383 36.781 9.742 2.104 3.452 19.659];
zu4=[17.131 -0.150 36.988 -.277 64.479 52.136 17.626 3.769 52.572 33.295];
xu5=[-1.766 1.590 -5.885 7.801 -5.68 -2.114 1.688 7.820 -9.489 11.704];
yu5=[64.859 31.316 60.270 52.77 13.246 35.453 9.829 1.738 4.073 19.844];
zu5=[17.191 -0.181 36.955 -0.313 64.516 50.793 17.599 3.663 52.491 33.122];
xu6=[-2.529 1.5793 -6.058 7.739 -6.568 -2.151 1.965 7.865 0.0783 11.554];  
yu6=[64.209 31.243 60.1323 52.373 13.175 35.134 9.861 1.861 3.391 19.587];
zu6=[17.755 -0.244 36.829 -.3012 64.004 50.475 17.374 3.723 52.378 32.328];
xu7=[-1.840 3.611 -6.052 7.287 -6.060 -2.213 1.886 7.399 -0.0987 11.864];
yu7=[64.304 30.883 60.216 52.6634 13.208 35.078 9.442 2.012 3.616 19.118];
zu7=[17.539 2.216 36.918 -0.057 64.368 50.725 17.461 3.878 52.130 32.899];
xu8=[-2.451 4.012 -6.330 7.957 -5.672 -1.918 1.406 7.592 -0.185 13.922];  
yu8=[64.172 30.922 60.446 52.840 13.258 35.093 9.745 1.666 3.799 19.360];
zu8=[17.334 1.462 36.898 -0.208 63.698 50.439 17.460 3.831 52.295 34.784];
xu9=[-2.146 4.043 -5.246 7.286 -6.278 -2.016 1.533 7.507 -0.4214 14.044];
yu9=[64.373 31.066 61.065 52.541 13.339 35.114 9.312 1.347 3.726 19.737];
zu9=[17.347 1.648 36.416 -0.291 64.564 50.745 17.533 3.994 52.660 35.142];
xu10=[-2.080 3.972 -3.89 7.828 -6.381 -2.184 1.591 7.782 -0.118 13.591];
yu10=[64.287 31.249 64.800 52.597 13.434 35.303 9.325 1.791 4.018 19.059];
zu10=[17.429 1.645 31.252 0.093 64.380 50.712 17.366 3.620 52.361 34.629];
xu=[xu1;xu2;xu3;xu4;xu5;xu6;xu7;xu8;xu9;xu10]
yu=[yu1;yu2;yu3;yu4;yu5;yu6;yu7;yu9;yu9;yu20]
zu=[zu1;zu2;zu3;zu4;zu5;zu6;zu7;zu8;zu9;zu10]
    end
  
for i=1:ntimes %number value sets after inytial
        parfor j=1:nh 
           % xzdiff(i,j)=xu(i,j)-zu(i,j);
            xuo=xu(1,j);
            yuo=yu(1,j);
            zuo=zu(1,j);
            msdx(i,j)=(xu(i,j)-xuo)^2;
            msdy(i,j)=(yu(i,j)-yuo)^2;
            msdz(i,j)=(zu(i,j)-zuo)^2;
               
        end
end

for ii=1:ntimes
      
        parfor j=1,nh ;%sum msd fr all atoms
            msdsx(ii)=msdsx(ii)+ msdx(ii,j);% sum for all atoms at each time i
            msdsy(ii)=msdsy(ii)+ msdy(ii,j);
            msdsz(ii)=msdsz(ii)+ msdz(ii,j);
            msds(ii)= msdsx(ii) + msdsy(ii) + msdsz(ii);
            
        end
clc       
    end
        [px,Sx,mux]=polyfit(msdsx,tp,1);
        [py,Sy,muy]=polyfit(msdsy,tp,1);
        [pz,Sz,muz]=polyfit(msdsz,tp,1);
        [ps,Ss,mus]=polyfit(msds,tp,1);
  
                  
end

        displ('finished overall')

Answer 1

[Walter Roberson]

xu1 = xu(1,:);
yu1 = yu(1,:);
zu1 = zu(1,:);
for i=1:ntimes %number value sets after inytial
    xui = xu(i,:);
    yui = yu(i,:);
    zui = zu(i,:);
    parfor j=1:nH
        % xzdiff(i,j)=xu(i,j)-zu(i,j);
        xuo = xu1(j);
        yuo = yu1(j);
        zuo = zu1(j);
        msdx(i,j) = (xui(j)-xuo)^2;
        msdy(i,j) = (yui(j)-yuo)^2;
        msdz(i,j) = (zui(j)-zuo)^2;
    end
end

The difficulty is that

xuo=xu(1,j);

and

msdx(i,j)=(xu(i,j)-xuo)^2;

index the same variable, xu, using a different combination of subscripts within the parfor, which is not permitted. The work-around I use here is to optimize out the xu(i,:) into an xui variable that is set outside the parfor according to the for i iteration.

You should be to vectorize the calculation without using parfor, something like

for i=1:ntimes %number value sets after inytial
  msdx(i,:) = (xu(i,:) - xu(1,:)).^2;
  msdy(i,:) = (yu(i,:) - yu(1,:)).^2;
  msdz(i,:) = (zu(i,:) = zu(1,:)).^2;
end

and in turn you could get rid of the loop,

msdx = (xu - xu(1,:)).^2;
msdy = (yu - yu(1,:)).^2;
msdz = (zu - zu(1,:)).^2;

and then instead of using explicit parfor, let MATLAB automatically delegate to MKL or LAPACK as appropriate.

[richg]

I have added sample executing code that reproduces the problem, but simplified. Can you tell if your response applies? I ran this on my laptop on R202a version.

How do you get Matlab ti delegate to MKL? According to information on Web, LAPACK linked into Matlab in R2020. The latest on the supercomputer at my university running the full program is R2018B. Does this mean I should ask them to install R2020a or R2020b to accomplish this? Is there a link to Matlab docucmention discussing this?

Answer 2

[richg]

I am running it on a supercomputer, version is R2018a. I also have R2020a on my laptop. I will create an executable example that will illustrate using my laptop and post it.

[richg]

I created the exicuable and posted it above.

[Walter Roberson]

Note that a number of fixes needed to be made, including in the constants. For example you had a place with a comma instead of a decimal place, and another place where you did not have a space between two adjacent elements.

I did not fully vectorize the code, as I thought it was important that you saw the techniques needed to get the parfor to work. You would very likely be better off vectorizing instead of using parfor.

ntimes=10;
nh=10;
tp=zeros(ntimes,1);
time_dt=0.001; %this time interval in picoseconds from LAMMPS
for i=0:ntimes-1
    tp(i+1)=i*time_dt; % should be 0,5000,10000 .. 95000
end

for ik=1:2
    
    msdx=zeros(ntimes,nh);
    msdy=zeros(ntimes,nh);
    msdz=zeros(ntimes,nh);
    msdsx=zeros(ntimes,1);
    msdsy=zeros(ntimes,1);
    msdsz=zeros(ntimes,1);
    
    switch(ik) % intialize data set 1 and datset 2 for test
        case 1
            xu1= [1.344 1.794 2.027 1.764 1.583 1.965 1.836 1.858 1.848 2.109];
            yu1 = [66.222 37.101 60.380 48.642 21.410 29.483 5.567 2.034 19.569 21.330];
            zu1 = [8.087  9.5287 40.727 1.599  68.364 56.351 9.78  1.899 46.770 36.80];
            xu2 = [1.794  1.986  1.712  1.478  1.314  2.022  2.253 2.302 1.825 1.847];
            yu2 = [66.404 37.063 60.583 48.333 21.54  29.446 5.785 1.861 19.182 21.365];
            zu2 = [7.349  9.521  40.82  2.192  68.43  56.460 9.919 1.424 46.80  37.048];
            xu3=[2.192  2.044  1.955  1.975  2.165  1.617  1.979 2.085 1.985  1.980];
            yu3=[65.822 36.843 60.282 48.532 21.673 29.551 5.768 1.553 19.732 21.339];
            zu3=[8.077 9.501 40.708 2.205 68.451 56.251 9.368 2.174 46.49 37.135];
            xu4=[2.092  1.853  1.595  2.065  1.915  1.855  1.769 1.599 2.202 2.135];
            yu4=[66.689 37.187 60.394 48.932 21.089 29.356 5.565 1.808 19.25 21.630];
            zu4=[7.544  9.655  40.930 1.831  68.155 56.654 9.737 1.872 46.47 36.18];
            xu5=[1.911 1.914 1.783 1.494 1.771 1.348 2.186 1.653 2.136 1.882];
            yu5=[66.162 36.799 60.356 48.547 21.494 29.258 5.617 1.685 19.556 21.402];
            zu5=[7.889 9.420 41.042 1.440 68.170 56.539 9.381 2.011 46.684 36.820];
            xu6=[1.582 1.759 1.409 1.840 1.914 1.977 2.315 1.797 2.079 1.842];
            yu6=[66.294 37.189 60.424 48.921 21.656 29.191 5.384 1.637 19.291 21.376];
            zu6=[7.890 9.452 40.833 1.788 68.256 56.597 9.616 1.730 46.90 37.182];
            xu7=[1.793 1.778 1.898 1.724 2.167 1.504 1.396 1.923 2.136 1.721];
            yu7=[66.378 36.750 60.517 48.925 21.223 29.80 5.442 1.81 19.071 21.230];
            zu7=[7.704 9.344 40.928 1.899 68.459 56.386 9.569 1.478 46.568 37.055];
            xu8=[1.626 1.782 1.964 1.851 1.922 1.498 2.260 1.589 3.325 2.289];
            yu8=[66.422 36.560 60.424 48.748 21.155 29.169 5.620 1.634 19.904 21.711];
            zu8=[7.778 9.310 40.652 1.892 69.405 56.745 9.696 2.066 45.049 37.014];
            xu9=[1.636 1.792 1.974 1.861 1.932 1.508 2.270 1.599 3.335 2.299];
            yu9=[66.432 36.570 60.434 48.758 21.165 29.179 5.630 1.644 19.914 21.721];
            zu9=[7.88 9.320 40.662 1.902 69.415 5.6755 9.706 2.076 45.059 37.029];
            xu10=[1.856 1.810 1.980 1.624 2.082 -0.442 2.039 1.880 3.783 1.804];
            yu10=[66.340 36.950 60.720 48.584 21.443 29.197 6.085 1.610 19.646 21.418];
            zu10=[7.804 9.664 40.93 2.757 60.291 54.544 9.696 2.799 44.495 37.426];
            xu=[xu1;xu2;xu3;xu4;xu5;xu6;xu7;xu8;xu9;xu10];
            yu=[yu1;yu2;yu3;yu4;yu5;yu6;yu7;yu9;yu9;yu10];
            zu=[zu1;zu2;zu3;zu4;zu5;zu6;zu7;zu8;zu9;zu10];
        otherwise
            xu1=[-1.848 1.822 -6.367 0.050 -4.131 -0.041 1.583 9.798 -2.056 11.801];
            yu1=[66.319 31.404 60.609 52.433 17.033 39.316 9.783 2.154 1.889 19.268];
            zu1=[19.434 0.0208 37.068 0.068 60.793 5.661 17.601 5.488 52.559 36.745];
            xu2= [-1.990 1.72113 -5.882 7.828 -4.491 -0.012 1.385 9.566 0.0174 11.760];
            yu2=[65.859 31.090 60.218 56.629 17.331 38.516 9.503 1.961 3.905 19.306];
            zu2=[19.111 0.068 37.102 -0.355 69.818 52.321 17.314 5.642 52.8-7 37.314];
            xu3=[-2.229 1.490 -6.213 7.823 -3.749 -2.986 1.800 9.683 -.430 13.591];
            yu3=[66.126 31.115 60.265 52.399 13.546 35.814 9.722 1.633 3.909 19.226];
            zu3=[19.301 -0.214 36.958 -0.081 66.214 53.645 17.292 5.719 52.124 35.105];
            xu4=[-2.454 1.583 -5.942 7.672 -5.775 -2.105 1.470 8.110 -0.0454 11.214];
            yu4=[64.070 32.435 60.623 52.671 13.383 36.781 9.742 2.104 3.452 19.659];
            zu4=[17.131 -0.150 36.988 -.277 64.479 52.136 17.626 3.769 52.572 33.295];
            xu5=[-1.766 1.590 -5.885 7.801 -5.68 -2.114 1.688 7.820 -9.489 11.704];
            yu5=[64.859 31.316 60.270 52.77 13.246 35.453 9.829 1.738 4.073 19.844];
            zu5=[17.191 -0.181 36.955 -0.313 64.516 50.793 17.599 3.663 52.491 33.122];
            xu6=[-2.529 1.5793 -6.058 7.739 -6.568 -2.151 1.965 7.865 0.0783 11.554];
            yu6=[64.209 31.243 60.1323 52.373 13.175 35.134 9.861 1.861 3.391 19.587];
            zu6=[17.755 -0.244 36.829 -.3012 64.004 50.475 17.374 3.723 52.378 32.328];
            xu7=[-1.840 3.611 -6.052 7.287 -6.060 -2.213 1.886 7.399 -0.0987 11.864];
            yu7=[64.304 30.883 60.216 52.6634 13.208 35.078 9.442 2.012 3.616 19.118];
            zu7=[17.539 2.216 36.918 -0.057 64.368 50.725 17.461 3.878 52.130 32.899];
            xu8=[-2.451 4.012 -6.330 7.957 -5.672 -1.918 1.406 7.592 -0.185 13.922];
            yu8=[64.172 30.922 60.446 52.840 13.258 35.093 9.745 1.666 3.799 19.360];
            zu8=[17.334 1.462 36.898 -0.208 63.698 50.439 17.460 3.831 52.295 34.784];
            xu9=[-2.146 4.043 -5.246 7.286 -6.278 -2.016 1.533 7.507 -0.4214 14.044];
            yu9=[64.373 31.066 61.065 52.541 13.339 35.114 9.312 1.347 3.726 19.737];
            zu9=[17.347 1.648 36.416 -0.291 64.564 50.745 17.533 3.994 52.660 35.142];
            xu10=[-2.080 3.972 -3.89 7.828 -6.381 -2.184 1.591 7.782 -0.118 13.591];
            yu10=[64.287 31.249 64.800 52.597 13.434 35.303 9.325 1.791 4.018 19.059];
            zu10=[17.429 1.645 31.252 0.093 64.380 50.712 17.366 3.620 52.361 34.629];
            xu=[xu1;xu2;xu3;xu4;xu5;xu6;xu7;xu8;xu9;xu10];
            yu=[yu1;yu2;yu3;yu4;yu5;yu6;yu7;yu9;yu9;yu10];
            zu=[zu1;zu2;zu3;zu4;zu5;zu6;zu7;zu8;zu9;zu10];
    end
    
    for i=1:ntimes %number value sets after initial
        xu1 = xu(1,:);
        yu1 = yu(1,:);
        zu1 = zu(1,:);
        xui = xu(i,:);
        yui = yu(i,:);
        zui = zu(i,:);
        parfor j=1:nh
            xuo=xu1(j);
            yuo=yu1(j);
            zuo=zu1(j);
            msdx(i,j)=(xui(j)-xuo)^2;
            msdy(i,j)=(yui(j)-yuo)^2;
            msdz(i,j)=(zui(j)-zuo)^2;
            
        end
    end
    
    msds = zeros(ntimes,1);
    for ii=1:ntimes
        msdsxii = 0;
        msdsyii = 0;
        msdszii = 0;
        parfor j=1:nh %sum msd fr all atoms
            msdsxii = msdsxii + msdx(ii,j);% sum for all atoms at each time i
            msdsyii = msdsyii + msdy(ii,j);
            msdszii = msdszii + msdz(ii,j);
        end
        msdsx(ii) = msdsxii;
        msdsy(ii) = msdsyii;
        msdsz(ii) = msdszii;
        msds(ii) = msdsxii + msdsyii + msdszii;
        clc
    end
    [px,Sx,mux]=polyfit(msdsx,tp,1);
    [py,Sy,muy]=polyfit(msdsy,tp,1);
    [pz,Sz,muz]=polyfit(msdsz,tp,1);
    [ps,Ss,mus]=polyfit(msds,tp,1);
    
    
end

disp('finished overall')