Matrix multiplication: Difference between revisions

m
→‎[[Matrix multiplication#ALGOL 68]]: use BOLD as per specification.
m (→‎[[Matrix multiplication#ELLA]]: move to E section.)
m (→‎[[Matrix multiplication#ALGOL 68]]: use BOLD as per specification.)
Line 142:
 
Alternatively, for multicore CPUs, the following recursive algorithm can be used:
<u>'''int</u>''' default upb := 3;
<u>'''mode</u>''' <u>'''field</u>''' = <u>'''long</u>''' <u>'''real</u>''';
<u>'''mode</u>''' <u>'''vector</u>''' = [default upb]<u>'''field</u>''';
<u>'''mode</u>''' <u>'''matrix</u>''' = [default upb, default upb]<u>'''field</u>''';
&cent; crude exception handling &cent;
<u>'''proc</u>''' <u>'''void</u>''' raise index error := <u>'''void</u>''': <u>'''goto</u>''' exception index error;
<u>'''sema</u>''' idle cpus = <u>'''level</u>''' ( 8 - 1 ); &cent; 8 = number of CPU cores minus parent CPU &cent;
&cent; define an operator to slice array into quarters &cent;
<u>'''op</u>''' <u>'''top</u>''' = (<u>'''matrix</u>''' m)<u>'''int</u>''': ( &lfloor;m + &lceil;m ) %2,
<u>'''bot</u>''' = (<u>'''matrix</u>''' m)<u>'''int</u>''': <u>'''top</u>''' m + 1,
<u>'''left</u>''' = (<u>'''matrix</u>''' m)<u>'''int</u>''': ( 2 &lfloor;m + 2 &lceil;m ) %2,
<u>'''right</u>''' = (<u>'''matrix</u>''' m)<u>'''int</u>''': <u>'''left</u>''' m + 1,
<u>'''left</u>''' = (<u>'''vector</u>''' v)<u>'''int</u>''': ( &lfloor;v + &lceil;v ) %2,
<u>'''right</u>''' = (<u>'''vector</u>''' v)<u>'''int</u>''': <u>'''left</u>''' v + 1;
<u>'''prio</u>''' <u>'''top</u>''' = 8, <u>'''bot</u>''' = 8, <u>'''left</u>''' = 8, <u>'''right</u>''' = 8; &cent; Operator priority - same as LWB & UPB &cent;
<u>'''op</u>''' &times; = (<u>'''vector</u>''' a, b)<u>'''field</u>''': ( &cent; dot product &cent;
<u>'''if</u>''' (&lfloor;a, &lceil;a) &ne; (&lfloor;b, &lceil;b) <u>'''then</u>''' raise index error <u>'''fi</u>''';
<u>'''if</u>''' &lfloor;a = &lceil;a <u>'''then</u>'''
a[&lceil;a] &times; b[&lceil;b]
<u>'''else</u>'''
<u>'''field</u>''' begin, end;
[]<u>'''proc</u>''' <u>'''void</u>''' schedule=(
<u>'''void</u>''': begin:=a[:<u>'''left</u>''' a] &times; b[:<u>'''left</u>''' b],
<u>'''void</u>''': end :=a[<u>'''right</u>''' a:] &times; b[<u>'''right</u>''' b:]
);
<u>'''if</u>''' <u>'''level</u>''' idle cpus = 0 <u>'''then</u>''' &cent; use current CPU &cent;
<u>'''for</u>''' thread <u>'''to</u>''' &lceil;schedule <u>'''do</u>''' schedule[thread] <u>'''od</u>'''
<u>'''else</u>'''
<u>'''par</u>''' ( &cent; run vector in parallel &cent;
schedule[1], &cent; assume parent CPU &cent;
( &darr;idle cpus; schedule[2]; &uarr;idle cpus)
)
<u>'''fi</u>''';
begin+end
<u>'''fi</u>'''
);
<u>'''op</u>''' &times; = (<u>'''matrix</u>''' a, b)<u>'''matrix</u>''': &cent; matrix multiply &cent;
<u>'''if</u>''' (&lfloor;a, 2 &lfloor;b) = (&lceil;a, 2 &lceil;b) <u>'''then</u>'''
a[&lfloor;a, ] &times; b[, 2 &lceil;b] &cent; dot product &cent;
<u>'''else</u>'''
[&lceil;a, 2 &lceil;b] <u>'''field</u>''' out;
<u>'''if</u>''' (2 &lfloor;a, 2 &lceil;a) &ne; (&lfloor;b, &lceil;b) <u>'''then</u>''' raise index error <u>'''fi</u>''';
[]<u>'''struct</u>'''(<u>'''bool</u>''' required, <u>'''proc</u>''' <u>'''void</u>''' thread) schedule = (
( <u>'''true</u>''', &cent; calculate top left corner &cent;
<u>'''void</u>''': out[:<u>'''top</u>''' a, :<u>'''left</u>''' b] := a[:<u>'''top</u>''' a, ] &times; b[, :<u>'''left</u>''' b]),
( &lfloor;a &ne; &lceil;a, &cent; calculate bottom left corner &cent;
<u>'''void</u>''': out[<u>'''bot</u>''' a:, :<u>'''left</u>''' b] := a[<u>'''bot</u>''' a:, ] &times; b[, :<u>'''left</u>''' b]),
( 2 &lfloor;b &ne; 2 &lceil;b, &cent; calculate top right corner &cent;
<u>'''void</u>''': out[:<u>'''top</u>''' a, <u>'''right</u>''' b:] := a[:<u>'''top</u>''' a, ] &times; b[, <u>'''right</u>''' b:]),
( (&lfloor;a, 2 &lfloor;b) &ne; (&lceil;a, 2 &lceil;b) , &cent; calculate bottom right corner &cent;
<u>'''void</u>''': out[<u>'''bot</u>''' a:, <u>'''right</u>''' b:] := a[<u>'''bot</u>''' a:, ] &times; b[, <u>'''right</u>''' b:])
);
<u>'''if</u>''' <u>'''level</u>''' idle cpus = 0 <u>'''then</u>''' &cent; use current CPU &cent;
<u>'''for</u>''' thread <u>'''to</u>''' &lceil;schedule <u>'''do</u>''' (required &rarr;schedule[thread] | thread &rarr;schedule[thread] ) <u>'''od</u>'''
<u>'''else</u>'''
<u>'''par</u>''' ( &cent; run vector in parallel &cent;
thread &rarr;schedule[1], &cent; thread is always required, and assume parent CPU &cent;
( required &rarr;schedule[4] | &darr;idle cpus; thread &rarr;schedule[4]; &uarr;idle cpus),
Line 209:
( required &rarr;schedule[2] | &darr;idle cpus; thread &rarr;schedule[2]; &uarr;idle cpus)
)
<u>'''fi</u>''';
out
<u>'''fi</u>''';
<u>'''format</u>''' real fmt = $g(-6,2)$; &cent; width of 6, with no '+' sign, 2 decimals &cent;
<u>'''proc</u>''' real matrix printf= (<u>'''format</u>''' real fmt, <u>'''matrix</u>''' m)<u>'''void</u>''':(
<u>'''format</u>''' vector fmt = $"("n(2 &lceil;m-1)(f(real fmt)",")f(real fmt)")"$;
<u>'''format</u>''' matrix fmt = $x"("n(&lceil;m-1)(f(vector fmt)","lxx)f(vector fmt)");"$;
&cent; finally print the result &cent;
printf((matrix fmt,m))
Line 222:
&cent; Some sample matrices to test &cent;
<u>'''matrix</u>''' a=((1, 1, 1, 1), &cent; matrix A &cent;
(2, 4, 8, 16),
(3, 9, 27, 81),
(4, 16, 64, 256));
<u>'''matrix</u>''' b=(( 4 , -3 , 4/3, -1/4 ), &cent; matrix B &cent;
(-13/3, 19/4, -7/3, 11/24),
( 3/2, -2 , 7/6, -1/4 ),
( -1/6, 1/4, -1/6, 1/24));
<u>'''matrix</u>''' c = a &times; b; &cent; actual multiplication example of A x B &cent;
print((" A x B =",new line));
real matrix printf(real fmt, c).
exception index error: