-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathresort_matrix_graph.R
131 lines (97 loc) · 3.79 KB
/
resort_matrix_graph.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
# Stephan Frickenhaus, Alfred Wegener Institute (2022)
#
# microbenchmark effects of matrix re-orderings on matrix vector products
# approach: use a precmputed 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.
#
# 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.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] } 221.6114 241.2720 260.9262 252.2645 267.5610 397.1865 150
#{ v.1 <- (m1 %*% v1)[, 1] } 120.0219 135.4608 147.8438 142.9620 150.6134 308.3793 150
# compare to benchmarks.R
#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
#-> 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.
#106120 393411 431442 429502 467814 645354
#> summary(apply(matrix(m1@ja,nr=n),2,sd))
#Min. 1st Qu. Median Mean 3rd Qu. Max.
#4.5 51.0 145.8 4089.4 607.1 692766.0
# comparison benchmarks.R:
#> 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
# -> less col-index locality