Re: SIMD layout again

"Yuri Davidovsky (as work at disclosure dot ie)" <[email protected]> Wed, 15 Jul 2026 17:13:15 +0200
Newsgroups gmane.lisp.lispworks.general
Message-ID <[email protected]>
--Apple-Mail=_DCD8F2CB-91F0-41DA-B717-7DBE47CE663E
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain;
	charset=utf-8


> On 15 Jul 2026, at 12:24, Tim Bradshaw (as tfb at tfeb dot org) =
<[email protected]> wrote:
>=20
> This is a question for people who understand SIMD things, which I =
don't.  There's a big preamble, sorry.
> ...
> The arrays have a rank one greater than the rank of the object rather =
than a 'stretched' last dimension because this makes handling rank-zero =
objects not have to be a weird special case.  Note that this doesn't =
alter the layout in memory of course.

I chose a different approach, I essentially implemented scalar objects =
and vectors of objects as separate entities. A scalar object is just a =
single typed-aref vector that accesses all its data linearly by

(sys:typed-aref type object field-offset)

An array of objects is also a single typed-aref vector, but its elements =
are accessed via the usual ay + x mechanics for a non-SIMD layout:

(sys:typed-aref type object (+ (* object-index object-size)
                               field-offset))

Essentially it is the same as accessing an element in a 2D array as in =
your approach, only here it is up to you to calculate its position in =
the memory. The reason why I kept both cases separate is because I =
actually do not need scalar objects at all, and added them just for =
completeness and a starting point when I began writing the code.

> =E2=80=A6 I can=E2=80=99t imagine a case where you would not just want =
them to be static, because you'll always know they're going to be passed =
to foreign code.

This is true, but if you declare things static, they had better persist =
throughout the lifetime of the application. If you allocate and =
deallocate things in the static area all the time, I would imagine you =
may end up wasting memory as the GC is likely limited in what it can do =
managing the static area of the heap.

> I think that SIMD things might well want the layout to be, for =
instance x0,x1,x2,x3,y0,y1,y2,y3,... =20

Yes, the homogenous fields need to be kept together since you typically =
apply the same operation to those, for example, when increasing some =
objects=E2=80=99 X-axis velocities by a specific delta.

> And note that you now have to make the length (or last dimension) be a =
multiple of 4, or you get a continuable error:

This is something to keep in mind when working with SIMD aligned data, =
ideally it should be a multiple of the SIMD register width. If it does =
not, you will need to mop up the tail of the array using scalar =
arithmetics (here is a C vector multiplication example for M platform):


#include <arm_neon.h> // SIMD intrinsics declarations
#include <stddef.h>   // size_t type declaration, used for memory =
addressing

void vector_multiply_neon(const float* a, // vector of floats A
                          const float* b, // vector of floats B
                          float* result,  // result vector
                          size_t size) {  // number of elements to =
process

    /// vector element index counter, we declare it separately
    /// because we will be using it in two loops
    size_t i =3D 0;=20

    /// pre-calculate the stopping point for the 4-way vector loop
    /// (size & ~3u) clears the lowest 2 bits, rounding size down to
    /// nearest multiple of 4, by flipping bits in a number 3 by ~3
    /// (0000 0011) -> (1111 1100) and ANDs it with the vector size
    size_t limit =3D size & ~(size_t)3;
   =20
    /// process 4 floats at a time (128-bit vectors)
    for (; i < limit; i +=3D 4) {

        /// load 4 floats from arrays A and B into SIMD registers VA and =
VB
	/// (two cycles)
        float32x4_t va =3D vld1q_f32(&a[i]); // vector load f32x4
        float32x4_t vb =3D vld1q_f32(&b[i]); // vector load f32x4
       =20
        /// multiply the SIMD registers VA and VB
        /// (single cycle)
        float32x4_t vres =3D vmulq_f32(va, vb);
       =20
        /// store the 4 result floats back to memory in RESULT vector
        /// (single cycle)
        vst1q_f32(&result[i], vres);
    }
   =20
    /// mop up the remaining elements (0 to 3 elements) using scalar =
math
    for (; i < size; ++i)
        result[i] =3D a[i] * b[i]; // four cycles per iteration
}

