-
Notifications
You must be signed in to change notification settings - Fork 0
/
benchmarks.R
130 lines (101 loc) · 3.59 KB
/
benchmarks.R
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
# Stephan Frickenhaus, Alfred Wegener Institute (2022)
#
# microbenchmark effects of matrix re-orderings on matrix vector products
# approach: use approx. Nearest Neighbor to create a
# pseudo-FE-linear-problem-matrix with random non-zeros;
# a reordering is generated by a SFC on the 2D random mesh nodes,
# and the matrix is reordered accordingly, to show identical numerics.
#
# SFC-code is based on:
# Rakowsky N. & Fuchs A. Efficient local resorting techniques with space filling curves applied
# to the tsunami simulation model TsunAWI. In IMUM 2011 - The 10th International Workshop
# on Multiscale (Un-)structured Mesh Numerical Modelling for coastal, shelf and global ocean
# dynamics, August 2011. http://hdl.handle.net/10013/epic.39576.d001
#install.packages("SparseM")
#install.packages("RANN")
require(RANN)
require(SparseM)
require(Rcpp)
require(microbenchmark)
sourceCpp("resortgrid_SFC_Rcpp.cpp",cacheDir = ".cpp-cache")
# run this under WSL or Linux
#require(parallel)
#mcaffinity(3) # choose core to bind R-execution to
N=1500000 # 2d nodes
n=15 # entries per row, non-symmetric
set.seed(10)
x=rnorm(N)
y=rnorm(N)
#plot(x,y,pch=".",type="l")
# create a matrix from n nearest neighbors
# typical for Differential equations linear problems
d=data.frame(x,y)
system.time(nn2(d,k=n)$nn.idx -> nn)
#for (i in sample(1:N,1000)) cat(i," ",length(intersect(order((x[i]-x)^2+(y[i]-y)^2)[1:n],nn[i,])),"\n")
#for (i in sample(1:N,1000)) cat(i," ",all(order((x[i]-x)^2+(y[i]-y)^2)[1:n]==nn[i,]),"\n")
dim(nn)
nn[1,]
nnz=prod(dim(nn))
ja=t(nn) # row-wise assembly in a vector
ja[1:(2*n)]
length(ja)-nnz
ia=(1:(N+1))*n-(n-1)
#ia[1:10]
set.seed(10)
ra=rnorm(nnz,mean=1.0)
m=new("matrix.csr",ra,as.integer(ja),as.integer(ia),dimension=as.integer(c(N,N)))
v=rnorm(N)
# resorting by p: p[i] is target index
mysfc(x,y) -> p
p[1:5]
# and the inverse
p.i<-(1:N)
p.i[p]<-1:N
p.i[p][1:10]
p[1:3]
# step 1
nnp=nn[p,] # reorder neighbor blocks (rows)
dim(nnp)
# step 2:
# transform indices in each block
nnp2=t(apply(nnp,1,function(x) p.i[x]))
dim(nnp2)
nnp2[1,]
nnp2[2,]
# nbrp2 is ready
jap=t(nnp2) # transpose to assemble in vector
rap=matrix(ra,ncol=n,byrow = TRUE)
# check correct assembly of values in first row
all(ra[1:n]==rap[1,])
rap1=rap[p,] # reorder rows of values
rap2=t(rap1)
dim(rap2)
#x[p] -> xp
#y[p] -> yp
m1=new("matrix.csr",as.numeric(rap2),as.integer(jap),as.integer(ia),dimension=as.integer(c(N,N)))
v1=v[p]
(mb0=microbenchmark({v.0<-(m%*%v)[,1]},{v.1<-(m1%*%v1)[,1]},times=150))
# N=1,5e6, n=15, notebook I7 , 12 core
#Unit: milliseconds
#expr min lq mean median uq max neval
#{ v.0 <- (m %*% v)[, 1] } 223.2206 236.2728 258.8170 247.0924 261.8256 462.8765 150
#{ v.1 <- (m1 %*% v1)[, 1] } 124.4876 133.4425 145.2893 138.7435 147.5078 316.1145 150
mb<-microbenchmark({v.0<-(m%*%v)[,1]},times=150)
mb1<-microbenchmark({v.1<-(m1%*%v1)[,1]},times=150)
save(m,file="Pseudo-FE-Matrix_m_rnd_1.5mio_n15.rdat")
save(m1,file="Pseudo-FE-Matrix_m1_sfc_1.5mio_n15.rdat")
sum(abs(v.0[p]-v.1))
plot(mb0)
rbind(mb,mb1)
# locality in ja index (standard dev of neighbour indices)
summary(apply(ja,2,sd))
summary(apply(jap,2,sd))
#N=1,5e6, n=15
#> summary(apply(ja,2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#106120 393411 431442 429502 467814 645354
#> summary(apply(jap,2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#4.5 11.4 28.1 2953.0 125.9 774489.3
plot((apply(ja,2,sd)),
(apply(jap,2,sd)),pch=".")