# Linear algebra performance issue

**URL:** <https://discourse.julialang.org/t/linear-algebra-performance-issue/42353>\
**Category:** General Usage\
**Created:** [July 1, 2020, 6:48am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353 "2020-07-01T06:48:09Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 1, 2020, 6:48am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/1 "2020-07-01T06:48:09Z")

</div>

Hello all,

I need to compute a 3x3 matrix many many times, it is given by P = J x sigma x inv(F)‘, where J is a scalar, sigma and F are 3x3 matrices; A’ =\> transpose of A.

What is the best way to compute this P as fast as possible in Julia?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 1, 2020, 7:07am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/2 "2020-07-01T07:07:37Z")

</div>

Since the matrix size is very small, consider using `SMatrix{3,3}` from StaticArrays.jl

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 1, 2020, 7:12am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/3 "2020-07-01T07:12:47Z")

</div>

Thanks, I did use that.

```julia
@timeit "2" P = J*sigma*inv(F)' # convert to Piola Kirchoof stress

```

But the above way of computing P involves many allocation, as measured by TimerOutput.jl.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 1, 2020, 7:14am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/4 "2020-07-01T07:14:39Z")

</div>

If you use StaticArrays it should not allocate at all, if it does, there might be some type-stability issue in your code. It’s hard to tell without having an example to run.

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 1, 2020, 7:20am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/5 "2020-07-01T07:20:23Z")

</div>

No temporary allocations at all? If so, I will double check with @code\_warntype. Thanks.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 1, 2020, 7:24am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/6 "2020-07-01T07:24:24Z")

</div>

Hmm, ~~it does appear to allocate a small amount, I’m not sure why~~

```julia
julia> a
3×3 SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
  0.873545 1.1805 2.00785
  1.26515 2.59622 0.654212
 -1.51945 -1.70419 0.517809

julia> b
3×3 SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
  1.5845 1.17958 1.72889
 -0.974253 -1.15078 -0.647302
  1.12172 0.0169034 0.0549179

julia> @time sum(2*a/b)
  0.000004 seconds (3 allocations: 176 bytes)
0.2072010681879124

```

EDIT: sorry, my benchmarking was flawed, it does not allocate anything

```julia
julia> @btime sum(2*$a/$b')
  26.254 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 1, 2020, 7:41am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/7 "2020-07-01T07:41:09Z")

</div>

Thanks. I confirm your answer. No allocation.

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 1, 2020, 9:35am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/8 "2020-07-01T09:35:55Z")

</div>

Actually, I got allocations because one of these matrices is MMatrix, not SMatrix. The reason is:

```julia
sigma::MMatrix{Float64,3,3}
function modify_stress!(sigma, epsilon::SMatrix{3,3,Float64})
  sigma .= (epsilon[1,1]+epsilon[2,2]+epsilon[3,3]) * UniformScaling(1.) 
end

```

Then, sigma got modified by function modify\_stress!().  
The thing is I cannot use SMatrix for sigma because it does not allow: sigma .= … By writing sigma = (epsilon[1,1]+epsilon[2,2]+epsilon[3,3]) \* UniformScaling(1.), sigma did not get modified!!!

How could I solve this dilemma?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 1, 2020, 9:41am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/9 "2020-07-01T09:41:29Z")

</div>

You do not need to update sigma in place when you use static arrays, simply return the new one

```julia
function stress(epsilon::SMatrix{3,3,Float64})
  (epsilon[1,1]+epsilon[2,2]+epsilon[3,3]) * UniformScaling(1.) 
end

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [July 1, 2020, 10:55am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/10 "2020-07-01T10:55:16Z")

</div>

> [@vinhphunguyen](#):
>
> ```julia
> sigma::MMatrix{Float64,3,3}
> function modify_stress!(sigma, epsilon::SMatrix{3,3,Float64})
> sigma .= (epsilon[1,1]+epsilon[2,2]+epsilon[3,3]) * UniformScaling(1.) 
> end
> 
> ```

I’m confused. Does this work for you? For me it gives an error:

```julia
ERROR: MethodError: no method matching length(::UniformScaling{Float64})

```

What’s the output?

(BTW, are you aware of the `tr` function? You can use `tr(epsilon)` instead of `epsilon[1,1]+epsilon[2,2]+epsilon[3,3]`)

This works for me:

```julia
sigma .= tr(epsilon) * I(3)

```

But you should probably use `SMatrix` and just replace `sigma`, as @baggepinnen suggests.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 1, 2020, 3:51pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/11 "2020-07-01T15:51:30Z")

