# User:NevilleDNZ/Matrix multiplication.a68

```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: exception index error;

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

¢ define an operator 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=(
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 ¢
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 = (
( 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 ¢
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 =",new line));
real matrix printf(real fmt, c) ∎

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