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)
采纳的回答
Bjorn Gustavsson
2020-12-2
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
