Matrix multiplication: Difference between revisions

Line 142:

Alternatively, for multicore CPUs, the following recursive algorithm can be used:

~~~~int~~~~ default upb := 3;

'''int''' default upb := 3;

~~~~mode~~~~ ~~~~field~~~~ = ~~~~long~~~~ ~~~~real~~~~;

'''mode''' '''field''' = '''long''' '''real''';

~~~~mode~~~~ ~~~~vector~~~~ = [default upb]~~~~field~~~~;

'''mode''' '''vector''' = [default upb]'''field''';

~~~~mode~~~~ ~~~~matrix~~~~ = [default upb, default upb]~~~~field~~~~;

'''mode''' '''matrix''' = [default upb, default upb]'''field''';

¢ crude exception handling ¢

~~~~proc~~~~ ~~~~void~~~~ raise index error := ~~~~void~~~~: ~~~~goto~~~~ exception index error;

'''proc''' '''void''' raise index error := '''void''': '''goto''' exception index error;

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

'''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~~~~: ( &lfloor;m + &lceil;m ) %2,

'''op''' '''top''' = ('''matrix''' m)'''int''': ( &lfloor;m + &lceil;m ) %2,

~~~~bot~~~~ = (~~~~matrix~~~~ m)~~~~int~~~~: ~~~~top~~~~ m + 1,

'''bot''' = ('''matrix''' m)'''int''': '''top''' m + 1,

~~~~left~~~~ = (~~~~matrix~~~~ m)~~~~int~~~~: ( 2 &lfloor;m + 2 &lceil;m ) %2,

'''left''' = ('''matrix''' m)'''int''': ( 2 &lfloor;m + 2 &lceil;m ) %2,

~~~~right~~~~ = (~~~~matrix~~~~ m)~~~~int~~~~: ~~~~left~~~~ m + 1,

'''right''' = ('''matrix''' m)'''int''': '''left''' m + 1,

~~~~left~~~~ = (~~~~vector~~~~ v)~~~~int~~~~: ( &lfloor;v + &lceil;v ) %2,

'''left''' = ('''vector''' v)'''int''': ( &lfloor;v + &lceil;v ) %2,

~~~~right~~~~ = (~~~~vector~~~~ v)~~~~int~~~~: ~~~~left~~~~ v + 1;

'''right''' = ('''vector''' v)'''int''': '''left''' v + 1;

~~~~prio~~~~ ~~~~top~~~~ = 8, ~~~~bot~~~~ = 8, ~~~~left~~~~ = 8, ~~~~right~~~~ = 8; ¢ Operator priority - same as LWB & UPB ¢

'''prio''' '''top''' = 8, '''bot''' = 8, '''left''' = 8, '''right''' = 8; ¢ Operator priority - same as LWB & UPB ¢

~~~~op~~~~ × = (~~~~vector~~~~ a, b)~~~~field~~~~: ( ¢ dot product ¢

'''op''' × = ('''vector''' a, b)'''field''': ( ¢ dot product ¢

~~~~if~~~~ (&lfloor;a, &lceil;a) ≠ (&lfloor;b, &lceil;b) ~~~~then~~~~ raise index error ~~~~fi~~~~;

'''if''' (&lfloor;a, &lceil;a) ≠ (&lfloor;b, &lceil;b) '''then''' raise index error '''fi''';

~~~~if~~~~ &lfloor;a = &lceil;a ~~~~then~~~~

'''if''' &lfloor;a = &lceil;a '''then'''

a[&lceil;a] × b[&lceil;b]

~~~~else~~~~

'''else'''

~~~~field~~~~ begin, end;

'''field''' begin, end;

[]~~~~proc~~~~ ~~~~void~~~~ schedule=(

[]'''proc''' '''void''' schedule=(

~~~~void~~~~: begin:=a[:~~~~left~~~~ a] × b[:~~~~left~~~~ b],

'''void''': begin:=a[:'''left''' a] × b[:'''left''' b],

~~~~void~~~~: end :=a[~~~~right~~~~ a:] × b[~~~~right~~~~ b:]

'''void''': end :=a['''right''' a:] × b['''right''' b:]

);

~~~~if~~~~ ~~~~level~~~~ idle cpus = 0 ~~~~then~~~~ ¢ use current CPU ¢

'''if''' '''level''' idle cpus = 0 '''then''' ¢ use current CPU ¢

~~~~for~~~~ thread ~~~~to~~~~ &lceil;schedule ~~~~do~~~~ schedule[thread] ~~~~od~~~~

'''for''' thread '''to''' &lceil;schedule '''do''' schedule[thread] '''od'''

~~~~else~~~~

'''else'''

~~~~par~~~~ ( ¢ run vector in parallel ¢

'''par''' ( ¢ run vector in parallel ¢

schedule[1], ¢ assume parent CPU ¢

( ↓idle cpus; schedule[2]; ↑idle cpus)

)

~~~~fi~~~~;

'''fi''';

begin+end

~~~~fi~~~~

'''fi'''

);

~~~~op~~~~ × = (~~~~matrix~~~~ a, b)~~~~matrix~~~~: ¢ matrix multiply ¢

'''op''' × = ('''matrix''' a, b)'''matrix''': ¢ matrix multiply ¢

~~~~if~~~~ (&lfloor;a, 2 &lfloor;b) = (&lceil;a, 2 &lceil;b) ~~~~then~~~~

'''if''' (&lfloor;a, 2 &lfloor;b) = (&lceil;a, 2 &lceil;b) '''then'''

a[&lfloor;a, ] × b[, 2 &lceil;b] ¢ dot product ¢

~~~~else~~~~

'''else'''

[&lceil;a, 2 &lceil;b] ~~~~field~~~~ out;

[&lceil;a, 2 &lceil;b] '''field''' out;

~~~~if~~~~ (2 &lfloor;a, 2 &lceil;a) ≠ (&lfloor;b, &lceil;b) ~~~~then~~~~ raise index error ~~~~fi~~~~;

'''if''' (2 &lfloor;a, 2 &lceil;a) ≠ (&lfloor;b, &lceil;b) '''then''' raise index error '''fi''';

[]~~~~struct~~~~(~~~~bool~~~~ required, ~~~~proc~~~~ ~~~~void~~~~ thread) schedule = (

[]'''struct'''('''bool''' required, '''proc''' '''void''' thread) schedule = (

( ~~~~true~~~~, ¢ calculate top left corner ¢

( '''true''', ¢ calculate top left corner ¢

~~~~void~~~~: out[:~~~~top~~~~ a, :~~~~left~~~~ b] := a[:~~~~top~~~~ a, ] × b[, :~~~~left~~~~ b]),

'''void''': out[:'''top''' a, :'''left''' b] := a[:'''top''' a, ] × b[, :'''left''' b]),

( &lfloor;a ≠ &lceil;a, ¢ calculate bottom left corner ¢

~~~~void~~~~: out[~~~~bot~~~~ a:, :~~~~left~~~~ b] := a[~~~~bot~~~~ a:, ] × b[, :~~~~left~~~~ b]),

'''void''': out['''bot''' a:, :'''left''' b] := a['''bot''' a:, ] × b[, :'''left''' b]),

( 2 &lfloor;b ≠ 2 &lceil;b, ¢ calculate top right corner ¢

~~~~void~~~~: out[:~~~~top~~~~ a, ~~~~right~~~~ b:] := a[:~~~~top~~~~ a, ] × b[, ~~~~right~~~~ b:]),

'''void''': out[:'''top''' a, '''right''' b:] := a[:'''top''' a, ] × b[, '''right''' b:]),

( (&lfloor;a, 2 &lfloor;b) ≠ (&lceil;a, 2 &lceil;b) , ¢ calculate bottom right corner ¢

~~~~void~~~~: out[~~~~bot~~~~ a:, ~~~~right~~~~ b:] := a[~~~~bot~~~~ a:, ] × b[, ~~~~right~~~~ b:])

'''void''': out['''bot''' a:, '''right''' b:] := a['''bot''' a:, ] × b[, '''right''' b:])

);

~~~~if~~~~ ~~~~level~~~~ idle cpus = 0 ~~~~then~~~~ ¢ use current CPU ¢

'''if''' '''level''' idle cpus = 0 '''then''' ¢ use current CPU ¢

~~~~for~~~~ thread ~~~~to~~~~ &lceil;schedule ~~~~do~~~~ (required →schedule[thread] | thread →schedule[thread] ) ~~~~od~~~~

'''for''' thread '''to''' &lceil;schedule '''do''' (required →schedule[thread] | thread →schedule[thread] ) '''od'''

~~~~else~~~~

'''else'''

~~~~par~~~~ ( ¢ run vector in parallel ¢

'''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),