</div>

Note that you probably want to solve `P=J*sigma \ (F')` as it will be likely both faster and more accurate.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [July 1, 2020, 3:59pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/12 "2020-07-01T15:59:36Z")

</div>

> [@baggepinnen](#):
>
> You do not need to update sigma in place when you use static arrays, simply return the new one

To elaborate on this, it is usually better to write things in a functional way when dealing with StaticArrays. So instead of having `modify_stress!` where you pass in a return value you simply have `compute_stress` and you only pass in the strain. This makes the code less stateful and you pass less arguments to function and it usually get significantly faster.

Since I see you are doing continuum mechanics, you could also take a look at [GitHub - Ferrite-FEM/Tensors.jl: Efficient computations with symmetric and non-symmetric tensors with support for automatic differentiation.](https://github.com/KristofferC/Tensors.jl).

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 3, 2020, 4:36pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/14 "2020-07-03T16:36:54Z")

</div>

Yes, it works for me, my old code is:

```julia
function update_stress!(sigma::SMatrix{3,3,Float64},
						epsilon::SMatrix{3,3,Float64})
  #sigma .= mat.lambda * (epsilon[1,1]+epsilon[2,2]+epsilon[3,3]) * UniformScaling(1.) + 2.0 * mat.mu * epsilon

```

I did not know of the tr() function. Thanks.

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 3, 2020, 4:38pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/15 "2020-07-03T16:38:02Z")

</div>

Thanks a lot Kristoffer.

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 3, 2020, 4:39pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/16 "2020-07-03T16:39:28Z")

</div>

Great, will try that. thanks.

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 4, 2020, 3:01am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/17 "2020-07-04T03:01:13Z")

</div>

HI Kristogger,

Tensor.jl is pretty cool. I have may be a stupid question, for 3x3 matrices, Tensor.jl is faster than StaticArrays.jl? I did some tests and the speed seems to be the same. But I want to ask your expert opinion.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [July 5, 2020, 5:25pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/18 "2020-07-05T17:25:43Z")

</div>

Should be pretty much the same overall. For 3x3 it might be faster in scalar products with other second order tensors because we use explicit SIMD instructions whereas StaticArrays.jl need help from LLVM (and LLVM is bad for the 3x3 case)

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 10, 2020, 9:45am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/19 "2020-07-10T09:45:47Z")

</div>

Carlsson, could u please help with this code?

```julia
struct Solid3D_1
        F :: Vector{SMatrix{3,3,Float64,9}} # F, 2x2 matrix
        ϵ :: Vector{SMatrix{3,3,Float64,9}} # F, 2x2 matrix
        σ :: Vector{SMatrix{3,3,Float64,9}} # F, 2x2 matrix
        parcount::Int64
        function Solid3D_1(count)
            strain = fill(zeros(3,3),count)
            defor = fill(zeros(3,3),count)
            stress = fill(zeros(3,3),count)
            return new(defor,strain,stress,count)
        end
end

struct Solid3D_2
        F :: Vector{Tensor{2,3,Float64}} # Tensor{order,dim,T<:Real} F, 2x2 matrix
        ϵ :: Vector{SymmetricTensor{2,3,Float64}} # Tensor{order,dim,T<:Real} F, 2x2 matrix
        σ :: Vector{SymmetricTensor{2,3,Float64}} # F, 2x2 matrix
        parcount::Int64
        function Solid3D_2(count)
            defor = zeros(Tensor{2,3}, count) #fill(zero(Tensor{2, 3}), count)
            strain = zeros(SymmetricTensor{2, 3},count)
            stress = zeros(SymmetricTensor{2, 3},count)
            return new(defor,strain,stress,count)
        end
end

function main_static_arrays(solid)
        defor = solid.F
        epsi = solid.ϵ
        stres = solid.σ
        for ip = 1:solid.parcount
            J = det(defor[ip])
            Finv = inv(defor[ip]) # no memory alloc
            P = J*stres[ip]*Finv' # convert to Piola Kirchoof stress, no memory alloc
        end
end

function main_tensor(solid)
        defor = solid.F
        epsi = solid.ϵ
        stres = solid.σ
        for ip = 1:solid.parcount
            #J = det(defor[ip])
            Finv = inv(defor[ip]) # no memory alloc
            #P = J*stres[ip]⋅Finv' # convert to Piola Kirchoof stress, no memory alloc
        end
end

function main()
parCount = 100000
s1 = Solid3D_1(parCount)
s2 = Solid3D_2(parCount)
@btime main_static_arrays(s1)
@btime main_tensor(s2)
end

main()

```

