On Wednesday, July 21, 2021 at 6:34:03 PM UTC-5, JimBrakefield wrote:
> On Wednesday, July 21, 2021 at 5:44:54 PM UTC-5, MitchAlsup wrote:
> > On Wednesday, July 21, 2021 at 4:48:32 PM UTC-5, luke.l...@gmail.com wrote:
> > > On Wednesday, July 21, 2021 at 9:20:21 PM UTC+1, Thomas Koenig wrote:
> > >
> > > > I agree that these instruction sets are badly designed and a major
> > > > headache for compiler writers (anybody who doubts that can read
> > > > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53947 and weep.)
> > > dang - 200+ dependent bugs.
> > > > For POWER, which you work on, there is also VFX.
> > > you may not be aware, this video helps explain it,
> > > https://youtu.be/gexy0z1YqFY?t=673
> > >
> > > i have made it very bluntly plain and clear that i am flatly refusing
> > > to implement any part of Power ISA SIMD in Libre-SOC. i consider it to
> > > be so harmful a paradigm that even though EABIv2.0 makes it mandatory
> > > i've still said "no f***** way". it's LITERALLY 700 extra instructions.
> > > Power ISA's scalar instruction count is just over 200 (214 or so).
> > > [which part of "RISC" was not clear, when VSX was added??]
> > >
> > > this will leave our processor "behind" for a number of years whilst the
> > > EABIv2.0 is sorted out and SIMD being mandatory is removed, replaced
> > > with something that makes it optional again (but at the same time makes
> > > it possible to run SIMD binaries). despite this, *i still said absolutely no
> > > way* - SIMD is that bad - and we will instead use illegal instruction
> > > trap-and-emulate.
> > >
> > > to give you some idea, Paul Mackerras has added only some of the VSX
> > > instructions to Microwatt, and it's increased the LUT4 usage in an FPGA
> > > from 20,000 to *50,000*.
> > > > Again, by far not
> > > > the best of design, but they also do something.
> > > yes - cause no end of problems.
> > > > All these instructions work well enough that, for example, you get
> > > > to 60-70% of the theoretical maximum floating point performance
> > > > for matrix multiplication
> > > SVP64's Matrix REMAP Schedule will get 100%.
> > > > if you take a 30-year old Fortran code
> > > > from Netlib, mechanically translate it to C via f2c, tune for more
> > > > recent cache sizes, sprinkle a few restrict qualifiers here and
> > > > there and translate for several versions for AVX 128, AVX 256 and
> > > > AVX 512 and select at runtime.
> > > no assembly optimisations are needed for SVP64 Matrix-Multiply.
> > > it's too short (and specifically designed for the job).
> > > > The original code is at
> > > > http://www.netlib.org/blas/gemm_based/ssgemmbased.tgz
> > > > and what I am referring to is at
> > > > https://gcc.gnu.org/git/?p=gcc.git;a=blob_plain;f=libgfortran/generated/matmul_r8.c;hb=refs/heads/trunk
> > > > .
> > > >
> > > > Now, that is something to beat.
> > > done, last week. replaced by three instructions.
> > >
> > > yes, that wasn't a typo: three (3) instructions.
> > >
> > > https://git.libre-soc.org/?p=openpower-isa.git;a=blob;f=src/openpower/decoder/isa/test_caller_svp64_matrix.py;hb=HEAD
> > >
> > > there are some caveats:
> > >
> > > 1) the size of the register file (Libre-SOC will be going for 128 FPs and 128 INTs)
> > > is a limiting factor.
> > > 2) the total number of FMACs needed must not exceed 127.
> > > 3) if you need the result to be zero'd out before beginning, it's 5 instructions not 3.
> > <
> > How many instructions get executed in a 10,000×10,000 matrix multiplication ?
> > In double precision. How many memory references are performed ?
> > >
> > > if you want to experiment and see what's going on, here's some stand-alone
> > > executable code:
> > > https://git.libre-soc.org/?p=openpower-isa.git;a=blob;f=src/openpower/decoder/isa/remapyield.py;hb=HEAD
> > >
> > > here's a unit test for it:
> > > https://git.libre-soc.org/?p=libreriscv.git;a=blob;f=openpower/sv/remapmatrixprint.py;hb=HEAD
> > >
> > > what that does is, it can be used to yield any of the three indices x y and z
> > > needed for performing matrix element lookups (assuming the data has
> > > been flattened from 2D to 1D).
> > >
> > > it should also be possible to specify that the data has been transposed
> > > (hmm, must add those as modes in the setup instruction, svshape)
> > > which is simply a matter of turning the sequence-order round of which
> > > of the three loops shall be inner and which outer.
> > >
> > > it should also be possible to set "inversion" (mirroring) mode, such that
> > > the input matrices may be considered mirrored (in any given dimension)
> > > but *without* actually needing to copy it.
> > >
> > > l.
> This is a trick question: It's a question of how to make best use of the caches.
> Or how do you manage cache churn?
> |> How many instructions get executed in a 10,000×10,000 matrix multiplication ?
> |> In double precision. How many memory references are performed ?
<
I think the real question is "Why not just use *GEMM() from BLAS level 3 ??
<
Knowing this is for scale-rotate,projections of 3D graphics is justification
for doing something other than DGEMM, but random matrix multiplication
one real arrays is not (that is use DGEMM--it will surprise you at how efficient
it is.)

