Parallel computation leaking memory and slowing a loop

I am doing a parallel computation, specifically a row by row integration of a symbolic matrix (sic I am integrating 4 rows at a time on my quadcore). 2 curious things occur:

*(1)* As the computation progresses, the free RAM steadily gets used up (see graph below). The initial surge in RAM is expected – it is the memory that is "actually" required by the computation. However, I am at a loss to explain why the remaining RAM then steadily continues to be used up. The issue is that eventually, after maybe an hour or so of computation and many rows having been integrated, the RAM is used up. MATLAB then starts to use virtual memory on the HDD, which decreases computation speeds to unacceptable levels (5 fold or so decrease in speed, with further slowing over time). I can circumvent this problem by dividing my matrix into smaller matrices (with fewer rows) that can individually be integrated before HDD switching is required. However, in my estimation this should not be necessary. I have plenty of RAM (16GB) and the actual computation only requires about 2900MB judging by the kink in the graph. Is MATLAB saving output to the RAM? What else could be messing up the memory allocation? The CPU effort never gets above ~65% and is typically around 40%.

<</matlabcentral/answers/uploaded_files/771/memory_allocation.png>>

*(2)* Much more concerning. Before the RAM is used up, a row integration typically takes about 60s per row. However, I just discovered that if I close the worker pool and effectively "reset" my processors (or instead simply restart matlab) and I manually instruct the same exact read/compute/write as is done in the loop, the same computation takes only fractions of a second (see command window dialogue below). I suspect this may be related to (1). This is completely unacceptable, why might this happening?

MATLAB command window output for integration of the first 4 rows:

Starting matlabpool using the ‘local’ profile … connected to 4 workers.
K_a_12
Elapsed time is 8.039642 seconds.
Commencing integration…
currently integrating row 1Elapsed time is 63.642951 seconds.
currently integrating row 3Elapsed time is 64.052436 seconds.
currently integrating row 4Elapsed time is 63.787564 seconds.
currently integrating row 2Elapsed time is 64.686729 seconds.

>> tic
K_a_12(2,:)=2*h_a*int(int(K_a_12_unintegrated(2,:),y,0,W),x);
toc
Elapsed time is 49.989623 seconds.

>> matlabpool close
Sending a stop signal to all the workers … stopped.
>> %Assign 4 cores
if matlabpool(‘size’) == 0
matlabpool open 4
end
Starting matlabpool using the ‘local’ profile … connected to 4 workers.
>> tic
K_a_12(2,:)=2*h_a*int(int(K_a_12_unintegrated(2,:),y,0,W),x);
toc
Elapsed time is 0.360260 seconds.
>>

