Skip to content

Commit 9009796

Browse files
JordiManyerJordiManyer
authored
Merge pull request #1059 from Antoinemarteau/master
optimized MonomialBases evaluations
2 parents b3cfd04 + c4ab198 commit 9009796

File tree

6 files changed

+279
-57
lines changed

6 files changed

+279
-57
lines changed

‎NEWS.md

+3
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
1717
- Fixed #974, an error when weak form is real but unknown vector is complex. Since PR[#1050](https://github.com/gridap/Gridap.jl/pull/1050).
1818
- Fixed issue where barycentric refinement rule in 3D would not produce oriented meshes. Since PR[#1055](https://github.com/gridap/Gridap.jl/pull/1055).
1919

20+
### Changed
21+
- Optimized MonomialBasis low-level functions. Since PR[#1059](https://github.com/gridap/Gridap.jl/pull/1059).
22+
2023
## [0.18.7] - 2024-10-8
2124

2225
### Added

‎benchmark/Project.toml

+1
Original file line numberDiff line numberDiff line change
@@ -2,3 +2,4 @@
22
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
33
Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
44
PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
5+
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

‎benchmark/benchmarks.jl

+1
Original file line numberDiff line numberDiff line change
@@ -12,3 +12,4 @@ end
1212
const SUITE = BenchmarkGroup()
1313

1414
@include_bm SUITE "bm_assembly"
15+
@include_bm SUITE "bm_monomial_basis"

‎benchmark/bm/bm_monomial_basis.jl

+190
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
module bm_monomial_basis
2+
3+
using PkgBenchmark, BenchmarkTools
4+
using Gridap
5+
using Gridap.Polynomials
6+
using Gridap.TensorValues
7+
using StaticArrays
8+
9+
################################################
10+
# src/Polynomials/MonomialBasis.jl: _set_value_!
11+
################################################
12+
13+
gradient_type = Gridap.Fields.gradient_type
14+
15+
_set_value! = Gridap.Polynomials._set_value!
16+
17+
function set_value_driver(f,T,D,x,n)
18+
k = 1
19+
s = one(T)
20+
for i in 1:n
21+
k = f(x,s,k)
22+
end
23+
end
24+
25+
function set_value_benchmarkable(D, T, V, n)
26+
C = num_indep_components(V)
27+
x = zeros(V,n*C)
28+
return @benchmarkable set_value_driver($_set_value!,$T,$D,$x,$n)
29+
end
30+
31+
##################################################
32+
# src/Polynomials/ModalC0Bases.jl: _set_value_mc0!
33+
##################################################
34+
35+
_set_value_mc0! = Gridap.Polynomials._set_value_mc0!
36+
37+
function set_value_mc0_driver(f,T,D,x,n)
38+
k = 1
39+
s = one(T)
40+
for i in 1:n
41+
k = f(x,s,k,2)
42+
end
43+
end
44+
45+
function set_value_mc0_benchmarkable(D, T, V, n)
46+
C = num_indep_components(V)
47+
x = zeros(V,2*n*C)
48+
return @benchmarkable set_value_mc0_driver($_set_value_mc0!,$T,$D,$x,$n)
49+
end
50+
51+
###################################################
52+
# src/Polynomials/MonomialBasis.jl: _set_gradient!
53+
###################################################
54+
55+
 _set_gradient! = Gridap.Polynomials. _set_gradient!
56+
57+
function set_gradient_driver(f,T,D,V,x,n)
58+
k = 1
59+
s = VectorValue{D,T}(ntuple(_->one(T),D))
60+
for i in 1:n
61+
k = f(x,s,k,V)
62+
end
63+
end
64+
65+
function set_gradient_benchmarkable(D, T, V, n)
66+
C = num_indep_components(V)
67+
G = gradient_type(V, zero(Point{D,T}))
68+
x = zeros(G,n*C);
69+
return @benchmarkable set_gradient_driver($_set_gradient!,$T,$D,$V,$x,$n)
70+
end
71+
72+
#####################################################
73+
# src/Polynomials/ModalC0Bases.jl: _set_gradient_mc0!
74+
#####################################################
75+
76+
 _set_gradient_mc0! = Gridap.Polynomials. _set_gradient_mc0!
77+
78+
function set_gradient_mc0_driver(f,T,D,V,x,n)
79+
k = 1
80+
s = VectorValue{D,T}(ntuple(_->one(T),D))
81+
for i in 1:n
82+
k = f(x,s,k,1,V)
83+
end
84+
end
85+
86+
function set_gradient_mc0_benchmarkable(D, T, V, n)
87+
C = num_indep_components(V)
88+
G = gradient_type(V, zero(Point{D,T}))
89+
x = zeros(G,n*C);
90+
return @benchmarkable set_gradient_mc0_driver($_set_gradient_mc0!,$T,$D,$V,$x,$n)
91+
end
92+
93+
#################################################
94+
# src/Polynomials/MonomialBasis.jl: _evaluate_1d!
95+
#################################################
96+
97+
_evaluate_1d! = Gridap.Polynomials._evaluate_1d!
98+
99+
function evaluate_1d_driver(f,order,D,v,x_vec)
100+
for x in x_vec
101+
f(v,x,order,D)
102+
end
103+
end
104+
105+
function evaluate_1d_benchmarkable(D, T, V, n)
106+
n = Integer(n/50)
107+
order = num_indep_components(V)
108+
v = zeros(D,order+1);
109+
x = rand(MVector{n,T})
110+
return @benchmarkable evaluate_1d_driver($_evaluate_1d!,$order,$D,$v,$x)
111+
end
112+
113+
################################################
114+
# src/Polynomials/MonomialBasis.jl:_gradient_1d!
115+
################################################
116+
117+
_gradient_1d! = Gridap.Polynomials._gradient_1d!
118+
119+
function gradient_1d_driver(f,order,D,v,x_vec)
120+
for x in x_vec
121+
f(v,x,order,D)
122+
end
123+
end
124+
125+
function gradient_1d_benchmarkable(D, T, V, n)
126+
n = Integer(n/10)
127+
order = num_indep_components(V)
128+
v = zeros(D,order+1);
129+
x = rand(MVector{n,T})
130+
return @benchmarkable gradient_1d_driver($_gradient_1d!,$order,$D,$v,$x)
131+
end
132+
133+
################################################
134+
# src/Polynomials/MonomialBasis.jl:_hessian_1d!
135+
################################################
136+
137+
_hessian_1d! = Gridap.Polynomials._hessian_1d!
138+
139+
function hessian_1d_driver(f,order,D,v,x_vec)
140+
for x in x_vec
141+
f(v,x,order,D)
142+
end
143+
end
144+
145+
function hessian_1d_benchmarkable(D, T, V, n)
146+
n = Integer(n/10)
147+
order = num_indep_components(V)
148+
v = zeros(D,order+1);
149+
x = rand(MVector{n,T})
150+
return @benchmarkable hessian_1d_driver($_hessian_1d!,$order,$D,$v,$x)
151+
end
152+
153+
#####################
154+
# benchmarkable suite
155+
#####################
156+
157+
const SUITE = BenchmarkGroup()
158+
159+
const benchmarkables = (
160+
set_value_benchmarkable,
161+
set_value_mc0_benchmarkable,
162+
set_gradient_benchmarkable,
163+
set_gradient_mc0_benchmarkable,
164+
evaluate_1d_benchmarkable,
165+
gradient_1d_benchmarkable,
166+
hessian_1d_benchmarkable
167+
)
168+
169+
const dims=(1, 2, 3, 5, 8)
170+
const n = 3000
171+
const T = Float64
172+
173+
for benchable in benchmarkables
174+
for D in dims
175+
TV = [
176+
VectorValue{D,T},
177+
TensorValue{D,D,T,D*D},
178+
SymTensorValue{D,T,Integer(D*(D+1)/2)},
179+
SymTracelessTensorValue{D,T,Integer(D*(D+1)/2)}
180+
]
181+
182+
for V in TV
183+
if V == SymTracelessTensorValue{1,T,1} continue end # no dofs
184+
name = "monomial_basis_$(D)D_$(V)_$(benchable)"
185+
SUITE[name] = benchable(D, T, V, n)
186+
end
187+
end
188+
end
189+
190+
end # module

‎src/Polynomials/ModalC0Bases.jl

+27-22
Original file line numberDiff line numberDiff line change
@@ -393,16 +393,10 @@ end
393393

394394
@inline function _set_value_mc0!(v::AbstractVector{V},s::T,k,l) where {V,T}
395395
ncomp = num_indep_components(V)
396-
m = zero(MVector{ncomp,T})
397396
z = zero(T)
398-
js = 1:ncomp
399-
for j in js
400-
for i in js
401-
@inbounds m[i] = z
402-
end
403-
@inbounds m[j] = s
404-
i = k+l*(j-1)
405-
@inbounds v[i] = Tuple(m)
397+
for j in 1:ncomp
398+
m = k+l*(j-1)
399+
@inbounds v[m] = ntuple(i -> ifelse(i == j, s, z),Val(ncomp))
406400
end
407401
k+1
408402
end
@@ -466,25 +460,36 @@ end
466460
# Indexing and m definition should be fixed if G contains symmetries, that is
467461
# if the code is optimized for symmetric tensor V valued FESpaces
468462
# (if gradient_type(V) returned a symmetric higher order tensor type G)
469-
@inline function _set_gradient_mc0!(
463+
@inline @generated function _set_gradient_mc0!(
470464
v::AbstractVector{G},s,k,l,::Type{V}) where {V,G}
465+
# Git blame me for readable non-generated version
471466
@notimplementedif num_indep_components(G) != num_components(G) "Not implemented for symmetric Jacobian or Hessian"
472-
473-
T = eltype(s)
474-
m = zero(Mutable(G))
475-
w = zero(V)
476-
z = zero(T)
477-
for (ij,j) in enumerate(CartesianIndices(w))
467+
468+
m = Array{String}(undef, size(G))
469+
N_val_dims = length(size(V))
470+
s_size = size(G)[1:end-N_val_dims]
471+
472+
body = "T = eltype(s); z = zero(T);"
473+
for ci in CartesianIndices(s_size)
474+
id = join(Tuple(ci))
475+
body *= "@inbounds s$id = s[$ci];"
476+
end
477+
478+
V_size = size(V)
479+
for (ij,j) in enumerate(CartesianIndices(V_size))
478480
for i in CartesianIndices(m)
479-
@inbounds m[i] = z
481+
m[i] = "z"
480482
end
481-
for i in CartesianIndices(s)
482-
@inbounds m[i,j] = s[i]
483+
for ci in CartesianIndices(s_size)
484+
id = join(Tuple(ci))
485+
m[ci,j] = "s$id"
483486
end
484-
i = k+l*(ij-1)
485-
@inbounds v[i] = m
487+
body *= "i = k + l*($ij-1);"
488+
body *= "@inbounds v[i] = ($(join(tuple(m...), ", ")));"
486489
end
487-
k+1
490+
491+
body = Meta.parse(string("begin ",body," end"))
492+
return Expr(:block, body ,:(return k+1))
488493
end
489494

490495
function _hessian_nd_mc0!(

0 commit comments

Comments
 (0)