内容简介:In this algebra, it’s possible to compute on a mesh of arbitrary 5 dimensionalAdditionally, inProgramming the
In this algebra, it’s possible to compute on a mesh of arbitrary 5 dimensional conformal geometric algebra simplices, which can be represented by a bundle of nested dyadic tensors.
julia> using Grassmann, StaticArrays; basis"+-+++" (⟨+-+++⟩, v, v₁, v₂, v₃, v₄, v₅, v₁₂, v₁₃, v₁₄, v₁₅, v₂₃, v₂₄, v₂₅, v₃₄, v₃₅, v₄₅, v₁₂₃, v₁₂₄, v₁₂₅, v₁₃₄, v₁₃₅, v₁₄₅, v₂₃₄, v₂₃₅, v₂₄₅, v₃₄₅, v₁₂₃₄, v₁₂₃₅, v₁₂₄₅, v₁₃₄₅, v₂₃₄₅, v₁₂₃₄₅) julia> value(rand(Chain{V,1,Chain{V,1}})) 5-element StaticArrays.SArray{Tuple{5},Chain{⟨+-+++⟩,1, ,253} where 253 where ,1,5} with indices SOneTo(5): -0.33459594357756073v₁ - 0.3920064467082769v₂ - 0.575474920388841v₃ + 0.6150287650825268v₄ - 0.7568209093000915v₅ -0.7402635950699139v₁ - 0.9303076362461833v₂ + 0.9729806462365271v₃ - 0.8514563480551867v₄ + 0.09906887873006287v₅ -0.7456570397821101v₁ - 0.6497560949330929v₂ + 0.4756585550844967v₃ - 0.31169948016530347v₄ - 0.9355928499340793v₅ -0.4555014543082292v₁ + 0.712268225360094v₂ - 0.7500443783398549v₃ - 0.36349628003234713v₄ + 0.5005769037801056v₅ -0.07402971220645727v₁ + 0.19911765119918146v₂ - 0.4980618129231722v₃ - 0.7728564279829087v₄ + 0.9165735719353756v₅ julia> A = Chain{V,1}(rand(SMatrix{5,5})) (0.9244801277294266v₁ + 0.029444337884018346v₂ + 0.745495522394158v₃ + 0.6695874677964055v₄ + 0.4998003712198389v₅)v₁ + (0.5423877012973404v₁ + 0.30112324458605655v₂ + 0.9530587650033631v₃ + 0.2706004745612134v₄ + 0.37762612797501616v₅)v₂ + (0.7730171467954035v₁ + 0.019660709510785912v₂ + 0.39119534821037494v₃ + 0.9403026278575068v₄ + 0.07545094732793833v₅)v₃ + (0.7184128110093908v₁ + 0.6295740775044767v₂ + 0.5179035493253021v₃ + 0.039081667453648716v₄ + 0.3719284661613145v₅)v₄ + (0.5033657705978616v₁ + 0.41183905359914386v₂ + 0.7761548051732969v₃ + 0.07635587137916744v₄ + 0.5582934197259402v₅)v₅
Additionally, in Grassmann.jl we prefer the nested usage of pure ChainBundle
parametric types for large re-usable global cell geometries, from which local dyadics can be selected.
Programming the A\b
method is straight forward with some Julia language metaprogramming and Grassmann.jl by first instantiating some Cramer symbols
Base.@pure function Grassmann.Cramer(N::Int) x,y = SVector{N}([Symbol(:x,i) for i ∈ 1:N]),SVector{N}([Symbol(:y,i) for i ∈ 1:N]) xy = [:(($(x[1+i]),$(y[1+i])) = ($(x[i])∧t[$(1+i)],t[end-$i]∧$(y[i]))) for i ∈ 1:N-1] return x,y,xy end
These are exterior product variants of the Cramer determinant symbols ( N! times N -simplex hypervolumes), which can be combined to directly solve a linear system:
@generated function Base.:\(t::Chain{V,1,<:Chain{V,1}},v::Chain{V,1}) where V N = ndims(V)-1 # paste this into the REPL for faster eval x,y,xy = Grassmann.Cramer(N) mid = [:($(x[i])∧v∧$(y[end-i])) for i ∈ 1:N-1] out = Expr(:call,:SVector,:(v∧$(y[end])),mid...,:($(x[end])∧v)) return Expr(:block,:((x1,y1)=(t[1],t[end])),xy..., :(Chain{V,1}(getindex.($(Expr(:call,:./,out,:(t[1]∧$(y[end])))),1)))) end
Which results in the following highly efficient @generated
code for solving the linear system,
(x1, y1) = (t[1], t[end]) (x2, y2) = (x1 ∧ t[2], t[end - 1] ∧ y1) (x3, y3) = (x2 ∧ t[3], t[end - 2] ∧ y2) (x4, y4) = (x3 ∧ t[4], t[end - 3] ∧ y3) Chain{V, 1}(getindex.(SVector(v ∧ y4, (x1 ∧ v) ∧ y3, (x2 ∧ v) ∧ y2, (x3 ∧ v) ∧ y1, x4 ∧ v) ./ (t[1] ∧ y4), 1))
Benchmarks with that algebra indicate a 3x faster performance than SMatrix
for applying A\b
to bundles of dyadic elements.
julia> @btime $(rand(SMatrix{5,5},10000)).\Ref($(SVector(1,2,3,4,5))); 2.588 ms (29496 allocations: 1.44 MiB) julia> @btime $(Chain{V,1}.(rand(SMatrix{5,5},10000))).\$(v1+2v2+3v3+4v4+5v5); 808.631 μs (2 allocations: 390.70 KiB) julia> @btime $(SMatrix(A))\$(SVector(1,2,3,4,5)) 150.663 ns (0 allocations: 0 bytes) 5-element SArray{Tuple{5},Float64,1,5} with indices SOneTo(5): -4.783720495603508 6.034887114999602 1.017847212237964 6.379374861538397 -4.158116538111051 julia> @btime $A\$(v1+2v2+3v3+4v4+5v5) 72.405 ns (0 allocations: 0 bytes) -4.783720495603519v₁ + 6.034887114999605v₂ + 1.017847212237964v₃ + 6.379374861538393v₄ - 4.1581165381110505v₅
Such a solution is not only more efficient than Julia’s StaticArrays.jl method for SMatrix
, but is also useful to minimize allocations in Grassmann.jl finite element assembly.
Similarly, the Cramer symbols can also be manipulated to invert the linear system or for determining whether a point is within a simplex.
julia> using Grassmann; basis"3" (⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃) julia> T = Chain{V,1}(Chain(v1),v1+v2,v1+v3) (1v₁ + 0v₂ + 0v₃)v₁ + (1v₁ + 1v₂ + 0v₃)v₂ + (1v₁ + 0v₂ + 1v₃)v₃ julia> barycenter(T) ∈ T, (v1+v2+v3) ∈ T (true, false)
There are multiple equivalent ways of computing the same results using the ⋅
and :
dyadic products.
julia> T\barycenter(T) == inv(T)⋅barycenter(T) true julia> sqrt(T:T) == norm(SMatrix(T)) true
The following Makie.jl streamplot was generated with the Grassmann.Cramer
method and interpolated from Nedelec edges of a Maxwell finite element solution.
More info about these examples is at https://grassmann.crucialflow.com/dev/tutorials/dyadic-tensors
Hermann Grassmann was the inventor of linear algebra as we know it today.
以上所述就是小编给大家介绍的《Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl》,希望对大家有所帮助,如果大家有任何疑问请给我留言,小编会及时回复大家的。在此也非常感谢大家对 码农网 的支持!
猜你喜欢:本站部分资源来源于网络,本站转载出于传递更多信息之目的,版权归原作者或者来源机构所有,如转载稿涉及版权问题,请联系我们。
算法导论
[美] Thomas H. Cormen、Charles E. Leiserson、Ronald L. Rivest、Clifford Stein / 高等教育出版社 / 2002-5 / 68.00元
《算法导论》自第一版出版以来,已经成为世界范围内广泛使用的大学教材和专业人员的标准参考手册。 这本书全面论述了算法的内容,从一定深度上涵盖了算法的诸多方面,同时其讲授和分析方法又兼顾了各个层次读者的接受能力。各章内容自成体系,可作为独立单元学习。所有算法都用英文和伪码描述,使具备初步编程经验的人也可读懂。全书讲解通俗易懂,且不失深度和数学上的严谨性。第二版增加了新的章节,如算法作用、概率分析......一起来看看 《算法导论》 这本书的介绍吧!