How to efficiently update a matrix in Parallel?

6 Set 2016
1 Risposta
Aggiornato 6 Set 2016
8 Visualizzazioni (30 giorni)

0 voti

I am looking to update a large sparse matrix (of size N^2 by M^2, currently about 9mil by 9mil, and looking to go larger).
What's the best way to leverage parallel computing to do this? (I'm trying to implement it using parallel cpu cores first, and eventually push it onto the GPU if possible). The matrix elements are defined by a smaller matrix of size NxM. Some of the elements of the NxM matrix will randomly change as i step through time, and when that happens, I want to update only the elements of the N^2xM^2 elements that change.
It seems like the most common ways of doing this involve breaking the matrix into blocks, or by column. However that seems like a lot of overhead, and because the matrix is so big i'm worried about memory issues, or wasting time on elements that don't need to be updated.
edit: Copied code below, and attached m file for ease of reading (formatting went wonky when i copy/pasted)

Accedi per rispondere a questa domanda.

Risposte (1)

James Tursa
James Tursa il 6 Set 2016

0 voti

In general, anytime you change elements in a sparse matrix to/from zero to non-zero, the entire matrix must be copied in memory. This is due to how the matrix elements are stored in memory. So your goal should be to limit the number of times this is done, which is going to depend on exactly how you are updating the elements. Not having seen your code for this, I can only offer the suggestion to gather all of your updates off to the side somewhere (e.g., save the row & column & value data in separate variables), and then apply all of these updates at once to your sparse matrix just prior to needing to use it. I.e., try to do the updates entirely in one function call so that (hopefully) only one data copy will take place.

3 Commenti
Mostra 1 commento meno recente Nascondi 1 commento meno recente