If you use SIMD libraries (like the previously mentioned Apple=E2=80=99s =
vDSP bundled with MacOS), you likely won=E2=80=99t need to worry about =
that, but you may need to worry about that when accessing things in =
Lisp. For example, typed-aref vectors in LW are always aligned on 8 =
bytes so you can=E2=80=99t rely on their lengths alone to identify how =
many objects you are holding there (unless they are a multiple of 8, or =
course). You will need to store the object count separately and pass it =
around.

> Lisp access to objects with permuted indices will be significantly =
slower, but that's not the point: they exist so that foreign SIMD code =
can eat them in a format it likes.

That is also the case with my approach, I estimated that accessing a =
field in a SIMD aligned struct vector requires 7 cycles (as opposed to =
just 2 in scalar alignment), here is an example for a single-float =
custom object property access:

1. Identify the SIMD frame in which the object lives (1 cycle):

(setq frame-index (ash object-index (- simd-frame-size-in-bits)))

2. Identify the offset of the SIMD frame in the vector (1 cycle):

(setq frame-offset (* frame-index simd-frame-size))

3. Identify the position of the object in the frame (1 cycle):

(setq object-in-frame-pos (- object-index frame-offset)

3. Identify the location of the field in its field sequence (1 cycle):

(setq field-pos (* object-in-frame-pos field-size)

4. Identify the actual field location in the SIMD frame by skipping over =
other field sequences before it  (1 cycle):

(setq field-offset (+ field-sequence-offset field-pos))

5. Now we can identify the actual location of the field in the vector  =
(1 cycle):

(setq field-index (+ field-offset frame-offset))

5. Finally we can get our single float field value (1 cycle):

(sys:typed-aref =E2=80=99single-float object-vector field-index)

If all this looks fairly complicated, it is because it is, do not worry =
if it goes over your head, it certainly did over mine. This is one of =
those things that you would need to write yourself from scratch to =
figure out.

> So the question is: before I redo the macro level (which is going to =
be fiddly), is something like this likely to be helpful for SIMD code =
(or any other foreign code purpose).

Yes, as long as the resulting data is aligned in a=20

((x=E2=82=80, x=E2=82=81, x=E2=82=82, x=E2=82=83)
 (y=E2=82=80, y=E2=82=81, y=E2=82=82, y=E2=82=83))

manner then there will be no issues with it on the foreign side. I did =
find the syntax

...
 ((x :array t :type fixnum :initform 12 :index-mappings ((2 1) (4 2) (6 =
3)))
  (y :array t :type fixnum :initform 13 :index-mappings ((1 4) (3 5) (5 =
6)))
...

somewhat confusing first, i.e.=20

x=E2=82=81 -> v[1], x=E2=82=82 -> v[2], x=E2=82=83 -> v[3]
y=E2=82=80 -> v[4], y=E2=82=81 -> v[5], y=E2=82=82 -> v[6]

wasn=E2=80=99t very obvious first, but since it is write once and can =
also be macro-automated if needed, I am not seeing any reasons why it =
wouldn=E2=80=99t work.








--Apple-Mail=_DCD8F2CB-91F0-41DA-B717-7DBE47CE663E
Content-Transfer-Encoding: quoted-printable
Content-Type: text/html;
	charset=utf-8

<html aria-label=3D"message body"><head><meta http-equiv=3D"content-type" =
content=3D"text/html; charset=3Dutf-8"></head><body =
style=3D"overflow-wrap: break-word; -webkit-nbsp-mode: space; =
line-break: after-white-space;"><br><blockquote type=3D"cite">On 15 Jul =
2026, at 12:24, Tim Bradshaw (as tfb at tfeb dot org) =
&lt;[email protected]&gt; wrote:<br><br =
class=3D"Apple-interchange-newline">This is a question for people who =
understand SIMD things, which I don't. &nbsp;There's a big preamble, =
sorry.</blockquote><blockquote type=3D"cite">...</blockquote><blockquote =
type=3D"cite">The arrays have a rank one greater than the rank of the =
object rather than a 'stretched' last dimension because this makes =
handling rank-zero objects not have to be a weird special case. =
&nbsp;Note that this doesn't alter the layout in memory of =
course.<br></blockquote><div><br></div><div>I chose a different =
approach, I essentially implemented scalar objects and vectors of =
objects as separate entities. A scalar object is just a single =
<i>typed-aref</i> vector that accesses all its data linearly =
by</div><div><br></div><div><b><font face=3D"Courier =
New">(sys:typed-aref type object&nbsp;</font></b><b><font face=3D"Courier =
New">field-offset)</font></b></div><div><br></div><div>An array of =
objects is also a single typed-aref vector, but its elements are =
accessed via the usual <i>ay + x</i> mechanics for a non-SIMD =
layout:</div><div><br></div><div><b><font face=3D"Courier =
New">(sys:typed-aref type object (+ (* object-index =
object-size)</font></b></div><div><b><font face=3D"Courier New">&nbsp; =
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; =
&nbsp; &nbsp; &nbsp; &nbsp; =
&nbsp;field-offset))</font></b></div><div><br></div><div>Essentially it =
is the same as accessing an element in a 2D array as in your approach, =
only here it is up to you to calculate its position in the memory. The =
reason why I kept both cases separate is because I actually do not need =
scalar objects at all, and added them just for completeness and a =
starting point when I began writing the code.</div><br><blockquote =
type=3D"cite">=E2=80=A6 I can=E2=80=99t imagine a case where you would =
not just want them to be static, because you'll always know they're =
going to be passed to foreign =
code.<br></blockquote><div><br></div><div>This is true, but if you =
declare things static, they had better persist throughout the lifetime =
of the application. If you allocate and deallocate things in the static =
area all the time, I would imagine you may end up wasting memory as the =
GC is likely limited in what it can do managing the static area of the =
heap.</div><br><blockquote type=3D"cite">I think that SIMD things might =
well want the layout to be, for instance x0,x1,x2,x3,y0,y1,y2,y3,... =
&nbsp;</blockquote><div><br></div><div>Yes, the homogenous fields need =
to be kept together since you typically apply the same operation to =
those, for example, when increasing some objects=E2=80=99 X-axis =
velocities by a specific delta.</div><br><blockquote type=3D"cite">And =
note that you now have to make the length (or last dimension) be a =
multiple of 4, or you get a continuable =
error:<br></blockquote><div><br></div><div>This is something to keep in =
mind when working with SIMD aligned data, ideally it should be a =
multiple of the SIMD register width. If it does not, you will need to =
mop up the tail of the array using scalar arithmetics (here is a C =
vector multiplication example for M =
platform):</div><div><br></div><div><br></div><div><div><font =
face=3D"Courier New"><b>#include &lt;arm_neon.h&gt; // SIMD intrinsics =
declarations</b></font></div><div><font face=3D"Courier New"><b>#include =
&lt;stddef.h&gt; &nbsp; // size_t type declaration, used for memory =
addressing</b></font></div><div><font face=3D"Courier =
New"><b><br></b></font></div><div><font face=3D"Courier New"><b>void =
vector_multiply_neon(const float* a, // vector of floats =
A</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; =
&nbsp; const float* b, // vector of floats B</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; =
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; float* result, &nbsp;// result =
vector</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; =
&nbsp; size_t size) { &nbsp;// number of elements to =
process</b></font></div><div><font face=3D"Courier =
New"><b><br></b></font></div><div><b style=3D"font-family: &quot;Courier =
New&quot;;">&nbsp; &nbsp; /// vector element index counter, we declare =
it separately</b></div><div><b style=3D"font-family: &quot;Courier =
New&quot;;">&nbsp; &nbsp; /// because we will be using it in two =
loops</b></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; size_t i =
=3D 0;&nbsp;</b></font></div><div><span class=3D"hljs-comment"><font =
face=3D"Courier New"><b><br></b></font></span></div><div><font =
face=3D"Courier New"><b><span class=3D"hljs-comment">&nbsp; &nbsp; /// =
pre-calculate the stopping point for the 4-way vector =
loop</span></b></font></div><div><font face=3D"Courier New"><b><span =
class=3D"hljs-comment">&nbsp; &nbsp; /// (size &amp; ~3u) clears the =
lowest 2 bits, rounding size down to</span></b></font></div><div><font =
face=3D"Courier New"><b><span class=3D"hljs-comment">&nbsp; &nbsp; /// =
nearest multiple of 4, by flipping bits in a number 3 by =
~3</span></b></font></div><div><font face=3D"Courier New"><b><span =
class=3D"hljs-comment">&nbsp; &nbsp; /// (0000 0011) -&gt; (1111 1100) =
and ANDs it with the vector size</span></b></font></div><div><font =
face=3D"Courier New"><b><span class=3D"hljs-keyword">&nbsp; &nbsp; =
size_t</span><span style=3D"font-size: medium;">&nbsp;limit =3D size =
&amp; ~(</span><span class=3D"hljs-keyword">size_t</span><span =
style=3D"font-size: medium;">)</span><span =
class=3D"hljs-number">3</span><span style=3D"font-size: =
medium;">;</span></b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp;&nbsp;</b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp; /// process 4 floats at a time (128-bit =
vectors)</b></font></div><div><font face=3D"Courier New"><b>&nbsp; =
&nbsp; for (; i &lt; limit; i +=3D 4) {</b></font></div><div><font =
face=3D"Courier New"><b><br></b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp; &nbsp; &nbsp; /// load 4 floats from arrays A and =
B into SIMD registers VA and VB</b></font></div><div><font face=3D"Courier=
 New"><b><span class=3D"Apple-tab-span" style=3D"white-space:pre">	=
</span>/// (two cycles)</b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp; &nbsp; &nbsp; float32x4_t va =3D =
vld1q_f32(&amp;a[i]); // vector load f32x4</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; &nbsp; float32x4_t vb =3D =
vld1q_f32(&amp;b[i]); // vector load f32x4</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; =
&nbsp;&nbsp;</b></font></div><div><font face=3D"Courier New"><b>&nbsp; =
&nbsp; &nbsp; &nbsp; /// multiply the SIMD registers VA and =
VB</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
&nbsp; &nbsp; /// (single cycle)</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; &nbsp; float32x4_t vres =3D =
vmulq_f32(va, vb);</b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; &nbsp; /// store the 4 =
result floats back to memory in RESULT vector</b></font></div><div><font =
face=3D"Courier New"><b>&nbsp; &nbsp; &nbsp; &nbsp; /// (single =
cycle)</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
&nbsp; &nbsp; vst1q_f32(&amp;result[i], =
vres);</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
}</b></font></div><div><font face=3D"Courier New"><b>&nbsp; =
&nbsp;&nbsp;</b></font></div><div><font face=3D"Courier New"><b>&nbsp; =
&nbsp; /// mop up the remaining elements (0 to 3 elements) using scalar =
math</b></font></div><div><font face=3D"Courier New"><b>&nbsp; &nbsp; =
for (; i &lt; size; ++i)</b></font></div><div><font face=3D"Courier =
New"><b>&nbsp; &nbsp; &nbsp; &nbsp; result[i] =3D a[i] * b[i]; // four =
cycles per iteration</b></font></div><div><font face=3D"Courier =
New"><b>}</b></font></div></div><div><br></div><div>If you use SIMD =
libraries (like the previously mentioned Apple=E2=80=99s vDSP bundled =
with MacOS), you likely won=E2=80=99t need to worry about that, but you =
may need to worry about that when accessing things in Lisp. For example, =
<i>typed-aref</i> vectors in LW are always aligned on <b>8 bytes</b> so =
you can=E2=80=99t rely on their lengths alone to identify how many =
objects you are holding there (unless they are a multiple of <b>8</b>, =
or course). You will need to store the object count separately and pass =
it around.</div><br><blockquote type=3D"cite">Lisp access to objects =
with permuted indices will be significantly slower, but that's not the =
point: they exist so that foreign SIMD code can eat them in a format it =
likes.<br></blockquote><div><br></div><div>That is also the case with my =
approach, I estimated that accessing a field in a SIMD aligned struct =
vector requires <b>7</b>&nbsp;cycles (as opposed to just <b>2</b> in =
scalar alignment), here is an example for a <i>single-float</i> custom =
object property access:</div><div><br></div><div>1. Identify the SIMD =
frame in which the object lives (<b>1</b> =
cycle):</div><div><br></div><div><b><font face=3D"Courier New">(setq =
frame-index (ash object-index (- =
simd-frame-size-in-bits)))</font></b></div><div><br></div><div>2. =
Identify the offset of the SIMD frame in the vector (<b>1</b> =
cycle):</div><div><br></div><div><b><font face=3D"Courier New">(setq =
frame-offset (* frame-index =
simd-frame-size))</font></b></div><div><br></div><div>3. Identify the =
position of the object in the frame (<b>1</b> =
cycle):</div><div><br></div><div><font face=3D"Courier New"><b>(setq =
object-in-frame-pos (- object-index =
frame-offset)</b></font></div><div><br></div><div>3. Identify the =
location of the field in its field sequence (<b>1</b> =
cycle):</div><div><br></div><div><b><font face=3D"Courier New">(setq =
field-pos (* object-in-frame-pos =
field-size)</font></b></div><div><br></div><div>4. Identify the actual =
field location in the SIMD frame by skipping over other field sequences =
before it &nbsp;(<b>1</b> cycle):</div><div><br></div><div><b><font =
face=3D"Courier New">(setq field-offset (+ field-sequence-offset =
field-pos))</font></b></div><div><br></div><div>5. Now we can identify =
the actual location of the field in the vector =
&nbsp;(<b>1</b>&nbsp;cycle):</div><div><br></div><div><b><font =
face=3D"Courier New">(setq field-index (+ field-offset =
frame-offset))</font></b></div><div><br></div><div>5. Finally we can get =
our single float field value (<b>1</b> =
cycle):</div><div><br></div><div><b><font face=3D"Courier =
New">(sys:typed-aref =E2=80=99single-float object-vector =
field-index)</font></b></div><div><br></div><div>If all this looks =
fairly complicated, it is because it is, do not worry if it goes over =
your head, it certainly did over mine. This is one of those things that =
you would need to write yourself from scratch to figure =
out.</div><br><blockquote type=3D"cite">So the question is: before I =
redo the macro level (which is going to be fiddly), is something like =
this likely to be helpful for SIMD code (or any other foreign code =
purpose).<br></blockquote><div><br></div><div>Yes, as long as the =
resulting data is aligned in a&nbsp;</div><div><br></div><div><b><font =
face=3D"Courier New">((x=E2=82=80, x=E2=82=81, x=E2=82=82, =
x=E2=82=83)</font></b></div><div><b><font face=3D"Courier =
New">&nbsp;(y</font></b><b><font face=3D"Courier =
New">=E2=82=80</font></b><b><font face=3D"Courier New">, =
y</font></b><b><font face=3D"Courier New">=E2=82=81</font></b><b><font =
face=3D"Courier New">, y</font></b><b><font face=3D"Courier =
New">=E2=82=82</font></b><b><font face=3D"Courier New">, =
y</font></b><b><font face=3D"Courier New">=E2=82=83</font></b><b><font =
face=3D"Courier New">))</font></b></div><div><br></div><div>manner then =
there will be no issues with it on the foreign side. I did find the =
syntax</div><div><br></div><div><b><font face=3D"Courier =
New">...</font></b></div><div><b><font face=3D"Courier New">&nbsp;((x =
:array t :type fixnum :initform 12 :index-mappings ((2 1) (4 2) (6 =
3)))<br>&nbsp;&nbsp;(y :array t :type fixnum :initform 13 =
:index-mappings ((1 4) (3 5) (5 6)))</font></b></div><div><b><font =
face=3D"Courier New">...</font></b></div><div><br></div><div>somewhat =
confusing first, i.e.&nbsp;</div><div><br></div><div><b><font =
face=3D"Courier New">x=E2=82=81 -&gt; v[1], x=E2=82=82 -&gt; v[2], =
x</font></b><b><font face=3D"Courier New">=E2=82=83</font></b><b><font =
face=3D"Courier New">&nbsp;-&gt; v[3]</font></b></div><b><font =
face=3D"Courier New">y=E2=82=80 -&gt; v[4], y=E2=82=81 -&gt; v[5], =
y=E2=82=82</font></b><b><font face=3D"Courier New">&nbsp;-&gt; =
v[6]</font></b><div><br></div><div>wasn=E2=80=99t very obvious first, =
but since it is write once and can also be macro-automated if needed, I =
am not seeing any reasons why it wouldn=E2=80=99t =
work.</div><div><br></div><div><br></div><div><br></div><div><br></div><di=
v><br></div><br><br></body></html>=

--Apple-Mail=_DCD8F2CB-91F0-41DA-B717-7DBE47CE663E--

_______________________________________________
Lisp Hug - the mailing list for LispWorks users
[email protected]
http://www.lispworks.com/support/lisp-hug.html