This package provides floating-point matrix multiplication. It runs on Forth-94 and Forth-2012 systems with the floating-point wordset.

include matmul.4th

to get the word

matmulr ( a b c n1 n2 n3 -- )

MATMULR computes C=AxB, where A, B, and C are matrixes: where A has n1 rows and n2 columns, B has n2 rows and n3 columns, and C has n1 rows and n3 columns. Matrixes are stored in row-major order, i.e., the elements of one row are adjacent, whereas the items of a column are separated by all the items in the first row.

Testing

To test whether this works as intended with your Forth system, type

create matmul-verbose include matmul.4th include test.4th

Performance

This may not be the best-performing matrix multiplication written in Forth and certainly not outside of Forth, but it has the advantage of just working (at least that's the intention). It's also not too badly performing: Gforth performs 0.8 FLOP/cycle on a Core i5-6600K using the faxpy primitive, and ideally Forth compilers will generate as good code for FAXPY-NOSTRIDE as that for the FAXPY primitive.

Here are some results (using unrolled versions):

cycles Core i5-6600K
Forth-2012 Forth-94 Forth-94 primitive
          Rees/Ertl   pahihu     faxpy
 *1048M       1202M    1517M           vfxlin Ndp387.fth
  2484M       2945M    2764M     *316M gforth-fast
 *3584M       5067M    5017M           SwiftForth HW stack, no WAIT

Details (irrelevant for most usage)

MATMULR calls a word FAXPY-NOSTRIDE, and the library contains several implementations of this word, and automatically selects the one expected to give the best performance. If you want to know which implementation is used, you can CREATE MATMUL-VERBOSE before including matmul.4th. You can also determine which implementation is used by defining a constant FAXPY-NOSTRIDE-VARIANT. The values for this constant are as follows:

unrolling
no yes
 1     use FAXPY (a primitive in Gforth)
 4  2  separate-stack (Forth-2012) implementation
 5  3  Forth-94 implementation by Joel Rees and Anton Ertl
15 13  Forth-94 implementation by pahihu

The Forth-94 implementations work on systems with separate FP stacks as well as on systems with a combined FP and data stack with any number of cells per FP item.

Use as Benchmark

For use as benchmark, there is benchmark.4th, which performs one 500x500 matrix multiplication. There is also the script runbenchs; by default, runbenchs will run the benchmark on gforth-fast, vfxlin (with various FP packages), and sf if available (and for VFX, if %vfxpath% works), and it will use perf (with a fallback to time) for timing. You can influence these choices and more by setting environment variables before calling runbenchs. E.g., if you want to benchmark iforth only, you can do (on a bourne-compatible shell):

SF=iforth GFORTH=xx VFXLIN=xx ./runbenchs

The environment variables used are:

TIMING  command used for measuring, e.g., `time`
RUN     file used for running benchmark (RUN=test.4th for testing)
GFORTH  how to invoke Gforth (default: gforth-fast)
SF      how to invoke SwiftForth (default: sf)
VFXLIN  How to invoke VFX (default: vfxlin)
VFXHOME where the VfxLin directory is (default: %vfxpath%)
VFXLIB  where the FP packages can be found (default: $VFXHOME/Lib)
VFXX86  where the x86-specific FP packages can be found

The defaults are reasonable, but you may need to define VFXHOME yourself.

You can also use this script to test all the different variants on several systems with

TIMING=command RUN=test.4th ./runbenchs