On Wednesday, July 21, 2021 at 9:20:21 PM UTC+1, Thomas Koenig wrote:

> The original code is at
> http://www.netlib.org/blas/gemm_based/ssgemmbased.tgz
> and what I am referring to is at
> https://gcc.gnu.org/git/?p=gcc.git;a=blob_plain;f=libgfortran/generated/matmul_r8.c;hb=refs/heads/trunk

Thomas: apologies, small screen, I originally looked at the
function header, and went, "naah, it'll be some tiny matmul function"
therefore comparative to SVP64 MM REMAP.

naah :)

Mitch later points out that the smaller matrices for 3D work are
quite different workloads. these are primarily lots of repeated
LOAD PROCESS STORE, nothing common between them, and regfiles
have to be made big enough to prevent spill (128, 256).

Much of what I sm doing is targetted at that style of workload
(8x8 DCT/IMDCT), no register spill.

I will have to look again at the larger DCT though because if it
is quite common then it is worth investigating. although, if
the better algorithm to use is RDFT i have that covered as well
in a different way.

l.

luke.l...@gmail.com <luke.leighton@gmail.com> schrieb:
> On Wednesday, July 21, 2021 at 9:20:21 PM UTC+1, Thomas Koenig wrote:

>> All these instructions work well enough that, for example, you get
>> to 60-70% of the theoretical maximum floating point performance
>> for matrix multiplication
>
> SVP64's Matrix REMAP Schedule will get 100%.
>
>> if you take a 30-year old Fortran code
>> from Netlib, mechanically translate it to C via f2c, tune for more
>> recent cache sizes, sprinkle a few restrict qualifiers here and
>> there and translate for several versions for AVX 128, AVX 256 and
>> AVX 512 and select at runtime.
>
> no assembly optimisations are needed for SVP64 Matrix-Multiply.
> it's too short (and specifically designed for the job).
>
>> The original code is at
>> http://www.netlib.org/blas/gemm_based/ssgemmbased.tgz
>> and what I am referring to is at
>> https://gcc.gnu.org/git/?p=gcc.git;a=blob_plain;f=libgfortran/generated/matmul_r8.c;hb=refs/heads/trunk
>> .
>>
>> Now, that is something to beat.
>
> done, last week. replaced by three instructions.

Do your three instructions also reorder loops, interchange strides
and copy stuff around in to a buffer? :-)

If you take a look at the source, this is where you'll find the
performance mostly comes from.

A simple looping construct, even fully vectorized, will not even
start come close for large matrices. It's how caching is dealt
with which allows whatever vectorization / SIMD / ... happens in
the central loop(s) to be used efficiently.

Or, if I may quote Terje, "Almost all programming can be viewed
as an exercise in caching."

> yes, that wasn't a typo: three (3) instructions.
>
> https://git.libre-soc.org/?p=openpower-isa.git;a=blob;f=src/openpower/decoder/isa/test_caller_svp64_matrix.py;hb=HEAD
>
> there are some caveats:
>
> 1) the size of the register file (Libre-SOC will be going for 128 FPs and 128 INTs)
> is a limiting factor.

