How to efficiently update a matrix in Parallel?

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)

Answer 1

[James Tursa]

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.

[Joshua L.]

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]

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.]

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