MATLAB Tutorial Chapter 4. Advanced matrix operations 4.1. Sparse matrices SPARSE MATRICES To show the efficiency gained by using sparse matrices, we will solve a PDE using finite differences twice. First, we will use the matrix commands that use the full matrix that we have learned so far. Second, we will use new commands that take advantage of the fact that most of the elements are zero to greatly reduce both the memory requirements and the number of floating point operations required to solve the PDE. clear all; remove all existing variables from memory num_pts = 100; number of grid points in simulation CALCULATION WITH FULL MATRIX FORMAT The following matrix is obtained from using central finite differences to discretize the Laplacian operator in 1-D. x = 1:num_pts; grid of x-values Set the matrix from discretizing the PDE with a 1-D grid containing num_pts points with a spacing of 1 between points. Afull=zeros(100,100); Afull(1,1) = 1; Afull(num_pts,num_pts) = 1; for i=2:(num_pts-1) sum over interior points Afull(i,i) = 2; Afull(i,i-1) = -1; Afull(i,i+1) = -1; end Dirichlet boundary conditions at x=1 and x=num_pts are set. BC1 = -10; value of f at x(1); BC2 = 10; value of f at x(num_pts); For the interior points, we have a source term. b_RHS = linspace(0,0,num_pts)'; create column vector of zeros b_RHS(1) = BC1; b_RHS(num_pts) = BC2; b_RHS(2:(num_pts-1)) = 0.05; for interior, b_RHS is source term We now use the standard MATLAB solver to obtain the solution of the PDE at the grid points. f = Afull\b_RHS; figure; plot(x,f); title('PDE solution from FD-CDS method (full matrix)'); xlabel('x'); ylabel('f(x)'); Let us now take a closer look at Afull. The command spy(A) makes a plot of the matrix A by writing a point wherever an element of A has a non-zero value. figure; spy(Afull); title('Structure of Afull');The number nz at the bottom is the number of non-zero elements. We see that only a small fraction of the matrix elements are non-zero. Since we numbered the grid points in a regular manner with the neighbors of each grid point stored in adjacent locations, the non-zero elements in this matrix are on the principal diagonal and the two diagonals immediately above and below. Even if we numbered the grid points irregularly, we would still have this small number of non-zero points. It is often the case, as it is here, that the matrices we encounter in the numerical simulation of PDE's are sparse; that is, only a small fraction of their points are non-zero. For this matrix, the total number of elements is num_elements = num_pts*num_pts; nzA = nnz(Afull); returns # of non-zero elements in Afull fraction_filled = nzA/num_elements This means that Afull is mostly empty space and we are wasting a lot of memory to store values we know are zero. Remove all variables from memory except Afull. clear x f b_RHS BC1 BC2 i num_elements nzA fraction_filled; SPARSE MATRIX We can convert a matrix to sparse format using the command "sparse". Asparse = sparse(Afull) MATLAB stores a sparse matrix as an NZ by 3 array where NZ is the number of non-zero elements. The first column is the row number and the second the column number of the non-zero element. The third column is the actual value of the non-zero element. The total memory usage is far smaller than with the full matrix format. whos Afull Asparse; clear Asparse; get rid of sparse matrix NOW WE WILL SOLVE USING SPARSE MATRIX FORMAT Next, we set the grid point values x = 1:num_pts; grid of x-values Now we declare the matrix A to have sparse matrix structure from the start. First, we calculate the number of non-zero elements (or an upper bound to this number). We see that for each row corresponding to an interior point, we have 3 values, whereas for the first and last row we only have one value. Therefore, the number of non-zero elements is nzA = 3*(num_pts-2) + 2; We now use "spalloc(m,n,nz)" that allocates memory for a m by n dimensioned sparse matrix with no more than nz non-zero elements. A = spalloc(num_pts,num_pts,nzA); We now set the values of the A matrix. A(1,1) = 1; A(num_pts,num_pts) = 1; for i=2:(num_pts-1) A(i,i) = 2; A(i,i-1) = -1; A(i,i+1) = -1; endDirichlet boundary conditions at x=1 and x=num_pts are set. BC1 = -10; value of f at x(1); BC2 = 10; value of f at x(num_pts); For the interior points, we have a source term. b_RHS = linspace(0,0,num_pts)'; create column vector of zeros b_RHS(1) = BC1; b_RHS(num_pts) = BC2; b_RHS(2:(num_pts-1)) = 0.05; for interior, b_RHS is source term Now, when we call the MATLAB standard solver, it automatically identifies that A is a sparse matrix, and uses solver algorithms that take advantage of this fact. f = A\b_RHS; figure; plot(x,f); title('PDE solution from FD-CDS method (sparse matrix)'); xlabel('x'); ylabel('f(x)'); whos A Afull; From the lines for A and Afull, we can see that the sparse matrix format requires far less memory that the full matrix format. Also, if N is the number of grid points, we see that the size of the full matrix is N^2; whereas, the size in memory of the sparse matrix is only approximately 3*N. Therefore, as N increases, the sparse matrix format becomes far more efficient than the full matrix format. For complex simulations with thousands of grid points, one cannot hope to solve these problems without taking advantage of sparsity. To see the increase in execution speed that can be obtained by using sparse matrices, examine the following two algorithms for multiplying two matrices. FULL MATRIX ALGORITHM FOR MATRIX MULTIPLICATION Bfull = 2*Afull; Cfull = 0*Afull; declare memory for C=A*B num_flops = 0; for i=1:num_pts for j=1:num_pts for k=1:num_pts Cfull(i,j) = Cfull(i,j) + Afull(i,k)*Bfull(k,j); num_flops = num_flops + 1; end end end disp(['# FLOPS with full matrix format = ', int2str(num_flops)]); SPARSE MATRIX ALGORITHM FOR MATRIX MULTIPLICATION B = 2*A; nzB = nnz(B); # of non-zero elements of B nzC_max = round(1.2*(nzA+nzB)); guess how much memory we'll need for C C = spalloc(num_pts,num_pts,nzC_max); [iA,jA] = find(A); find (i,j) elements that are non-zero in A [iB,jB] = find(B); find (i,j) elements that are non-zero in B num_flops = 0; for ielA = 1:nzA iterate over A non-zero elements for ielB = 1:nzB iterate over B non-zero elements if(iB(ielB)==jA(ielA)) the pair contributes to C i = iA(ielA); k = jA(ielA);j = jB(ielB); C(i,j) = C(i,j) + A(i,k)*B(k,j); num_flops = num_flops + 1; end