The Burroughs Scientific Processor was a supercomputer design
that attempted to solve the problem of accessing memory with
stride through a unique memory design.
There were seventeen memory banks, each with its own data bus
and address bus, and consecutive words in memory rotated through
the memory banks.
So if one was retrieving a large block of consecutive words from the
memory, one could always use the full width of the bus on every
cycle (except the last). And if one wanted to access words occurring
at regular intervals - with a stride - as long as the displacement
between consecutive words wasn't a multiple of 17, one could also
always use the full width.
So unlike a conventional memory organization, doing matrix
multiplication on a 2048 by 2048 matrix would be no problem;
2048, as a power of two, isn't relatively prime to 16, but it is
relatively prime to 17.
An easier fix would be to store the matrix in a (2048,2049) matrix -
or (2049,2048) depending on whether your language is row-major
or column-major - but, hey, it was a thought.
Now, I've noted that a 36-bit single-precision float, using the IEEE-754
hidden first bit, could combine the precisiion of IBM 7090 floats with
the exponent range of IBM 360 floats. And that the Control Data
6600 showed that 60-bit double precision is accurate.
But there's no way to generalize the BSP design to handle both of
these widths with the same advantage for stride, unless one accepts
a serious inefficiency (i.e. have 17 memory banks, each with a 12-bit
data bus, store successive words of each float in the same memory
bank).
But if one used different lengths of numbers, maybe something is
possible.
37 and 61 are both prime numbers. So one just has a 37 times 61-bit
data bus chopped both ways!
But that's awfully wide. Can we do better?
38 and 62 are 2*19 and 2*31. The 2 is embedded in both widths,
so it doesn't interfere with avoiding issues - have a memory bank that
looks like 31 38-bit wide memory banks, or 19 62-bit wide memory
banks. That cuts the width almost in half!
We can go for a third the width - 39 is 3 times 13, and 63 is 3 times 21;
21 isn't a prime, but 3 and 7 are also tolerable as values that can be
excluded as factors of efficient strides - it's really only 2 and 5 that
we're worried about. So now the width is 13 times 63, or 21 times 39.
I was hoping for something that would also allow an
intermediate-length float, with a precision of about 10 decimal digits,
to work as well. I could see an arrangement handling intermediates
and doubles: 49 bits and 63 bits, so now we have 9 times 49, or
7 times 63.

John Savard

Quadibloc <jsavard@ecn.ab.ca> writes:
>So unlike a conventional memory organization, doing matrix
>multiplication on a 2048 by 2048 matrix would be no problem;

A correct statement would be: "Like a conventional memory
organization, ...".

Here's a way to implement matrix multiplication that uses stride-1
accesses on the c and b arrays and always accesses the same element of
a in the inner loop:

for (i=0; i<n; i++)
for (k=0; k<m; k++)
for (j=0; j<p; j++)
c[i*p+j]+=a[i*m+k]*b[k*p+j];

Here the two-dimensional array accesses are explicitly converted into
one-dimensional accesses to make it more obvious what is going on.

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

Anton Ertl wrote:
> Quadibloc <jsavard@ecn.ab.ca> writes:
>> So unlike a conventional memory organization, doing matrix
>> multiplication on a 2048 by 2048 matrix would be no problem;
>
> A correct statement would be: "Like a conventional memory
> organization, ...".
>
> Here's a way to implement matrix multiplication that uses stride-1
> accesses on the c and b arrays and always accesses the same element of
> a in the inner loop:
>
> for (i=0; i<n; i++)
> for (k=0; k<m; k++)
> for (j=0; j<p; j++)
> c[i*p+j]+=a[i*m+k]*b[k*p+j];
>
> Here the two-dimensional array accesses are explicitly converted into
> one-dimensional accesses to make it more obvious what is going on.

That is a nice way to organize it.

A compiler would probably convert/optimize it to something similar to this?