Line 209:

( required →schedule[2] | ↓idle cpus; thread →schedule[2]; ↑idle cpus)

)

~~~~fi~~~~;

'''fi''';

out

~~~~fi~~~~;

'''fi''';

~~~~format~~~~ real fmt = $g(-6,2)$; ¢ width of 6, with no '+' sign, 2 decimals ¢

'''format''' real fmt = $g(-6,2)$; ¢ width of 6, with no '+' sign, 2 decimals ¢

~~~~proc~~~~ real matrix printf= (~~~~format~~~~ real fmt, ~~~~matrix~~~~ m)~~~~void~~~~:(

'''proc''' real matrix printf= ('''format''' real fmt, '''matrix''' m)'''void''':(

~~~~format~~~~ vector fmt = $"("n(2 &lceil;m-1)(f(real fmt)",")f(real fmt)")"$;

'''format''' vector fmt = $"("n(2 &lceil;m-1)(f(real fmt)",")f(real fmt)")"$;

~~~~format~~~~ matrix fmt = $x"("n(&lceil;m-1)(f(vector fmt)","lxx)f(vector fmt)");"$;

'''format''' matrix fmt = $x"("n(&lceil;m-1)(f(vector fmt)","lxx)f(vector fmt)");"$;

¢ finally print the result ¢

printf((matrix fmt,m))

Line 222:

¢ Some sample matrices to test ¢

~~~~matrix~~~~ a=((1, 1, 1, 1), ¢ matrix A ¢

'''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 ¢

'''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 ¢

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

print((" A x B =",new line));

real matrix printf(real fmt, c) ∎

real matrix printf(real fmt, c).

exception index error: