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 nonzero 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 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 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 …