MATLAB: How to apply ‘parfor’ in this for-loop code

parfor

      mu_11=zeros(Z_N,N_c-1,Z_N,N_c-1);
      nu_11=zeros(Z_N,N_c-1,Z_N,N_c-1);
      syms s_1 s_2
      for n=1:Z_N
          for k=1:N_c-1
              for m=0:Z_N-1 % Defining Functions
                  f=cot(pi*(z_s1_cavityside(n,k)-l(m+1)*s_1-zeta(m+1))*exp(-i*(pi/2-beta*pi/180))/Z_N/h-m*pi/Z_N);
                  g=s_1*cot(pi*(z_s1_cavityside(n,k)-l(m+1)*s_1-zeta(m+1))*exp(-i*(pi/2-beta*pi/180))/Z_N/h-m*pi/Z_N);
                  for p=1:N_c-1 % Perform Numerical Integration
                      mu_11(n,k,(m+1),p)=vpaintegral(f,s_1,S_1(p),S_1(p+1),'RelTol',1e-4,'AbsTol',0,'MaxFunctionCalls',Inf);
                      nu_11(n,k,(m+1),p)=vpaintegral(g,s_1,S_1(p),S_1(p+1),'RelTol',1e-4,'AbsTol',0,'MaxFunctionCalls',Inf);
                  end
              end
          end
          progress_11=100*n/Z_N
      end
Why can't I apply the 'parfor' at the first 'for-loop'? I thought it would be naturally possible because the matrix index is all independent of the variables…
The error message says that in the lowest for-loop, the usage of mu_11 and nu_11 4D matrices is inappropriate for applying the parfor statement

Best Answer

You need to rewrite the code in approximately the form:
      mu_11=zeros(Z_N,N_c-1,Z_N,N_c-1);
      nu_11=zeros(Z_N,N_c-1,Z_N,N_c-1);
      syms s_1 s_2
      for n=1:Z_N
          mu_11n = zeros(N_c-1,Z_N,N_c-1);
          nu_11n = zeros(N_c-1,Z_N,N_c-1);
          z_s1_cavitysiden = z_s1_cavityside(n,:);
          for k=1:N_c-1
              for m=0:Z_N-1 % Defining Functions
                  f=cot(pi*(z_s1_cavitysiden(k)-l(m+1)*s_1-zeta(m+1))*exp(-i*(pi/2-beta*pi/180))/Z_N/h-m*pi/Z_N);
                  g=s_1*cot(pi*(z_s1_cavitysiden(k)-l(m+1)*s_1-zeta(m+1))*exp(-i*(pi/2-beta*pi/180))/Z_N/h-m*pi/Z_N);
                  for p=1:N_c-1 % Perform Numerical Integration
                      mu_11n(k,(m+1),p)=vpaintegral(f,s_1,S_1(p),S_1(p+1),'RelTol',1e-4,'AbsTol',0,'MaxFunctionCalls',Inf);
                      nu_11n(k,(m+1),p)=vpaintegral(g,s_1,S_1(p),S_1(p+1),'RelTol',1e-4,'AbsTol',0,'MaxFunctionCalls',Inf);
                  end
              end
          end
          mu_11(n,:,:,:) = mu_11n;
          nu_11(n,:,:,:) = nu_11n;
          progress_11=100*n/Z_N
      end
It is also recommended to move calculations out of loops when possible. For example, -i*(pi/2-beta*pi/180))/Z_N/h does not change with n, k, or m, and could be pre-computed. Likewise, z_s1_cavitysiden(k) does not depend upon m, so the value could be computed inside the k loop instead of inside the m loop.
Related Question