Joshua L.
Joshua L. il 6 Set 2016
Modificato: James Tursa il 6 Set 2016
The style of the matrix is such that i'm only ever updating non zero elements to another non zero, I shouldn't have to change anything to/from 0.
I'm hoping to have to avoid calling spdiags to rebuild it each time, but that is a possibility (For some reason, it seems to take more memory than the matrix, so i could push to larger grid sizes if i updated after creating it with spdiags. I'm not sure why)
I've attached my current code below. Nothing is parallelized yet, since i was struggling to make it work. I attached the full thing for clarity, but the part I'm attempting to optimise starts at "Trying to parallelize here " (equivalent to line 148). edit: attempting to fix the formatting, apparently my commenting interfered with the web versions commenting. Everything below is 1 program (I also attached the m file)
parpool(4);
g=@(x,y) ((node_col(x,1)+node_col(y,1))/2 );
R_metal=single(1);
R_ins=single (1000);
True_I= 1*10^(-5);
cond_metal= single (1/R_metal);
cond_ins=single(1/R_ins);
X=1000;
Y=1000;
disp 'building matrix..'
tic
A=ones(X*Y+2,5);
A(:,1)=-single(cond_ins);
A(:,2)=-single(cond_ins);
A(:,3)=single(5*cond_ins);
A(1,3)=single(20*cond_ins);
A(X*Y+2,3)=single(20*cond_ins);
A(2,3)=single(4*cond_ins);
A(1+X:X*Y+1-X,3)=single(4*cond_ins);
A(X*Y+1,3)=single(4*cond_ins);
A(X*Y+1-X,3)=single(4*cond_ins);
A(:,4)=-single(cond_ins);
A(:,5)=-(cond_ins);
A(X+2,3)=single(3*cond_ins);
A(X+1,3)=single(4*cond_ins);
A(X+1,2)=single(0);
A(X+2,4)=single(0);
for j=2:Y-1
A(j*X+2,3)=single(3*cond_ins);
if j==Y-1
A(j*X+2,3)=single(4*cond_ins);
end
A(j*X+1,3)=single(3*cond_ins);
A(j*X+1,2)=single(0);
A(j*X+2,4)=single(0);
end
d=[-X -1 0 1 X];
G=spdiags(A,d,X*Y+2,X*Y+2);
clear A
for i=2:X+1
G(1,i)=-single ( 2* cond_ins);
G(i,1)=-single ( 2* cond_ins);
G(X*Y+2,X*Y+3-i)=-single ( 2* cond_ins);
G(X*Y+3-i,X*Y+2)=-single ( 2* cond_ins);
end
toc
disp 'finished building original matrix'
I= sparse(X*Y+1,1);
I(1,1)=single(1);
I(X*Y+2,1)=single(-1);
V=G \ I;
R_eq1=(V(1,1)-V(X*Y+1,1)) * 100
rng(0,'twister');
T_start=365;
T_end=370;
T_stepsize=0.1;
iter=(T_end-T_start)/(T_stepsize);
sigmasq = 10;
Tc_avg = 335.5;
temp=T_start;
node_grid= R_ins .* ones(X,Y);
node_col=R_ins .* ones(X*Y,1);
switch_grid=zeros(X,Y);
switch_col=zeros(X*Y,1);
Tc_grid = sigmasq.*randn(X,Y) + Tc_avg;
R_eq_vec=zeros(iter,1);
for temp_step=1:iter
temp=temp+ T_stepsize
disp 'time to update G'
tic
for ii=1:X
for jj=1:Y
if node_grid(ii,jj) == R_ins
if temp > Tc_grid(ii,jj)
node_grid(ii,jj)=R_metal;
node_col(jj+X*(ii-1),1)=node_grid(ii,jj);
switch_grid(ii,jj)=1;
switch_col(jj+X*(ii-1),1)=switch_grid(ii,jj);
end
end
end
end
n=0;
nodeindex2=0;
spmd;
node_col=R_ins .* ones(X*Y,1);
for i=1:X
for j=1:Y
if switch_col(j+X*(i-1),1)==1
n=(j+(i-1)*X);
%check if first row
if n <= X
%check if left edge
if n==1
G(n+1,n+1)=1./g(n,n+1)+1./g(n,n+X)+2./node_col(n,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n+X+1)+2./node_col(n+1,1);
G(n+1+X,n+1+X)=1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
G(1,n+1)=-2./node_col(n,1);
G(n+1,1)=G(1,n+1);
G(n+1,n+2)=-1./g(n,n+1);
G(n+2,n+1)=-1./g(n+1,n);
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
end
%check right edge
if n==X
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n+X)+2./node_col(n,1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n+X-1)+2./node_col(n-1,1);
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
G(1,n+1)=-2./node_col(n,1);
G(n+1,1)=G(1,n+1);
G(n+1,n+2)=0;
G(n+2,n+1)=0;
G(n+1,n)=-1./g(n,n-1);
G(n,n+1)=-1./g(n-1,n);
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
end
%neither edge
if n>1 && n < X
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n+1)+1./g(n,n+X)+2./node_col(n,1);
if n == 2
G(n,n)=1./g(n-1,n)+1./g(n-1,n+X-1)+2./node_col(n-1,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n+X+1)+2./node_col(n+1,1);
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n == X-1
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+X+1)+2./node_col(n+1,1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n+X-1)+2./node_col(n-1,1);
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n > 2 && n < X-1
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n+X-1)+2./node_col(n-1,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n+X+1)+2./node_col(n+1,1);
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
G(1,n+1)=-2./node_col(n,1);
G(n+1,1)=G(1,n+1);
G(n+1,n+2)=-1./g(n,n+1);
G(n+2,n+1)=-1./g(n+1,n);
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
end
updatenum=0;
for iupdate=2:X+1
updatenum=updatenum+G(iupdate,1);
end
G(1,1)=-updatenum;
end
%check if last row
if n > X*Y-X
%check if left edge
if n == X*Y-X+1
G(n+1,n+1)=1./g(n,n+1)+1./g(n,n-X)+2./node_col(n,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+2./node_col(n+1,1);
G(n+1-X,n+1-X)=1./g(n-X,n+1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
G(X*Y+2,n+1)=-2./node_col(n,1);
G(n+1,X*Y+2)=G(n+1,X*Y+2);
G(n+1,n+2)=0;
G(n+2,n+1)=0;
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);
G(n,n+1)=-1./g(n+1,n);
G(n+1,n)=-1./g(n+1,n);
end
%check if right edge
if n == X*Y
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n-X)+2./node_col(n,1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+2./node_col(n-1,1);
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
G(X*Y+2,n+1)=-2./node_col(n,1);
G(n+1,X*Y+2)=G(X*Y+2,n+1);
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
end
%not edges
if n > X*Y-X+1 && n < X*Y
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n+1)+1./g(n,n-X)+2./node_col(n,1);
if n == X*Y-X+2
G(n,n)=1./g(n-1,n)+1./g(n-1,n-X-1)+2./node_col(n-1,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+2./node_col(n+1,1);
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
end
if n == X*Y-1
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n-X+1)+2./node_col(n+1,1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+2./node_col(n-1,1);
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
end
if n > X*Y-X+2 && n < X*Y-1
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+2./node_col(n-1,1);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+2./node_col(n+1,1);
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
end
G(X*Y+2,n+1)=-2./node_col(n,1);
G(n+1,X*Y+2)=-2./node_col(n,1);
G(n+1,n+2)=-1./g(n,n+1);
G(n+2,n+1)=-1./g(n+1,n);
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);;
end
updatenum2=0;
for iupdate2=2:X+1
updatenum2=updatenum2+G(iupdate2,1);
end
G(X*Y+2,X*Y+2)=-updatenum2;
end
%not first or last rows
if n > X && n < X*Y-X+1
%check if a left edge
if mod(n,X)==1
G(n+1,n+1)=1./g(n,n+1)+1./g(n,n-X)+1./g(n,n+X);
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+1./g(n+1,n+X+1);
if n < X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n == X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n+1+X)+1./g(n+X,n)+2./node_col(n+X,1);
end
if n > 2*X
G(n+1-X,n+1-X)=1./g(n-X,n+1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
end
if n ==2*X
G(n+1-X,n+1-X)=1./g(n-X,n+1-X)+1./g(n-X,n)+2./node_col(n-X,1);
end
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
G(n+1,n+2)=-1./g(n,n+1);
G(n+2,n+1)=-1./g(n+1,n);
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);
G(n,n+1)=0;
G(n+1,n)=0;
end
%check if right edge
if mod(n,X)== 0
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n-X)+1./g(n,n+X);
if n < X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n == X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n)+2./node_col(n+X,1);
end
if n > 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n)+1./g(n-X,n-X-X);
end
if n ==2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n)+2./node_col(n-X,1);
end
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+1./g(n-1,n+X-1);
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
G(n+1,n+2)=0;
G(n+2,n+1)=0;
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);
end
%not edges
if mod(n,X)~= 0 && mod(n,X) ~=1
G(n+1,n+1)=1./g(n,n-1)+1./g(n,n+1)+1./g(n,n-X)+1./g(n,n+X);
if mod(n,X)==2
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+1./g(n+1,n+X+1);
G(n,n)=1./g(n-1,n)+1./g(n-1,n-X-1)+1./g(n-1,n+X-1);
if n < X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n == X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+2./node_col(n+X,1);
end
if n > 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n-X-X)+1./g(n-X,n);
end
if n == 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+2./node_col(n-X,1);
end
end
if mod(n,X)==X-1
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n-X+1)+1./g(n+1,n+X+1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+1./g(n-1,n+X-1);
if n < X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n <= X*Y-X && n >= X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+2./node_col(n+X,1);
end
if n > 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n-X-X)+1./g(n-X,n);
end
if n > X && n <= 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+2./node_col(n-X,1);
end
end
if mod(n,X) > 2 && mod(n,X)< X-1
G(n+2,n+2)=1./g(n+1,n)+1./g(n+1,n+2)+1./g(n+1,n-X+1)+1./g(n+1,n+X+1);
G(n,n)=1./g(n-1,n-2)+1./g(n-1,n)+1./g(n-1,n-X-1)+1./g(n-1,n+X-1);
if n < X*Y-2*X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+1./g(n+X,n+X+X);
end
if n > X*Y-2*X && n <= X*Y-X
G(n+1+X,n+1+X)=1./g(n+X,n-1+X)+1./g(n+X,n+1+X)+1./g(n+X,n)+2./node_col(n+X,1);
end
if n > 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n-X-X)+1./g(n-X,n);
end
if n > X && n < 2*X
G(n+1-X,n+1-X)=1./g(n-X,n-1-X)+1./g(n-X,n+1-X)+1./g(n-X,n)+2./node_col(n-X,1);
end
end
G(n,n+1)=-1./g(n-1,n);
G(n+1,n)=-1./g(n,n-1);
G(n+1,n+2)=-1./g(n,n+1);
G(n+2,n+1)=-1./g(n+1,n);
G(n+1+X,n+1)=-1./g(n+X,n);
G(n+1,n+1+X)=-1./g(n,n+X);
G(n+1-X,n+1)=-1./g(n-X,n);
G(n+1,n+1-X)=-1./g(n,n-X);
end
end
end
end
end
end
toc
for ii=1:X
for jj=1:Y
if switch_grid(ii,jj) == 1;
switch_grid(ii,jj)=2;
end
end
end
disp 'time to invert'
tic
V= G \ I;
R_eq= (V(1,1)-V(X*Y+1,1)) * 100;
R_eq_vec(temp_step,1)=R_eq;
toc
end
tempaxis=[T_start:T_stepsize:T_end]';
R_eq_vec=[R_eq1,R_eq_vec'];
semilogy(tempaxis,R_eq_vec);
delete(gcp('nocreate'));
James Tursa
James Tursa il 6 Set 2016
Modificato: James Tursa il 6 Set 2016
I can't find "Trying to parallelize here" anywhere in your code so I do not know where to look.
Also, MATLAB does not have a sparse single type. It only has sparse double and sparse logical types. So this code probably doesn't do what you think it does:
A=ones(X*Y+2,5); % <-- A is double
A(:,1)=-single(cond_ins); % <-- A is still double
:
G=spdiags(A,d,X*Y+2,X*Y+2); % <-- G is sparse double, not sparse single
All of those assignments to A elements with single(etc) stuff on the rhs of the assignment do nothing to affect the type of A. All of those single(etc) values on the rhs simply get converted to double as they are assigned to elements of A, which is still double. (Same comment applies to element assignments to G and I).
You could do this to make A single:
A=ones(X*Y+2,5,'single'); % <-- A is single
But then you would run into problems converting it directly to sparse. E.g.,
>> sparse(ones(2,3,'single'))
Undefined function 'sparse' for input arguments of type 'single'.
Joshua L.
Joshua L. il 6 Set 2016
Modificato: Joshua L. il 6 Set 2016
Thank you for the warning, i realized afterwards that it wasn't doing anything. I ended up just leaving it in for now, but Matlab's sparse matrix requires double, I need to go through and remove those
Sorry about the "Trying to parallelize here", i put it in the comment, but not the m file (oops). I added it now, it's lines 148-489. It's long because there are a bunch of edge cases (Where i'd go out of bounds), but really all it is doing is updating G using values in node_col, which is fixed for each step in temperature

Accedi per commentare.

Categorie

Scopri di più su Matrix Indexing in Centro assistenza e File Exchange

Prodotti

Tag

Richiesto:

Joshua L.
Joshua L.
il 6 Set 2016

Modificato:

Joshua L.
Joshua L.
il 6 Set 2016

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by