The sun rising and Algol68's parallel matrix multi[plication...

The sun rose again this morning and presented me with the urge to multiply a matrix with its own inverse... ;-)

I posted the below set of Algol68 Linear Algebra Operators sometime back, code could do with an update to handle sparse matrices... but as it stands it is a cute example of how you can use those extra CPUs and FPUs on your modern laptop...  The next step (ignoring stability) might be to augment with Strassen's O(n^log2(7)) recursive matrix multiplication algorithm, but again using Algol68's parallel processing... O(n^log2(7))/n(cpus) ... Enjoy!

int default upb := 3;
mode field = long real;
mode vector = [default upb]field;
mode matrix = [default upb, default upb]field;

¢ crude exception handling ¢
proc void raise index error := void: goto exception index error;

sema idle cpus = level ( 8 - 1 ); ¢ 8 = number of CPU cores minus parent CPU ¢

¢ define a operators to slice array into quarters ¢
op top = (matrix m)int: ( ⌊m + ⌈m ) %2,
   bot = (matrix m)int: top m + 1,
   left = (matrix m)int: ( 2 ⌊m + 2 ⌈m ) %2,
   right = (matrix m)int: left m + 1,
   left = (vector v)int: ( ⌊v + ⌈v ) %2,
   right = (vector v)int: left v + 1; 
prio top = 8, bot = 8, left = 8, right = 8; ¢ Operator priority - same as LWB & UPB ¢

op × = (vector a, b)field: ( ¢ dot product ¢
  if (⌊a, ⌈a) ≠ (⌊b, ⌈b) then raise index error fi;
  if ⌊a = ⌈a then
    a[⌈a] × b[⌈b]
  else
    field begin, end;
    []proc void schedule=( ¢ of threads ¢
      void: begin:=a[:left a] × b[:left b], 
      void: end  :=a[right a:] × b[right b:]
    );
    if level idle cpus = 0 then ¢ use current CPU ¢
      for thread to ⌈schedule do schedule[thread] od
    else 
      par ( ¢ run vector in parallel ¢
        schedule[1], ¢ assume parent CPU ¢
        ( ↓idle cpus; schedule[2]; ↑idle cpus)
      ) 
    fi;
    begin+end
  fi
);

op × = (matrix a, b)matrix: ¢ matrix multiply ¢
  if (⌊a, 2 ⌊b) = (⌈a, 2 ⌈b) then
    a[⌊a, ] × b[, 2 ⌈b] ¢ dot product ¢
  else
    [⌈a, 2 ⌈b] field out;
    if (2 ⌊a, 2 ⌈a) ≠ (⌊b, ⌈b) then raise index error fi;
    []struct(bool required, proc void thread) schedule = ( ¢ of threads ¢
      ( true, ¢ calculate top left corner ¢
        void: out[:top a, :left b] := a[:top a, ] × b[, :left b]), 
      ( ⌊a ≠ ⌈a, ¢ calculate bottom left corner ¢
        void: out[bot a:, :left b] := a[bot a:, ] × b[, :left b]), 
      ( 2 ⌊b ≠ 2 ⌈b, ¢ calculate top right corner ¢
        void: out[:top a, right b:] := a[:top a, ] × b[, right b:]), 
      ( (⌊a, 2 ⌊b) ≠ (⌈a, 2 ⌈b) , ¢ calculate bottom right corner ¢
        void: out[bot a:, right b:] := a[bot a:, ] × b[, right b:])
    );
    if level idle cpus = 0 then ¢ use current CPU ¢
      for thread to ⌈schedule do (required →schedule[thread] | thread →schedule[thread] ) od
    else 
      par ( ¢ run vector in parallel ¢
        thread →schedule[1], ¢ thread is always required, and assume parent CPU ¢
        ( required →schedule[4] | ↓idle cpus; thread →schedule[4]; ↑idle cpus),
           ¢ try to do opposite corners of matrix in parallel if CPUs are limited ¢
        ( required →schedule[3] | ↓idle cpus; thread →schedule[3]; ↑idle cpus),
        ( required →schedule[2] | ↓idle cpus; thread →schedule[2]; ↑idle cpus)
      )
    fi;
    out
  fi;

format real fmt = $g(-6,2)$; ¢ width of 6, with no '+' sign, 2 decimals ¢
proc real matrix printf= (format real fmt, matrix m)void:(
  format vector fmt = $"("n(2 ⌈m-1)(f(real fmt)",")f(real fmt)")"$;
  format matrix fmt = $x"("n(⌈m-1)(f(vector fmt)","lxx)f(vector fmt)");"$;
  ¢ finally print the result ¢
  printf((matrix fmt,m))
);

¢ Some sample matrices to test ¢
matrix a=((1,  1,  1,   1), ¢ matrix A ¢
          (2,  4,  8,  16),
          (3,  9, 27,  81),
          (4, 16, 64, 256));

matrix b=((  4  , -3  ,  4/3,  -1/4 ), ¢ matrix B ¢
          (-13/3, 19/4, -7/3,  11/24),
          (  3/2, -2  ,  7/6,  -1/4 ),
          ( -1/6,  1/4, -1/6,   1/24));

matrix c = a × b; ¢ actual multiplication example of A x B ¢

print((" A x B = Product of a and b",new line));
real matrix printf(real fmt, c).

exception index error: 
  putf(stand error, $x"Exception: index error."l$)

Output:

 A x B = Product of a and b
  (( 1.00, -0.00, -0.00, -0.00),
   ( -0.00, 1.00, -0.00, -0.00),
   ( -0.00, -0.00, 1.00, -0.00),
   ( -0.00, -0.00, -0.00, 1.00));