How to more efficiently populate a sparse matrix

I have a code that solves the laplace equation in axisymmetric coordinates in a cone-shaped domain using a finite difference method. The solution is then of the form

phi = A\b;

where A is my operator matrix and b are the source terms. My issue is concerned with creating the matrix A which is sparse and diagonal. My current approach is:

Nz = someConstant;                               % Number of intervals in z
Nr = someConstant;                               % Number of intervals in r
N = Nr*Nz;
A = sparse(N,N);                                 % N rows, N columns
b = sparse(N,1);                                 % N rows, 1 column

for ii = 2:Nr-1
    for jj = 2:Nz-1
        n = ii + (jj - 1)*Nr;                    % Chosen indexing convention
        r = R(ii,jj);                            % Array containing radial coordinates
        dr = DR(jj);                             % Vector containing dr terms
        dz = someConstant;
          
        A(n, n     ) = -2*(1/dz^2+1/dr^2);       % \phi_{i,j}
        A(n, n - 1 ) = 1/dr^2 - 1/(2*dr*r);      % \phi_(i-1,j)
        A(n, n + 1 ) = 1/dr^2 + 1/(2*dr*r);      % \phi_(i+1,j)
        A(n, n - Nr) = 1/dz^2;                   % \phi_(i, j-1)
        A(n, n + Nr) = 1/dz^2;                   % \phi_(i, j+1)
        
        b(n, 1     ) = 0;                        % Source term 
    end
end

Matlab warns me that 'this sparse indexing expression is likely to be slow' when filling in the terms in A but I don't know how to improve this. The particular concern is that because my domain is a cone shape dr depends on the z index.

(This code snippit does not include the boundary conditions so b = zeros(N,1) which will not give a solution)

Answer 1

[Bjorn Gustavsson]

The general advice is to calculate and save the non-zero elements, and their array-indices, in 3 arrays and then once they've been calculated create your sparse array. Something like this (not your example converted):

Avals = zeros(5*n^2,1);
idx1 = Avals;
idx2 = idx1;
for i1 = 2:n
    for i2 = 2:n
        Avals(i2+n*(i1-1)+0) = -4; % Centre-point of Laplacian differentiation-operator
        idx1(i2+n*(i1-1)+0) = i1;
        idx2(i2+n*(i1-1)+0) = i2;
        Avals(i2+n*(i1-1)+1) = 1; % One of the nearest neighbors
        idx1(i2+n*(i1-1)+0) = i1-1; % ...etc
        idx2(i2+n*(i1-1)+0) = i2;
    end
end
A = sparse(idx1,idx2,Avals,n,n);

That way you only create space for your sparse array once, and the index and magnitude-arrays are also pre-allocated.

(In my experience I should test these generation-snippets with small enough arrays to manually inspect them.)

HTH

[Toysey]

Thank you very much - that is exactly what I needed to know