Ah.

I thought that the aim of that project was to stay compatible with
POWER, to make use of the existing infrastructure, but looking at
the web site, that does not appear to be the case.

So, if you want to move away from POWER anyway, what is the reason
to start with it? I can see why you would not want to use RISC-V
(whose drawbacks we recently discussed :-), but what about Mitch's
My 66000? It is certainly a much cleaner design than POWER,
which has grown over x generations.

[...]

> it should also be possible to specify that the data has been transposed
> (hmm, must add those as modes in the setup instruction, svshape)
> which is simply a matter of turning the sequence-order round of which
> of the three loops shall be inner and which outer.
>
> it should also be possible to set "inversion" (mirroring) mode, such that
> the input matrices may be considered mirrored (in any given dimension)
> but *without* actually needing to copy it.

For large matrices, that is of questionable use because you need
to copy into cache-friendly buffers anyway.

For matrices much smaller than relevant for caches (3*3 or 4*4
matrices are also often used) that can of course be a win.

Thomas Koenig wrote:
> luke.l...@gmail.com <luke.leighton@gmail.com> schrieb:
>> On Wednesday, July 21, 2021 at 7:20:29 PM UTC+1, MitchAlsup wrote:
>>> On Wednesday, July 21, 2021 at 11:01:05 AM UTC-5, luke.l...@gmail.com wrote:
>>>>> short vector SIMD does not count ??
>>>> i've not heard of that one. where can i find out about it?
>>> MMX
>>> SSE
>>> SSE2
>>> AVX
>>
>> oh right, those. yes. god-awful messes, all of them. no, they don't count:
>> they're distinctly part of the problem, not the solution.
>
> I agree that these instruction sets are badly designed and a major
> headache for compiler writers (anybody who doubts that can read
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53947 and weep.)

I like the status ("verified") and the "Depends on:" list. :-)
>
> For POWER, which you work on, there is also VFX. Again, by far not
> the best of design, but they also do something.
>
> All these instructions work well enough that, for example, you get
> to 60-70% of the theoretical maximum floating point performance

My rule of thumb for optimization work has been "is this within a factor
of two of speed of light?", i.e. theoretical maximum?

If it is, I'll usually go on to look for other opportunities unless I've
already done so and this is the final iteration of a competition.

> for matrix multiplication if you take a 30-year old Fortran code
> from Netlib, mechanically translate it to C via f2c, tune for more
> recent cache sizes, sprinkle a few restrict qualifiers here and
> there and translate for several versions for AVX 128, AVX 256 and
> AVX 512 and select at runtime.

That's usually the way to the fastest code, I like to use runtime
version selection by benchmarking the various applicable versions and
picking the fastest on the current HW.

Terje

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

On Thursday, July 22, 2021 at 6:17:22 AM UTC+1, Thomas Koenig wrote:

> > done, last week. replaced by three instructions.
> Do your three instructions also reorder loops,

yes, that's exactly what they're designed to do. it's termed "REMAP".

> interchange strides

yes.

> and copy stuff around in to a buffer? :-)

that's an "action". so you would apply loop-reordering (REMAP)
*to* a MV or a LD or a ST.

> If you take a look at the source, this is where you'll find the
> performance mostly comes from.
>
> A simple looping construct, even fully vectorized, will not even
> start come close for large matrices.

no, correct: i didn't read the source in detail (small screen),
agreed 100%.

however, we're not _doing_ a "simple" looping construct....
but, i have to focus unfortunately: i'd love to go down the
additional rabbithole of coping with 10k+ matrix dimensions,
i'll have to settle for 3D workloads (3, 4, 5. 8 dimensions) at
the moment.

> Or, if I may quote Terje, "Almost all programming can be viewed
> as an exercise in caching."

:)

> > yes, that wasn't a typo: three (3) instructions.
> >
> > https://git.libre-soc.org/?p=openpower-isa.git;a=blob;f=src/openpower/decoder/isa/test_caller_svp64_matrix.py;hb=HEAD
> >
> > there are some caveats:
> >
> > 1) the size of the register file (Libre-SOC will be going for 128 FPs and 128 INTs)
> > is a limiting factor.
> Ah.
>
> I thought that the aim of that project was to stay compatible with
> POWER,

