\ mmul.4th
\
\ Double precision floating point matrix multiplication
\ for an integrated data/fp stack.
\
\ Krishna Myneni, Creative Consulting for Research
\ and Education
\
\ Usage: a1 a2 a3 nr1 nc1 nc2 df_mmul
\
\ Requires:
\   ans-words.4th
\
\ Revisions:
\   2017-05-17  km; created.
\   2017-05-18  km; precompute row offsets for a1 and a2.
\
\ Notes:
\   0. Matrix data is assumed to be stored in row order
\
\   1. Only the word DF_MMUL is specific for an integrated
\      data/fp stack. Other words work with separate fp
\      stack.

variable nc1
variable nc2
variable a1
variable a2
variable roffs1
variable roffs2

\ Convert row of a1 and column of a2 to
\ corresponding addresses
: df_r1c2>a1a2 ( row1 col2 -- arow1 acol2 )
    floats a2 @ + >r
    roffs1 @ * a1 @ + r>   \ -- arow1 acol2
;

0 [if]
\ Multiply row of a1 with col of a2, element by element,
\ and accumulate the sum.
: df_mul_r1c2 ( row1 col2 -- rsum )
    df_r1c2>a1a2
    2>r 0e 2r>  \ rsum a1 a2
    nc1 @ 0 DO
      2dup 2>r >r f@ r> f@ f* f+
      2r> roffs2 @ + >r
      float+ r>
    LOOP
    2drop
;
[then]

1 [if]
\ Multiply row of a1 with col of a2, element by element,
\ and accumulate the sum.
: df_mul_r1c2 ( row1 col2 -- rsum )
    df_r1c2>a1a2
    2>r 0e 2r>  \ rsum a1 a2
    nc1 @ 0 DO
      2dup f@ f@ f* f+
      roffs2 @ + swap
      float+ swap
    LOOP
    2drop
;
[then]

0 [if]
: df_mul_r1c2 ( row1 col2 -- rsum )
    df_r1c2>a1a2 [ 1 floats ] literal swap roffs2 @ nc1 @ v* ;
[then]

\ Multiply two double-precision matrices with data beginning at
\ a1 and a2, and store at a3. Proper memory allocation is
\ assumed, as are the dimensions for a2, i.e. nr2 = nc1 is
\ assumed. This word assumes an integrated data/fp stack.
: df_mmul ( a1 a2 a3 nr1 nc1 nc2 -- )
    nc2 ! nc1 ! 2>r a2 ! a1 !

    \ offsets to next row for a1 and a2
    nc1 @ floats roffs1 !
    nc2 @ floats roffs2 !

    2r>  \ a3 nr1
    0 DO
      nc2 @ 0 DO
        J I rot >r df_mul_r1c2 r@ f!
        r> float+
      LOOP
    LOOP
    drop
;

: matmulr df_mmul ;

\ Test code
\
\ Requires the additional files:
\   ttester.4th
\   fsl/fsl-util.4th

false [IF]
base @
decimal

\ Allot and initialize three 2x2 matrices

2 2 double matrix a{{
2 2 double matrix b{{
2 2 double matrix c{{

t{ 1e a{{ 0 0 }} f! ->  }t
t{ 2e a{{ 0 1 }} f! ->  }t
t{ 3e a{{ 1 0 }} f! ->  }t
t{ 4e a{{ 1 1 }} f! ->  }t

t{ 5e b{{ 0 0 }} f! ->  }t
t{ 6e b{{ 0 1 }} f! ->  }t
t{ 7e b{{ 1 0 }} f! ->  }t
t{ 8e b{{ 1 1 }} f! ->  }t

TESTING df_mmul
t{ a{{ 0 0 }} b{{ 0 0 }} c{{ 0 0 }} 2 2 2 df_mmul -> }t
t{ c{{ 0 0 }} f@  ->  19e }t
t{ c{{ 0 1 }} f@  ->  22e }t
t{ c{{ 1 0 }} f@  ->  43e }t
t{ c{{ 1 1 }} f@  ->  50e }t

base !
[THEN]