I got:

```julia
julia> include("Main_Test_Tensor.jl")
  48.933 μs (0 allocations: 0 bytes)
  3.037 ms (100000 allocations: 7.63 MiB)

```

Why?

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [July 10, 2020, 10:08am UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/20 "2020-07-10T10:08:51Z")

</div>

From the docs example ([https://kristofferc.github.io/Tensors.jl/stable/man/constructing\_tensors/#Constructing-tensors-1](https://kristofferc.github.io/Tensors.jl/stable/man/constructing_tensors/#Constructing-tensors-1)), it looks like `SymmetricTensor{order,dim,T}` is incomplete type specification, so that the elements of `Solid3D_2` are vectors with abstract element types. The complete type signature must be `F :: Vector{Tensor{2,3,Float64,9}}` and `ϵ :: Vector{SymmetricTensor{2,3,Float64,6}}`, I guess.  
Try the following definition:

```julia
struct Solid3D_2{T<:Tensor{2,3,<:AbstractFloat}, ST<:SymmetricTensor{2,3,<:AbstractFloat}, IT<:Integer}
        F :: Vector{T} # Tensor{order,dim,T<:Real} F, 2x2 matrix
        ϵ :: Vector{ST} # Tensor{order,dim,T<:Real} F, 2x2 matrix
        σ :: Vector{ST} # F, 2x2 matrix
        parcount::IT
        function Solid3D_2(count)
            defor = zeros(Tensor{2,3}, count) #fill(zero(Tensor{2, 3}), count)
            strain = zeros(SymmetricTensor{2, 3},count)
            stress = zeros(SymmetricTensor{2, 3},count)
            T, ST = eltype.(defor, strain) 
            return new{T,ST,Int}(defor,strain,stress,count)
        end
end

```

---

<div class="post-metadata">

**Author:** ![vinhphunguyen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vinhphunguyen/32/16316_2.png) [@vinhphunguyen](https://discourse.julialang.org/u/vinhphunguyen)\
**Post date:** [July 10, 2020, 12:46pm UTC](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353/21 "2020-07-10T12:46:44Z")

</div>

Thanks Vasily, I used your suggestion, but not a parameterised struct:

```julia
function main_static_arrays(solid)
        defor = solid.F
        epsi = solid.ϵ
        stres = solid.σ
        for ip = 1:solid.parcount
            J = det(defor[ip])
            Finv = inv(defor[ip]) # no memory alloc
            P = J*stres[ip]*Finv' # convert to Piola Kirchoof stress, no memory alloc
            L = Finv*Finv 
            D = 0.5 * (L + L') # no memory alloc
        end
end

function main_tensor(solid)
        defor = solid.F
        epsi = solid.ϵ
        stres = solid.σ
        for ip = 1:solid.parcount
            J = det(defor[ip])
            Finv = inv(defor[ip]) # no memory alloc
            P = J*stres[ip]⋅Finv' # convert to Piola Kirchoof stress, no memory alloc
            L = Finv⋅Finv 
            D = symmetric(L) # no memory alloc
        end
end

function main1()
    parCount = 100000
    s1 = Solid3D_1(parCount)
    @btime main_static_arrays($s1)
end

function main2()
    parCount = 100000
    s2 = Solid3D_2(parCount)
    @btime main_tensor($s2)
end

main1()
main2()

```

```julia
julia> include("Main_Test_Tensor.jl")
  51.429 μs (0 allocations: 0 bytes)
  77.122 μs (0 allocations: 0 bytes)

```

Seems to me symmetric(A) is slow as I removed it, Tensors.jl was faster slightly.

[Next page](https://discourse.julialang.org/t/linear-algebra-performance-issue/42353.md?page=2)