with the newly-created "Embedded Floating Point Compliancy Subset",
yes. what we will then do is take advantage of the fact that it is OPTIONAL
to add say an MMU, and anything else.

Henriok and i worked to update the wikipedia page a couple months back:
https://en.wikipedia.org/wiki/Power_ISA#Compliancy

> to make use of the existing infrastructure, but looking at
> the web site, that does not appear to be the case.

at the low level (basic Scalar level), yes. we can execute basic non-accelerated
programs.

that is good enough to get started.

however a SIMD ISA as a front-end paradigm is just hopelessly inadequate
for 3D VPU GPU work.

it would be commercially suicidal to even attempt to use SIMD (designed
around 2003 and great for the PS3), in the year 2021, for today's GPU and
VPU workloads.

> So, if you want to move away from POWER anyway, what is the reason
> to start with it?

because of the existing infrastructure.

> I can see why you would not want to use RISC-V
> (whose drawbacks we recently discussed :-),

this says everything that needs to be said... technically:
https://news.ycombinator.com/item?id=24459314

> but what about Mitch's
> My 66000? It is certainly a much cleaner design than POWER,
> which has grown over x generations.

llvm gcc binutils linux kernel u-boot libc6 OSes packages libraries mindshare
zephyr RTOS openocd

add up the cost in time alone of developing ports of those and you start
to see why OpenRISC 1k took around 12 years to reach maturity.

> For large matrices, that is of questionable use because you need
> to copy into cache-friendly buffers anyway.
>
> For matrices much smaller than relevant for caches (3*3 or 4*4
> matrices are also often used) that can of course be a win.

yes, that's the plan, that's what SVP64 REMAP is primarily designed for.

l.

On Thursday, July 22, 2021 at 12:34:03 AM UTC+1, JimBrakefield wrote:

> This is a trick question: It's a question of how to make best use of the caches.
> Or how do you manage cache churn?

fascinatingly, then, the algorithm i developed err... 30 years ago would still be relevant
if adapted to L1/L2 caches rather than the original exercise, Real memory plus much
slower external storage.

the last row (middle loop) can be kept in memory before moving on to the next column.

for i in 0 to n
start, end, direction = end, start, -direction
for j in start to end by direction
for k in 0 to p
madds

except, err, offset the middle access by one so that the values loaded from the middle row
can be left in memory for one inner loop.

efficient arge matrix multiply is a whooole rabbithole area of computer science research :)

l.

now implemented, it's 5 instructions to perform the inner and outer
butterfly schedules:

"svremap 27, 1, 0, 2, 0, 1, 1",
"svshape 8, 1, 1, 2, 0",
"sv.fdmadds 0.v, 0.v, 0.v, 8.v",
"svshape 8, 1, 1, 3, 0",
"sv.fadds 0.v, 0.v, 0.v"

the first instruction should technically speaking come after the
second. svshape 8,1,1,2,0 establishes a DCT of length 8,
where "2" indicates "outer butterfly mode". this sets up
the REMAP schedule but does not say which actual registers
the REMAPing should apply to: that's svremap's job.

in binary, 27 - "0b11011" - applies to each of the 5 arguments,
which represent 1st src, 2nd src, 3rd src, 1st dest and 2nd dest
respectively, saying which register should be REMAPed to which
schedule (0-3)

sv.fdmadds is the core in-place (twin) MUL-ADD-SUB, it does this:

FRS <- FPADD32(FRA, FRB)
tmp <- FPSUB32(FRA, FRB)
FRT <- FPMUL32(FRC, tmp)

so it is *three* input registers and *two* outputs - this is key
to achieving an in-place DCT because the registers which are
going to be overwritten are "in-flight" in the same pipeline
stage.

fourth instruction is the "inner" butterfly of the Lee DCT, and
this performs the usual iterative shifted-by-one adds. nothing
special, here: simple register RaW / WaR hazard avoidance
can ensure that the chain of adds operate correctly:
ADD x[i], x[i-1]

where it will get very interesting is to do in-place DCT without
having the COS tables pre-computed, using Vertical-First
Mode (similar to VVM)

l.