# How to cook with SIMD?
**URL:** https://discourse.haskell.org/t/how-to-cook-with-simd/10224
**Category:** Learn
**Created:** [August 25, 2024, 2:47pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224 "2024-08-25T14:47:09Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@wiz](https://discourse.haskell.org/u/wiz)
#### Post date: [August 25, 2024, 2:47pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/1 "2024-08-25T14:47:09Z")
I’ve been trying to make my rendering code a bit faster with wide ops, but everything I try ends up slower than the dumbest straightforward variant.
E.g. pointwise (+) for `Vec4 ByteArray#` is a tiny bit faster with FFIng, but loses by a third (!) to `Vec4 Float Float Float Float` and then using regular (+) per component.
Replacing scalar maths with the new FMA primops made it slower too.
I’ve found what appears to be an ideal target for this - 4-way ray-AABB intersection.
The thing is ~6.5x slower than a Haskell binary tree of AABBs with ordinary math for intersections.
I would be happy with the “explanation” that I simply suck at this or the partitioning is suboptimal. However…
Adding insult to injury, the profiler tells me that the 4-way version is indeed faster than 2-way:
| | 2-way | 4-way + SIMD | |
| --- | --- | --- | --- |
| Entries | 3085m | 951m | Yay, wide trees, less work |
| Time | 95% | 68.7% | Smaller portion of CPU resources, should be faster, right? |
| Alloc | 60.1% | 3.4% | 1/20 of allocations, gonna be good |
| Real time | 8s | 55s | …excuse me? |
---
### Author:  [@atravers](https://discourse.haskell.org/u/atravers)
#### Post date: [August 25, 2024, 9:09pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/2 "2024-08-25T21:09:52Z")
Out of a (web)search for `haskell ghc simd poor performance`, this:
> [@Optimiser performance problems](http://discourse.haskell.org/t/optimiser-performance-problems/8906):
>
> I’ve written the infinity norm version of the Adam gradient descent optimiser below. It should work with most multivariate convex problems. It seems to work ok except it is epically slow. Optimising x^2 - x with a bad starting point (e.g. -2000) as shown below in the test1 function takes about 9 seconds. Is there something I can do to make this orders of magnitude faster? I don’t see any straightforward, SIMD options, but I could be wrong? Perhaps a more Haskell friendly way to restate the comp…
seems to be the most recent of the first six _“vaguely-relevant”_ results - does anything there help?
---
### Author:  [@Bodigrim](https://discourse.haskell.org/u/Bodigrim)
#### Post date: [August 25, 2024, 9:31pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/3 "2024-08-25T21:31:19Z")
@wiz could you show the code? What kind of SIMD are we talking about?
---
### Author:  [@eldritch-cookie](https://discourse.haskell.org/u/eldritch-cookie)
#### Post date: [August 26, 2024, 12:23pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/4 "2024-08-26T12:23:04Z")
my understanding is that aside from LLVM pigbacking we have no SIMD primops without them any “SIMD” is actually just normal instructions simulating SIMD which will at best be as fast as not wrapping the primitives
---
### Author:  [@wiz](https://discourse.haskell.org/u/wiz)
#### Post date: [August 26, 2024, 2:52pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/5 "2024-08-26T14:52:07Z")
Exhibit 1: The basics.
```haskell
data Vec4 = Vec4 ByteArray#
```
```haskell
void Vec4_plus_Vec4_SIMD(float *x, float *y, float *o) {
_m128 xv = _mm_load_ps(x);
_m128 yv = _mm_load_ps(y);
xv = _mm_add_ps(xv, yv);
_mm_store_ps(o, xv);
}
```
Didn’t had much hope for this one, but it did made things tiiiiiiny bit faster (15.5ns → 15.0ns).
Now, behold the power of Haskell!
```haskell
data Vec4' = Vec4' {-# UNPACK #-} Float {-# UNPACK #-} Float {-# UNPACK #-} Float {-# UNPACK #-} Float
plus (Vec4' a0 a1 a2 a3) (Vec4 b0 b1 b2 b3) = Vec4' (a0 + b0) (a1 + b1) (a2 + b2) (a3 + b3)
```
This only takes 10ns.
Exhibit 2, something more normal.
```haskell
data Vec3 = Vec3 Float Float Float -- unpacked ofc
normalize v@(Vec3 x y z) =
if nearZero q || nearZero (abs (1 - q))
then v
else Vec3 (x * r) (y * r) (z * r)
where
q = dot v v
r = 1 / sqrt q
nearZero a = a <= 1e-6
```
12.5ns
```haskell
void Vec3_normalize_SIMD(float x, float y, float z, float *O) {
simde__m128 v4 = simde_mm_set_ps(x, y, z, 0); // init
simde__m128 q4 = simde_mm_dp_ps(v4, v4, 0xFF); // dot
float q = simde_mm_cvtss_f32(q4); // extract to scalar for testing
if (q <= 1e-6 || fabsf(1.0f - q) <= 1e-6) {
simde_mm_store_ps(O, v4); // bail out
} else {
simde__m128 n4 = simde_mm_rsqrt_ss(q4); // r
simde__m128 x = simde_mm_shuffle_ps(n4, n4, 0); // broadcast r
simde_mm_store_ps(O, simde_mm_mul_ps(v4, x)); // multiply and store
}
}
```
Getting multiple floats out required some wrapping/unwrapping
```haskell
normalize3 :: Vec3 -> Vec3
normalize3 v = withVec4 out \a' b' c' _ -> vec3 a' b' c'
where
out = unsafePerformIO do
result@(Vec4 o) <- unsafeNewVec4
withVec3 v \a b c -> vec3_normalize a b c o
pure result
```
35.2ns. Okay, with so much dancing around the representation it is kinda expected.
Exhibit 3: brb… gotta collect more data
---
### Author:  [@wiz](https://discourse.haskell.org/u/wiz)
#### Post date: [August 26, 2024, 2:56pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/6 "2024-08-26T14:56:03Z")
Well, my data is already unboxed and stuff. The SIMD was dismissed early in the thread and I want to know specifically how to actually _use_ it.
---
### Author:  [@Bodigrim](https://discourse.haskell.org/u/Bodigrim)
#### Post date: [August 26, 2024, 3:00pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/7 "2024-08-26T15:00:55Z")
@wiz does “ps” stand for picoseconds? If yes then something is wrong with your benchmarks, 10 ps is unrealistically fast.
---
### Author:  [@wiz](https://discourse.haskell.org/u/wiz)
#### Post date: [August 26, 2024, 7:40pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/8 "2024-08-26T19:40:15Z")
Those are nanoseconds indeed.
---
### Author:  [@Bodigrim](https://discourse.haskell.org/u/Bodigrim)
#### Post date: [August 26, 2024, 7:58pm UTC](https://discourse.haskell.org/t/how-to-cook-with-simd/10224/9 "2024-08-26T19:58:31Z")
My understanding is that to call opaque foreign code Haskell has to save / restore all CPU registers. This is lots of work, more than addition of four `Float`s, explaining Exhibit 1 and, at least partially, Exhibit 2.
A loop over an array of `Float`s should reside in C code for vectorised instructions to pay off.