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