1function ex2 (n) 2%EX2: create an n-by-n 2D mesh, four different ways 3 4% Example: 5% ex2 6% See also: cs_demo 7 8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com 9 10if (nargin < 1) 11 n = 30 ; 12end 13 14subplot (1,2,1) ; 15 16% method 1: create an n-by-n 2D mesh for the 2nd difference operator 17tic 18ii = zeros (5*n^2, 1) ; 19jj = zeros (5*n^2, 1) ; 20xx = zeros (5*n^2, 1) ; 21k = 1 ; 22for j = 0:n-1 23 for i = 0:n-1 24 s = j*n+i + 1 ; 25 ii (k:k+4) = [(j-1)*n+i j*n+(i-1) j*n+i j*n+(i+1) (j+1)*n+i ] + 1 ; 26 jj (k:k+4) = [s s s s s] ; 27 xx (k:k+4) = [-1 -1 4 -1 -1] ; 28 k = k + 5 ; 29 end 30end 31 32% remove entries beyond the boundary 33keep = find (ii >= 1 & ii <= n^2 & jj >= 1 & jj <= n^2) ; 34ii = ii (keep) ; 35jj = jj (keep) ; 36xx = xx (keep) ; 37A = sparse (ii,jj,xx) ; 38t1 = toc ; disp (t1) ; 39% subplot (2,2,1) ; 40spy (A) 41title (sprintf ('%d-by-%d 2D mesh\n', n, n)) ; 42 43% method 2, using no for loops 44tic 45nn = 1:n^2 ; 46i2 = [nn-n ; nn-1 ; nn ; nn+1 ; nn+n] ; 47j2 = repmat (nn, 5, 1) ; 48x2 = repmat ([-1 -1 4 -1 -1]', 1, n^2) ; 49keep = find (i2 >= 1 & i2 <= n^2 & j2 >= 1 & j2 <= n^2) ; 50i2 = i2 (keep) ; 51j2 = j2 (keep) ; 52x2 = x2 (keep) ; 53C = sparse (i2,j2,x2) ; 54t2 = toc ; disp (t2) ; 55 56% subplot (2,2,2) ; plot (j2) ; 57% title ('2D fast j2') ; 58disp (A-C) ; 59 60any (ii-i2) 61any (jj-jj) 62 63% method 3: create an n-by-n-by-n 3D mesh for the 2nd difference operator 64tic 65[A, keep, ii, jj, xx] = mesh3d1 (n) ; 66ii = ii (keep) ; 67jj = jj (keep) ; 68xx = xx (keep) ; 69t3 = toc ; disp (t3) ; 70tic 71E = sparse (ii,jj,xx) ; 72t3b = toc ; disp (t3b) ; 73subplot (1,2,2) ; spy (E) ; 74title (sprintf ('%d-by-%d-by-%d 3D mesh\n', n, n, n)) ; 75 76% method 4, using no for loops 77tic 78nn = 1:n^3 ; 79i2 = [nn-n^2 ; nn-n ; nn-1 ; nn ; nn+1 ; nn+n ; nn+n^2] ; 80j2 = repmat (nn, 7, 1) ; 81x2 = repmat ([-1 -1 -1 6 -1 -1 -1]', 1, n^3) ; 82keep = find (i2 >= 1 & i2 <= n^3 & j2 >= 1 & j2 <= n^3) ; 83i2 = i2 (keep) ; 84j2 = j2 (keep) ; 85x2 = x2 (keep) ; 86t4 = toc ; disp (t4) ; 87tic 88F = sparse (i2,j2,x2) ; 89t4b = toc ; disp (t4b) ; 90disp (E-F) ; 91