ai = a
ci = c
for (i=0; i<n; i++) {
bk = b
for (k=0; k<m; k++) {
aik = ai[k];
for (j=0; j<p; j++) {
ci[j] += aik * bk[j]
}
bk += p;
}
ai += m;
ci += p;
}

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Thursday, August 24, 2023 at 11:04:35 AM UTC+3, Anton Ertl wrote:
> Quadibloc <jsa...@ecn.ab.ca> writes:
> >So unlike a conventional memory organization, doing matrix
> >multiplication on a 2048 by 2048 matrix would be no problem;
> A correct statement would be: "Like a conventional memory
> organization, ...".
>
> Here's a way to implement matrix multiplication that uses stride-1
> accesses on the c and b arrays and always accesses the same element of
> a in the inner loop:
>
> for (i=0; i<n; i++)
> for (k=0; k<m; k++)
> for (j=0; j<p; j++)
> c[i*p+j]+=a[i*m+k]*b[k*p+j];
>
> Here the two-dimensional array accesses are explicitly converted into
> one-dimensional accesses to make it more obvious what is going on.
>
> - anton
> --
> 'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
> Mitch Alsup, <c17fcd89-f024-40e7...@googlegroups.com>

Daxpy-based matrix multiplication is simple and performance is relatively
robust, but it is not likely to reach top performance on any of modern
general-purpose CPUs. There are too many stores in the inner loop :(
Last time I tried it on Skylake(C), even more sophisticated variations of
this idea, like updating 2 lines of c[] by 3-4 lines of b[] in the common
inner loop were 20-30% slower than the best schemes.

For non-small matrices all best schemes were variations of multiple dot
product - m rows of a[] multiplied by n columns of b[] in the inner loop;
where n is typically 4 or 5 multiplied by number of lanes in SIMD register,
i.e. n=16 or 20 for DGEMM on AVX2.
These schemes are relatively fragile due to two factors:
- they do contain non-unit-stride access to b[] in the inner loop
- one has to have enough of accumulators to hide the latency of FMA.
Despite fragility with proper tuning they are fastest.
There are two main approaches for making non-unit-stride access to be
bearable:
1. Create compact copy of n columns of b[] and reuse it over many
rows of a[] and c[].
2. Divide a[] into vertical stripes, divide b[] into horizontal stripes.
Select the size of the stripe Ks such that n*Ks easily fit in L1D cache.

The 1st approach is best when all dimensions of a[] and b[] are big.
The 2nd approach can be preferable for relatively smaller matrices,
but it is sensitive to physical width of matrix b[]. It would fail badly
when width of b[] is an exact multiple of size of L1D way.

I'd imagine that on something like POWER 7/8/9/10, because of write through
policy of their L1D cache, the disadvantage of Daxpy vs dot is even bigger
than on Skylake(C).

Michael S <already5chosen@yahoo.com> writes:
>On Thursday, August 24, 2023 at 11:04:35=E2=80=AFAM UTC+3, Anton Ertl wrote=
>:
>> Quadibloc <jsa...@ecn.ab.ca> writes:=20
>> >So unlike a conventional memory organization, doing matrix=20
>> >multiplication on a 2048 by 2048 matrix would be no problem;
>> A correct statement would be: "Like a conventional memory=20
>> organization, ...".=20
....
>Daxpy-based matrix multiplication is simple and performance is relatively
>robust, but it is not likely to reach top performance on any of modern
>general-purpose CPUs.

My point was that you don't need funny memory organizations in order
to get good matrix multiplication performance.

"DAXPY-based" matrix multiplication shows that in four lines of code.

If I use five lines of code to implement a "DDOT-based" matrix
multiplication (in both cases I call neither DAXPY nor DDOT), I fall
prey to the limitations of the memory subsystem:

for (i=0; i<n; i++)
for (j=0; j<p; j++) {
for (k=0, r=0.0; k<m; k++)
r += a[i*m+k]*b[k*p+j];
c[i*p+j]=r;
}

On an Ivy Bridge this takes 23.3 cycles per inner-loop iteration with
n,p,m=700 (probably similar for n,p,m=2048 without transparent huge
pages. THP can be hit or miss; when they hit, it's 4.4 cycles per
iteration.

By contrast, the "DAXPY-based" variant I showed earlier takes
1.4cycles per original iteration with auto-vectorization from gcc -O3.

OpenBLAS with one thread takes 0.36cycles/original iteration, so yes,
there is room for improvement beyond the 4-line variant shown above,
but OpenBLAS is much bigger and more complex than the 4-line variant.

0.36 cycles/original iteration means 11.11 DP FLOPS/cycle.
<https://en.wikichip.org/wiki/flops> reports that Ivy Bridge can
perform no more than 8 DP FLOPS per cycle, so this makes me wonder
what's wrong.

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

On Thursday, August 24, 2023 at 3:12:33 PM UTC+3, Anton Ertl wrote:
> Michael S <already...@yahoo.com> writes:
> >On Thursday, August 24, 2023 at 11:04:35=E2=80=AFAM UTC+3, Anton Ertl wrote=
> >:
> >> Quadibloc <jsa...@ecn.ab.ca> writes:=20
> >> >So unlike a conventional memory organization, doing matrix=20
> >> >multiplication on a 2048 by 2048 matrix would be no problem;
> >> A correct statement would be: "Like a conventional memory=20
> >> organization, ...".=20
> ...
> >Daxpy-based matrix multiplication is simple and performance is relatively
> >robust, but it is not likely to reach top performance on any of modern
> >general-purpose CPUs.
> My point was that you don't need funny memory organizations in order
> to get good matrix multiplication performance.
>
> "DAXPY-based" matrix multiplication shows that in four lines of code.
>
> If I use five lines of code to implement a "DDOT-based" matrix
> multiplication (in both cases I call neither DAXPY nor DDOT), I fall
> prey to the limitations of the memory subsystem:
> for (i=0; i<n; i++)
> for (j=0; j<p; j++) {
> for (k=0, r=0.0; k<m; k++)
> r += a[i*m+k]*b[k*p+j];
> c[i*p+j]=r;
> }
>
> On an Ivy Bridge this takes 23.3 cycles per inner-loop iteration with
> n,p,m=700 (probably similar for n,p,m=2048 without transparent huge
> pages. THP can be hit or miss; when they hit, it's 4.4 cycles per
> iteration.
>
> By contrast, the "DAXPY-based" variant I showed earlier takes
> 1.4cycles per original iteration with auto-vectorization from gcc -O3.
>
> OpenBLAS with one thread takes 0.36cycles/original iteration, so yes,
> there is room for improvement beyond the 4-line variant shown above,
> but OpenBLAS is much bigger and more complex than the 4-line variant.
>
> 0.36 cycles/original iteration means 11.11 DP FLOPS/cycle.
> <https://en.wikichip.org/wiki/flops> reports that Ivy Bridge can
> perform no more than 8 DP FLOPS per cycle, so this makes me wonder
> what's wrong.

2 FLOP / 0.36 = 5.56 DP FLOPS/cycle. Rmax/Rpeak = 69%
I didn't try, but would think that for n,p,m=700 I can do a little better,
may be Rmax/Rpeak = 75%.
Of course, it's much easier to optimize for specific matrix size than to write
DGEMM that is not too bad at all sizes and shapes.

To your bigger point, I fully agree that achieving reasonable performance by
DAXPY approach is much easier. Just saying that it does not cut for very
best. And what does cut for best is somewhat sensitive to stride aliasing.
Not horribly so, but you can't ignore the issue completely.

> - anton
> --
> 'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
> Mitch Alsup, <c17fcd89-f024-40e7...@googlegroups.com>

Michael S <already5chosen@yahoo.com> schrieb:

> Of course, it's much easier to optimize for specific matrix size than to write
> DGEMM that is not too bad at all sizes and shapes.

For small sizes (2,3,4,5) a straightfoward loop is usually better,
becaue the highly-tuned BLAS libraries have a big startup time.
This is especially true if the size and strides are known at
compile time.

Example:

subroutine mm(a,b,c)
double precision, dimension(2,2) :: a, b, c
c = matmul(a,b)
end subroutine mm

will be translated, on my Rocket Lake home system, with -Ofast
-march=native, to

mm_:
..LFB0:
.cfi_startproc
vmovupd (%rdi), %ymm0
vmovupd (%rsi), %ymm1
vpermpd $68, %ymm0, %ymm3
vpermilpd $0, %ymm1, %ymm2
vpermpd $238, %ymm0, %ymm0
vpermilpd $15, %ymm1, %ymm1
vmulpd %ymm1, %ymm0, %ymm0
vfmadd231pd %ymm3, %ymm2, %ymm0
vmovupd %ymm0, (%rdx)
vzeroupper
ret

(I am confident that you or Terje will point out numerous glaring
inefficiencies in the compiler-generated code :-)

On Thursday, August 24, 2023 at 6:22:28 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > Of course, it's much easier to optimize for specific matrix size than to write
> > DGEMM that is not too bad at all sizes and shapes.
> For small sizes (2,3,4,5) a straightfoward loop is usually better,
> becaue the highly-tuned BLAS libraries have a big startup time.
> This is especially true if the size and strides are known at
> compile time.
>
> Example:
>
> subroutine mm(a,b,c)
> double precision, dimension(2,2) :: a, b, c
> c = matmul(a,b)
> end subroutine mm
>
> will be translated, on my Rocket Lake home system, with -Ofast
> -march=native, to
>
> mm_:
> .LFB0:
> .cfi_startproc
> vmovupd (%rdi), %ymm0
> vmovupd (%rsi), %ymm1
> vpermpd $68, %ymm0, %ymm3
> vpermilpd $0, %ymm1, %ymm2
> vpermpd $238, %ymm0, %ymm0
> vpermilpd $15, %ymm1, %ymm1
> vmulpd %ymm1, %ymm0, %ymm0
> vfmadd231pd %ymm3, %ymm2, %ymm0
> vmovupd %ymm0, (%rdx)
> vzeroupper
> ret
>
> (I am confident that you or Terje will point out numerous glaring
> inefficiencies in the compiler-generated code :-)

It is actually not bad at all.
I think that I can improve it a little by replacing
vmovupd (%rdi), %ymm0
vpermpd $68, %ymm0, %ymm3
vpermpd $238, %ymm0, %ymm0
by
vbroadcastf128 (%rdi), %ymm3
vbroadcastf128 16(%rdi), %ymm0

But it would be very close and in some less common situations
my variant could end up slower.

Also instead of
vpermilpd $0, %ymm1, %ymm2
vpermilpd $15, %ymm1, %ymm1
I'd write
vunpcklpd %ymm1,%ymm1,%ymm2
vunpckhpd %ymm1,%ymm1,%ymm1
It is not faster, but 4 bytes shorter, so can be of advantage
when function is used rarely.

But good work by compiler nevertheless.

Michael S wrote:
> On Thursday, August 24, 2023 at 6:22:28 PM UTC+3, Thomas Koenig wrote:
>> Michael S <already...@yahoo.com> schrieb:
>>> Of course, it's much easier to optimize for specific matrix size than to write
>>> DGEMM that is not too bad at all sizes and shapes.
>> For small sizes (2,3,4,5) a straightfoward loop is usually better,
>> becaue the highly-tuned BLAS libraries have a big startup time.
>> This is especially true if the size and strides are known at
>> compile time.
>>
>> Example:
>>
>> subroutine mm(a,b,c)
>> double precision, dimension(2,2) :: a, b, c
>> c = matmul(a,b)
>> end subroutine mm
>>
>> will be translated, on my Rocket Lake home system, with -Ofast
>> -march=native, to
>>
>> mm_:
>> .LFB0:
>> .cfi_startproc
>> vmovupd (%rdi), %ymm0
>> vmovupd (%rsi), %ymm1
>> vpermpd $68, %ymm0, %ymm3
>> vpermilpd $0, %ymm1, %ymm2
>> vpermpd $238, %ymm0, %ymm0
>> vpermilpd $15, %ymm1, %ymm1
>> vmulpd %ymm1, %ymm0, %ymm0
>> vfmadd231pd %ymm3, %ymm2, %ymm0
>> vmovupd %ymm0, (%rdx)
>> vzeroupper
>> ret
>>
>> (I am confident that you or Terje will point out numerous glaring
>> inefficiencies in the compiler-generated code :-)
>
>
> It is actually not bad at all.
> I think that I can improve it a little by replacing
> vmovupd (%rdi), %ymm0
> vpermpd $68, %ymm0, %ymm3
> vpermpd $238, %ymm0, %ymm0
> by
> vbroadcastf128 (%rdi), %ymm3
> vbroadcastf128 16(%rdi), %ymm0
>
> But it would be very close and in some less common situations
> my variant could end up slower.
>
> Also instead of
> vpermilpd $0, %ymm1, %ymm2
> vpermilpd $15, %ymm1, %ymm1
> I'd write
> vunpcklpd %ymm1,%ymm1,%ymm2
> vunpckhpd %ymm1,%ymm1,%ymm1
> It is not faster, but 4 bytes shorter, so can be of advantage
> when function is used rarely.
>
> But good work by compiler nevertheless.
>
I agree. Any time a compiler gets within 80-90% of my best hand-tuned
effort, I stop worrying about that particular algorithm and simple move
on to more useful challenges, leaving the compiler to do what it does so
well.(1*)

On a few (rare, but they do happen!) occations, I have even learned new
approaches from looking at compiler output.

1) The only exception is for progrmaming competitions where every single
cycle matters. For real code, even 50% of speed of light/theoretical max
is usually fine.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Thursday, August 24, 2023 at 4:00:35 AM UTC-6, Michael S wrote:

> Daxpy-based matrix multiplication is simple and performance is relatively
> robust, but it is not likely to reach top performance on any of modern
> general-purpose CPUs. There are too many stores in the inner loop :(

True, but my proposed memory architecture also doesn't work on modern
CPUs, because it depends on accessing individual array elements from
memory, so it's mismatched to DRAM. In any event, on further reflection,
the idea proposed here is simply too elaborate to be useful. If a derivative
of the BSP is to be used, it will have to be a much simpler one.

Basically, this kind of memory structure would need to be applied to
static RAM inside a module, using a technique such as HBM.

My first thought of a simplification was to allow some memory wastage.

19*5 = 95, 23*4 = 92, and 31*3 = 93, so one could have a memory array
with 95 12-bit data lanes as its width that would have good non-interfering
sizes for double (multiples of 19 excluded), intermediate (multiples of 23
excluded), and single (multiples of 31 excluded).

However, it makes more sense to just offer 60 lanes, each 12-bits wide, and
simply add one to the row (or column, if the language is row-major, as
FORTRAN usually is) size when it's a round number that would have problems.

John Savard

On Friday, August 25, 2023 at 8:28:10 AM UTC-5, Quadibloc wrote:
> On Thursday, August 24, 2023 at 4:00:35 AM UTC-6, Michael S wrote:
>
> > Daxpy-based matrix multiplication is simple and performance is relatively
> > robust, but it is not likely to reach top performance on any of modern
> > general-purpose CPUs. There are too many stores in the inner loop :(
<
> True, but my proposed memory architecture also doesn't work on modern
> CPUs, because it depends on accessing individual array elements from
> memory, so it's mismatched to DRAM.
<
And especially mismatched to Cache line sizes, sets, and sizes that make
DRAM accesses unpredictable. Also note:: BSP was inherently 16 CPUS
making 1 data access each per cycle, whereas a modern chip with 16
CPUs is making at most 1 line access every 4-8 beats.
<
> In any event, on further reflection,
> the idea proposed here is simply too elaborate to be useful. If a derivative
> of the BSP is to be used, it will have to be a much simpler one.
>
> Basically, this kind of memory structure would need to be applied to
> static RAM inside a module, using a technique such as HBM.
>
> My first thought of a simplification was to allow some memory wastage.
>
> 19*5 = 95, 23*4 = 92, and 31*3 = 93, so one could have a memory array
> with 95 12-bit data lanes as its width that would have good non-interfering
> sizes for double (multiples of 19 excluded), intermediate (multiples of 23
> excluded), and single (multiples of 31 excluded).
>
> However, it makes more sense to just offer 60 lanes, each 12-bits wide, and
> simply add one to the row (or column, if the language is row-major, as
> FORTRAN usually is) size when it's a round number that would have problems.
>
> John Savard