My source code is attached below (note that symbolic toolbox is required). If you do not have a quadcore processor, you can change _matlabpool open 4_ to _matlabpool open_ [# of processors on your machine]:

%% Initialization

%Clear variables
clear;
clc;

%Assign 4 cores
if matlabpool(‘size’) == 0
matlabpool open 4
end

%% Problem description
syms G_a E_a h_a L1 L2 La W h_1 h_2 A1 B1 D1 G1 a11 a12 a16 a22 a26 a66 b11 b12 b16 b22 b26 b66 d11 d12 d16 d22 d26 d66 a44 a45 a55

% Laminate 1
A1=[a11,a12,a16;a12,a22,a26;a16,a26,a66];
B1=zeros(3);
D1=[d11,d12,d16;d12,d22,d26;d16,d26,d66];
G1=[a44,a45;a45,a55];

%Laminate 2
A2=A1;
B2=B1;
D2=D1;
G2=G1;

%% Assumed Displacement Field
m=10;
n=5;

syms x y z
B11 = sym(zeros(m+1, 1));
B21 = sym(zeros(m+1, 1));
B12 = sym(zeros(n+1, 1));
B22 = sym(zeros(n+1, 1));

for i = 0:m
A=factorial(m)/(factorial(i)*factorial(m-i))*((x-(L2-La))/L1)^i*(1-(x-(L2-La))/L1)^(m-i);
B11(i+1)=A;
end

for i = 0:n
B=factorial(n)/(factorial(i)*factorial(n-i))*(y/W)^i*(1-y/W)^(n-i);
B12(i+1)=B;
end

for i = 0:m
A=factorial(m)/(factorial(i)*factorial(m-i))*(x/L2)^i*(1-x/L2)^(m-i);
B21(i+1)=A;
end

for i = 0:n
B=factorial(n)/(factorial(i)*factorial(n-i))*(y/W)^i*(1-y/W)^(n-i);
B22(i+1)=B;
end

zero_vector=zeros((m+1)*(n+1),1);
counter=0;
for i = 0:m
for j=0:n
counter=counter+1;
C=B11(i+1)*B12(j+1);
Vx1(counter,1)=C; %x-displacement
Vy1(counter,1)=C; %y-displacement
Vz1(counter,1)=C; %z-displacement
Vrx1(counter,1)=C; %x-rotation
Vry1(counter,1)=C; %y-rotation
end
end
Vx1=[Vx1;zero_vector;zero_vector;zero_vector;zero_vector];
Vy1=[zero_vector;Vy1;zero_vector;zero_vector;zero_vector];
Vz1=[zero_vector;zero_vector;Vz1;zero_vector;zero_vector];
Vrx1=[zero_vector;zero_vector;zero_vector;Vrx1;zero_vector];
Vry1=[zero_vector;zero_vector;zero_vector;zero_vector;Vry1];

counter=0;
for i = 0:m
for j=0:n
counter=counter+1;
C=B21(i+1)*B22(j+1);
Vx2(counter,1)=C; %x-displacement
Vy2(counter,1)=C; %y-displacement
Vz2(counter,1)=C; %z-displacement
Vrx2(counter,1)=C; %x-rotation
Vry2(counter,1)=C; %y-rotation
end
end
Vx2=[Vx2;zero_vector;zero_vector;zero_vector;zero_vector];
Vy2=[zero_vector;Vy2;zero_vector;zero_vector;zero_vector];
Vz2=[zero_vector;zero_vector;Vz2;zero_vector;zero_vector];
Vrx2=[zero_vector;zero_vector;zero_vector;Vrx2;zero_vector];
Vry2=[zero_vector;zero_vector;zero_vector;zero_vector;Vry2];

Vxx1=diff(Vx1,x);
Vyy1=diff(Vy1,y);
Vxy1=diff(Vx1,y);
Vyx1=diff(Vy1,x);
Vzx1=diff(Vz1,x);
Vzy1=diff(Vz1,y);
Vrxx1=diff(Vrx1,x);
Vryy1=diff(Vry1,y);
Vrxy1=diff(Vrx1,y);
Vryx1=diff(Vry1,x);

Vxx2=diff(Vx2,x);
Vyy2=diff(Vy2,y);
Vxy2=diff(Vx2,y);
Vyx2=diff(Vy2,x);
Vzx2=diff(Vz2,x);
Vzy2=diff(Vz2,y);
Vrxx2=diff(Vrx2,x);
Vryy2=diff(Vry2,y);
Vrxy2=diff(Vrx2,y);
Vryx2=diff(Vry2,x);

%Strain matrices
B_epsilon_1=[Vxx1.’;Vyy1.’;Vxy1.’+Vyx1.’];
B_kappa_1=[Vrxx1.’;Vryy1.’;Vrxy1.’+Vryx1.’];
B_gamma_1=[Vzy1.’-Vry1.’;Vzx1.’-Vrx1.’];
B_epsilon_2=[Vxx2.’;Vyy2.’;Vxy2.’+Vyx2.’];
B_kappa_2=[Vrxx2.’;Vryy2.’;Vrxy2.’+Vryx2.’];
B_gamma_2=[Vzy2.’-Vry2.’;Vzx2.’-Vrx2.’];
B_a_1=1/(2*h_a)*[(-Vz1).’;(-Vx1).’+h_1*Vrx1.’;(-Vy1).’+h_1*Vry1.’];
B_a_2=1/(2*h_a)*[Vz2.’;Vx2.’+h_2*Vrx2.’;Vy2.’+h_2*Vry2.’];

B_1=[B_epsilon_1;B_kappa_1;B_gamma_1];
B_2=[B_epsilon_2;B_kappa_2;B_gamma_2];

C1=[A1,B1,zeros(3,2);B1,D1,zeros(3,2);zeros(2,6),G1];
C2=[A2,B2,zeros(3,2);B2,D2,zeros(3,2);zeros(2,6),G2];
Ca=[E_a,0,0;0,G_a,0;0,0,G_a];

%————————————
%!!! THIS IS WHERE THE PROBLEMS OCCUR
%————————————

%K_a_12
disp(‘K_a_12’)
tic
K_a_12_unintegrated=B_a_1.’*Ca*B_a_2;
K_a_12_unintegrated = reshape(K_a_12_unintegrated,1089,100);
K_a_12=sym(zeros(1089,100));
toc

fprintf(‘Commencing integration…n’);

parfor i=1:1089
tic
fprintf(‘currently integrating row %d’,i);
K_a_12(i,:)=2*h_a*int(int(K_a_12_unintegrated(i,:),y,0,W),x);
toc
end

return;I am doing a parallel computation, specifically a row by row integration of a symbolic matrix (sic I am integrating 4 rows at a time on my quadcore). 2 curious things occur: