Donald Knuth's 64-bit MMIX architecture includes several novel instructions that operate matrixwise on an 8x8 square matrix (MOR, MXOR).

(MMIX also has instructions like BDIF that operate vectorwise on 8x 8-bit integers, which is less novel; x86 does that just fine.)

Treating your 64-bit word as a square matrix is possible only when the word size is a perfect square.

Back in the day, there were computers with a 36-bit word length. That's famously enough to store 6x 6-bit (pre-ASCII) characters in a single word. Did any historical 36-bit computer ever offer operations that were best understood by reading machine words as 6x6 matrices?

How about any historical 16-bit computer?

Besides complicated instructions like MMIX's MOR and MXOR, the simplest example of a "matrix operation" in the sense I mean would be a unary TRANSPOSE operation:

TRANSPOSE(1001'1010'1010'1011) = 1111'0000'0111'1001
TRANSPOSE(1001'0010'1011'0100) = 1010'0001'0110'1010
  • Does MAC instruction count? These simplify matrix-matrix or matrix-vector multiplication. In general, the trade-off for having a very specific "6x6" etc. instruction against using general indexing or field extraction (which many machines did have) is just to bad. You want instructions you can use everywhere, not instructions with a single purpose. – dirkt Sep 27 at 4:00
  • And if you look at how e.g. the Cray-1 handles vector operations (which includes matrix operations): It does so by having complete vector unit with vector registers, which operating on floating point, because that's where you needed computing power. 6-bit integers are pretty worthless for the applications where you wanted to use a computer back then. – dirkt Sep 27 at 4:02
  • Some Crays also had a powerful BMM instruction for bit matrix multiplication. See: Yedidya Hilewitz, Cedric Lauradoux and Ruby B. Lee, "Bit Matrix Multiplication in Commodity Processors". In: 2008 International Conference on Application-Specific Systems, Architectures and Processors, Leuven, 2008, pp. 7-12 – njuffa Sep 27 at 6:25
  • I'm vaguely recalling that the IBM 7094 had some weird instructions that might have fit your criteria. – Hot Licks Sep 28 at 0:40
  • 1
    @njuffa: That's a very relevant paper indeed! Link: "Bit Matrix Multiplication in Commodity Processors" (Yedidya Hilewitz, Cédric Lauradoux, Ruby B. Lee; 2008) It sounds like the Cray X1 has a 64x64=4096-bit bmm instruction whose details I don't understand. The paper gives several possible applications of a hypothetical 8x8 BMM instruction, such as "reverse a string of two-bit symbols (e.g. ACGT)," LFSRs, and somehow the FFT. – Quuxplusone Sep 28 at 1:40

I don't know if it counts as a 16-bit historical computer, but there were a number of 80287 clones. One of them, from IIT, added a F4x4 instruction that operated on the entire FP stack to perform matrix operations (the IIT chip had not one, but four stacks of eight FPU registers each).

A friend of mine had one of these coprocessors back in the day, and it came with a demo program showing a speed test of this F4x4 instruction. The leafsheet that served as fitting and user guide also mentioned this feature.

The only information I have found so far is from the COPRO16A.TXT file, by Norbert Juffa. Excerpts from that document include:

The IIT 2C87 provides extra functions not available on any other 287 chip [38]. It has 24 user-accessible floating-point registers organized into three register banks. Additional instructions (FSBP0, FSBP1, FSBP2) allow switching from one bank to another. (Transfers between registers in different banks are not supported, however, so this feature by itself is of limited usefulness. Also, there seems to be only one status register (containing the stack top pointer), so it has to be manually loaded and stored when switching between banks with a different number of registers in use [40]). The register bank's main purpose is to aid the fourth additional instruction the 2C87 has (F4X4), which does a full multiply of a 4x4 matrix by a 4x1 vector, an operation common in 3D- graphics applications [39]. The built-in matrix multiply speeds this operation up by a factor of 6 to 8 when compared to a programmed solution according to the manufacturer [38]. Tests show the speed-up to be indeed in this range [40]. For the 3C87, I measured the execution time of F4X4 to be about 280 clock cycles; the execution time on the 2C87 should be somewhat larger - I estimate it to be around 310 clock cycles due to the higher CPU-NDP communication overhead in instruction execution in 286/287 systems (~45-50 clock cycles) compared with 386/387 systems (~16-20 clock cycles). As desirable as the F4X4 instruction may seem, however, there are very few applications that make use of it when an IIT coprocessor is detected at run time (among them Schroff Development's Silver Screen and Evolution Computing's Fast-CAD 3-D [25]).

And a routine showing how to use the F4x4 instruction:

  ; IIT_MUL_4x4 multiplicates a four-by-four matrix by an array of four
  ; dimensional vectors. This operation is needed for 3D transformations
  ; in graphics data processing. There are arrays for each component of
  ; a vector.  Thus there is an array containing all the x components,
  ; another containing all the y components and so on. Each component is
  ; an 8 byte IEEE floating-point number. Two indices into the array of
  ; vectors are given. The first is the index of the vector that will be
  ; processed first, the second is the index of the vector processed
  ; last. This subroutine uses the special instructions only available
  ; on IIT coprocessors to provide fast matrix multiply capabilities.
  ; So make sure to use it only on IIT coprocessors.

  IIT_MUL_4x4   PROC    NEAR

        AddrX   EQU DWORD PTR [BP+24] ; address of X component array
        AddrY   EQU DWORD PTR [BP+20] ; address of Y component array
        AddrZ   EQU DWORD PTR [BP+16] ; address of Z component array
        AddrW   EQU DWORD PTR [BP+12] ; address of W component array
        AddrT   EQU DWORD PTR [BP+8]  ; addr. of 4x4 transf. matrix
        F       EQU WORD  PTR [BP+6]  ; first vector to process
        K       EQU WORD  PTR [BP+4]  ; last vector to process
        RetAddr EQU WORD  PTR [BP+2]  ; return address saved by call
        SavdBP  EQU WORD  PTR [BP+0]  ; saved frame pointer
        SavdDS  EQU WORD  PTR [BP-2]  ; caller's data segment
        Ctrl87  EQU WORD  PTR [BP-4]  ; caller's 80x87 control word

        PUSH    BP                    ; save TURBO-Pascal frame ptr
        MOV     BP, SP                ; new frame pointer
        PUSH    DS                    ; save TURBO-Pascal data seg.
        SUB     SP, 2                 ; make local variabe
        FSTCW   [Ctrl87]              ; save 80x87 ctrl word
        LES     SI, AddrT             ; ptr to transformation matrix
        FINIT                         ; initialize coprocessor
        FSBP2                         ; set register bank 2
        FLD     QWORD PTR ES:[SI]     ; load a[0,0]
        FLD     QWORD PTR ES:[SI+32]  ; load a[1,0]
        FLD     QWORD PTR ES:[SI+64]  ; load a[2,0]
        FLD     QWORD PTR ES:[SI+96]  ; load a[3,0]
        FLD     QWORD PTR ES:[SI+8]   ; load a[0,1]
        FLD     QWORD PTR ES:[SI+40]  ; load a[1,1]
        FLD     QWORD PTR ES:[SI+72]  ; load a[2,1]
        FLD     QWORD PTR ES:[SI+104] ; load a[3,1]
        FINIT                         ; initialize coprocessor
        FSBP1                         ; set register bank 1
        FLD     QWORD PTR ES:[SI+16]  ; load a[0,2]
        FLD     QWORD PTR ES:[SI+48]  ; load a[1,2]
        FLD     QWORD PTR ES:[SI+80]  ; load a[2,2]
        FLD     QWORD PTR ES:[SI+112] ; load a[3,2]
        FLD     QWORD PTR ES:[SI+24]  ; load a[0,3]
        FLD     QWORD PTR ES:[SI+56]  ; load a[1,3]
        FLD     QWORD PTR ES:[SI+88]  ; load a[2,3]
        FLD     QWORD PTR ES:[SI+120] ; load a[3,3]

                                      ; transformation matrix loaded

        MOV     AX, F                 ; index of first vector
        MOV     DX, K                 ; index of last vector

        MOV     BX, AX                ; index 1st vector to process
        MOV     CL, 3                 ; component has 8 (2**3) bytes
        SHL     BX, CL                ; compute offset into arrays

        FINIT                         ; initialize coprocessor
        FSBP0                         ; set register bank 0

  $mat_loop:LES     SI, AddrW             ; addr. of W component array
        FLD     QWORD PTR ES:[SI+BX]  ; W component current vector
        LES     SI, AddrZ             ; addr. of Z component array
        FLD     QWORD PTR ES:[SI+BX]  ; Z component current vector
        LES     SI, AddrY             ; addr. of Y component array
        FLD     QWORD PTR ES:[SI+BX]  ; Y component current vector
        LES     SI, AddrX             ; addr. of X component array
        FLD     QWORD PTR ES:[SI+BX]  ; X component current vector
        F4X4                          ; mul 4x4 matrix by 4x1 vector
        INC     AX                    ; next vector
        MOV     DI, AX                ; next vector
        SHL     DI, CL                ; offset of vector into arrays

        FSTP    QWORD PTR ES:[SI+BX]  ; store X comp. of curr. vect.
        LES     SI, AddrY             ; address of Y component array
        FSTP    QWORD PTR ES:[SI+BX]  ; store Y comp. of curr. vect.
        LES     SI, AddrZ             ; address of Z component array
        FSTP    QWORD PTR ES:[SI+BX]  ; store Z comp. of curr. vect.
        LES     SI, AddrW             ; address of W component array
        FSTP    QWORD PTR ES:[SI+BX]  ; store W comp. of curr. vect.

        MOV     BX, DI                ; ofs nxt vect. in comp. arrays
        CMP     AX, DX                ; nxt vector past upper bound?
        JLE     $mat_loop             ; no, transform next vector
        FLDCW   [Ctrl87]              ; restore orig 80x87 ctrl word

        ADD      SP, 2                ; get rid of local variable
        POP      DS                   ; restore TP data segment
        POP      BP                   ; restore TP frame pointer
        RET      24                   ; pop parameters and return
  IIT_MUL_4x4   ENDP

  CODE      ENDS

| improve this answer | |
  • 1
    I think I get the general idea, but FWIW this answer could still be improved by including some links (and/or quotations therefrom) explaining the details. – Quuxplusone Sep 26 at 21:38
  • 3
    [Full disclosure: I worked for IIT from 1993 to 1995; on FPU of a "486" design, not on the math coprocessors] Four instructions fsbp0, fsbp1, fsbp2, fsbp3 allowed selecting the currently active bank. Typically, the 4x4 matrix would be loaded into banks 1 and 2. 4-dimensional vectors would then be loaded into bank 0, and f4x4 would perform a matrix-vector operation, with the transformed vector going back into bank0, and bank 3 being used for temp registers. The whole mechanism was described only in an IIT application note / white paper; I don't think I have a copy of that anymore. – njuffa Sep 27 at 5:14
  • @mcleod_ideafix Good find! I had totally forgotten that I had put this 4x4 example code in the document (which I created as a student, ca. 1990-1992, prior to joining IIT as an employee). – njuffa Sep 28 at 7:46

I think the big question you have to ask here is: What would such an instruction have been useful for? Single-bit data types don't usually go together with matrix operations.

I've looked at some examples where the MMIX instructions are used, and they don't seem to be very compelling, given the hardware complexity of doing so. For example, if you want a ROL instruction, you can add a few gates to a SHL instruction's implementation and avoid all the gymnastics needed to emulate it with a MOR instruction. Another example given involves spacing out a narrow character string (eg. ASCII) into a wider one (eg. UTF-16), which strikes me as not a common enough operation to justify using anything more clever than the obvious scalar loop.

The best examples of this type of operation I know of come from relatively recent CPUs with SIMD units, such as ARM NEON (a 32-bit architecture with 64-bit and 128-bit SIMD registers). These often include both "horizontal" (eg. sum all elements in a single register) and "vertical" (eg. multiply two registers element-wise and accumulate into a third) operations that can be used to efficiently support matrix arithmetic by combining a relatively small number of instructions. There are also instructions for interleaving and de-interleaving data (aka. zip and unzip operations) which can be very useful for transposing matrices.

| improve this answer | |
  • I don't know how to research this, but it's possible that the IBM 7094 had an instruction, CRQ (Convert by replacement from the MQ) that might have been used to preform matrix style operations. – Walter Mitty Sep 27 at 0:51
  • 3
    I wouldn’t be so dismissive. The RISC-V bit-manipulation extension, version 0.92 includes bmator, bmatxor and bmatflip; the spec (§2.8) cites not only MMIX, but also Cray XMT and even x86 as precedent. If a number of architectures extend their instruction sets to include this operation, there clearly seems to be some practical demand for it; it’s not like x86 is particularly known for pandering to theorists. – user3840170 Sep 27 at 9:22
  • The Cray BMM capabilities (bit matrix multiply; standard on the Cray SV1, optional on some other models) were used for bioinformatics (mostly sequence alignment, I think) and one could reasonably speculate that there were also applications to cryptography, although nothing is published on that to my knowledge. – njuffa Sep 27 at 9:58
  • @user3840170 - I strongly suspect that the primary application for modern implementations of this are for cryptographic purposes, particularly calculating secure hashes. The requirement to do this quickly is a relatively modern application. – occipita Sep 27 at 10:03
  • 2
    Single-bit data types go great with matrix operations such as bitwise OR, XOR, and AND, and I can easily imagine kinds of prioritization or optimization algorithms benefiting from operations that would set the three lower bits of the destination register to the larger bottom-three-bits value of the two source operands, and do likewise for the other eleven groups of three bits. Such needs probably don't arise so often that they can't be met by using combinations of bitwise operators and shifts (e.g. identify all the bits that are set in A but not B, shift right, and with octal 333333333333... – supercat Sep 27 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.