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