-
Notifications
You must be signed in to change notification settings - Fork 0
/
resort_matrix_graph_column_sort.R
135 lines (102 loc) · 4 KB
/
resort_matrix_graph_column_sort.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
131
132
133
134
135
# Stephan Frickenhaus, Alfred Wegener Institute (2022)
#
# microbenchmark effects of matrix re-orderings on matrix vector products
# approach: use a precomputed pseudo-FE-linear-problem-matrix with random non-zeros;
# a reordering is generated by a SFC on the 2D layout of the adjacency matrix,
# and the matrix is reordered accordingly, to show identical numerics.
#
# HERE : with column index sorting in rows
#
# 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
rm(list=ls())
require(Rcpp)
sourceCpp("resortgrid_SFC_Rcpp.cpp",cacheDir = ".cpp-cache/")
require(SparseM)
load("Pseudo-FE-Matrix_m_rnd_1.5mio_n15_colsort.rdat")
require(Matrix)
# number of nodes
N=m@dimension[1]
# with static number of neighbors n::
n=m@ia[2]-m@ia[1]
#the matrix from benchars.R is transformed to a graph structure;
# we need the coordinate format; probably there is a faster way without coo-format
as.matrix.coo(m)->m.coo
TN <- sparseMatrix(dims=dim(m),i=m.coo@ia,j=m.coo@ja,x=1,use.last.ij = TRUE)
# now make an adjacency graph from the matrix
require(igraph)
system.time(ig <- graph.adjacency(TN,mode="undirected"))
# now the pivotMDS layout
require(graphlayouts)
system.time(g.layout<-layout_with_pmds(ig,dim = 2,pivot=15))
# plot(g.layout,pch=".")
# SFC from graph-layout-coordinates
x<-g.layout[,1]
y<-g.layout[,2]
mysfc(x,y) -> p
#plot(x,y,pch=19,cex=0.6)
#lines(x[p],y[p],col=3,lwd=0)
#points(x,y,pch=19,cex=0.6)
##
# benchmark this graph-layout SFC-reordering
#
mysfc(x,y) -> p
# need the the inverse
p.i<-(1:N)
p.i[p]<-1:N
p.i[p][1:10]
p[1:3]
# step 1
nn.m=t(matrix(m@ja,nr=n)) # retrieve Nxn nearest neighbors from original matrix
nn.m[1,]
nnp=nn.m[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(m@ra,ncol=n,byrow = TRUE)
# check correct assembly of values in first row
all(m@ra[1:n]==rap[1,])
rap1=rap[p,] # reorder rows of values
rap2=t(rap1)
dim(rap2)
v=rnorm(N)
m1=new("matrix.csr",as.numeric(rap2),as.integer(jap),as.integer(m@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
#Unit: milliseconds
#expr min lq mean median uq max neval
# { v.0 <- (m %*% v)[, 1] } 220.3754 239.8412 284.6610 253.1899 279.5407 868.7194 150
#{ v.1 <- (m1 %*% v1)[, 1] } 126.0684 136.0200 156.3108 144.6116 160.3212 423.0101 150
# compare with benchmark on mesh-resorting
#Unit: milliseconds
#expr min lq mean median uq max neval
#{ v.0 <- (m %*% v)[, 1] } 230.4816 257.3719 288.5977 275.8483 297.0142 499.4441 50
#{ v.1 <- (m1 %*% v1)[, 1] } 129.0277 140.7763 151.1027 150.9353 158.1047 187.3399 50
#-> similar speedup
#report on index locality by statistics of SD(col indices)
summary(apply(matrix(m@ja,nr=n),2,sd))
summary(apply(matrix(m1@ja,nr=n),2,sd))
#> summary(apply(matrix(m@ja,nr=n),2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#126100 393481 431423 429500 467697 645868
#> summary(apply(matrix(m1@ja,nr=n),2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#4.6 48.5 136.2 4060.2 553.6 654082.1
#compare to mesh-resorting benchmark
#> summary(apply(ja,2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#126100 393481 431423 429500 467697 645868
#> summary(apply(jap,2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#4.5 11.4 28.2 2911.7 126.2 774387.3
# -> less col-index locality