( MatrixMult)
( 17th March 2016)

( A: 2 x 3; B: 3 x 4; C: 2 x 4   -  C is calculated result of multiplying matrices A and B)
( Need to have #rows & #columns embedded in 2D array - to use as variables later)
( Integer arithmetic used - see eg calc_det4 for using floating point)


0 VARIABLE n   ( handy variable)

( For <BUILDS: #rows,#columns.....) ( word array)
: 2Darray  ( row,column....#columns,#rows,actual addr of {row,column} element)
  <BUILDS 2DUP SWAP , , * 4 * ALLOT DOES>   ( row,col,addr...)
  DUP @ SWAP 4+ DUP @ SWAP 4+      ( row,col,#rows,#columns,start addr of array...)
  >R SWAP DROP                     ( row,col,#columns...)
  ROT * +                          ( byte offset...)
  4*                               ( word offset...)
  R> + ;                           ( actual addr of row,column element...)


: array_dims  4- DUP @ SWAP 4- @ ;  ( eg 0 0 fred ... #columns,#rows)

( define sizes of 2D arrays A,B and C)
2 3 2Darray A
3 4 2Darray B
2 4 2Darray C

( load A using stored #rows,#columns)
: load_A  ( a0->an,row)
  0 0 A array_dims  DROP ( #columns...)
  0 SWAP 1- DO DUP >R   I A ! R>
   -1 +LOOP DROP ;

5 6 1   0 load_A    ( define and load row 0 for A)
7 1 3   1 load_A    ( define and load row 1 for A)

( load B using stored #rows,#columns)
: load_B  ( a0->an,row)
  0 0 B array_dims  DROP ( #columns...)
  0 SWAP 1- DO DUP >R   I B ! R>
   -1 +LOOP DROP ;

2  3  1  2  0 load_B
1  1  0  0  1 load_B
0  0 -1 -1  2 load_B

( load C using stored #rows,#columns)
: load_C  ( a0->an,row)
  0 0 C array_dims  DROP ( #columns...)
  0 SWAP 1- DO DUP >R   I C ! R>
   -1 +LOOP DROP ;

( clear matrix C prior to calculation)
: zero_C
 0 0 0 0  0 load_C
 0 0 0 0  1 load_C ;

( words to display each matrix)
: dsp_A 0 0 A array_dims  0 DO CR DUP 0 DO J I A @ . LOOP LOOP DROP ;
: dsp_B 0 0 B array_dims  0 DO CR DUP 0 DO J I B @ . LOOP LOOP DROP ;
: dsp_C 0 0 C array_dims  0 DO CR DUP 0 DO J I C @ . LOOP LOOP DROP ;
: dsp_ABC CR ." Matrix A:" dsp_A CR
          CR ." Matrix B:" dsp_B CR
          CR ." Resultant Matrix C:" dsp_C CR ;

: A*Brow             ( #cB,#rB,#rB,#cB...)
     0 DO DUP
         0 DO          
          I    J B @   ( DUP CR ." B: " .  I . J .)  ( used for checking)
          n @  I A @   ( DUP CR ." A: " .  ."  n=" n ?)
          *            ( DUP CR ." A*B: " . )
          n @  J C +!  ( n @ J C ?  KEY DROP CR)
        LOOP
      LOOP DROP ;

: A*B
  CR ." Matrix Multiplication:  C = A*B" CR
  zero_C
  0 0 A array_dims  ( #columns,#rows...) ( for A) 
  0 0 B array_dims  ( #columns,#rows...) ( for B) 
  2SWAP  SWAP DROP  ( #cB,#rB,#rA...)
  0 DO I n !        ( n is the row number in A)
     2DUP           ( #cB,#rB,#cB,#rB...)
     SWAP A*Brow
    LOOP 2DROP
    dsp_ABC ;


