Matrix multiplication: Difference between revisions
Content added Content deleted
m (→[[ALGOL 68]]: use the bold form of the if ~ then ~ else ~ fi syntax. - easier to read blocks.) |
m (→[[ALGOL 68]]) |
||
Line 37: | Line 37: | ||
MODE FIELD = LONG REAL; # field type is LONG REAL # |
MODE FIELD = LONG REAL; # field type is LONG REAL # |
||
INT default upb:=3; |
INT default upb:=3; |
||
MODE |
MODE VECTOR = [default upb]FIELD; |
||
MODE |
MODE MATRIX = [default upb,default upb]FIELD; |
||
# crude exception handling # |
# crude exception handling # |
||
Line 44: | Line 44: | ||
# define the vector/matrix operators # |
# define the vector/matrix operators # |
||
OP * = ( |
OP * = (VECTOR a,b)FIELD: ( # basically the dot product # |
||
FIELD result:=0; |
FIELD result:=0; |
||
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
||
Line 51: | Line 51: | ||
); |
); |
||
OP * = ( |
OP * = (VECTOR a, MATRIX b)VECTOR: ( # overload vector times matrix # |
||
[2 LWB b:2 UPB b]FIELD result; |
[2 LWB b:2 UPB b]FIELD result; |
||
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
||
Line 59: | Line 59: | ||
<pre style="background-color:#ffe"> |
<pre style="background-color:#ffe"> |
||
# this is the task portion # |
# this is the task portion # |
||
OP * = ( |
OP * = (MATRIX a, b)MATRIX: ( # overload matrix times matrix # |
||
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result; |
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result; |
||
IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI; |
IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI; |
||
Line 67: | Line 67: | ||
</pre> |
</pre> |
||
# Some sample matrices to test # |
# Some sample matrices to test # |
||
MATRIX a=((1, 1, 1, 1), # matrix A # |
|||
(2, 4, 8, 16), |
(2, 4, 8, 16), |
||
(3, 9, 27, 81), |
(3, 9, 27, 81), |
||
(4, 16, 64, 256)); |
(4, 16, 64, 256)); |
||
MATRIX b=(( 4 , -3 , 4/3, -1/4 ), # matrix B # |
|||
(-13/3, 19/4, -7/3, 11/24), |
(-13/3, 19/4, -7/3, 11/24), |
||
( 3/2, -2 , 7/6, -1/4 ), |
( 3/2, -2 , 7/6, -1/4 ), |
||
( -1/6, 1/4, -1/6, 1/24)); |
( -1/6, 1/4, -1/6, 1/24)); |
||
MATRIX prod = a * b; # actual multiplication example of A x B # |
|||
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, |
PROC real matrix printf= (FORMAT real fmt, MATRIX m)VOID:( |
||
FORMAT |
FORMAT vector fmt = $"("n(2 UPB m-1)(f(real fmt)",")f(real fmt)")"$; |
||
FORMAT matrix fmt = $x"("n(UPB m-1)(f( |
FORMAT matrix fmt = $x"("n(UPB m-1)(f(vector fmt)","lxx)f(vector fmt)");"$; |
||
# finally print the result # |
# finally print the result # |
||
printf((matrix fmt,m)) |
printf((matrix fmt,m)) |
||
Line 102: | Line 102: | ||
Alternatively, for multicore CPUs, the following recursive algorithm can be used: |
Alternatively, for multicore CPUs, the following recursive algorithm can be used: |
||
<u>int</u> default upb := |
<u>int</u> default upb := 3; |
||
<u>mode</u> <u>field</u> = <u>long</u> <u>real</u>; |
<u>mode</u> <u>field</u> = <u>long</u> <u>real</u>; |
||
<u>mode</u> <u> |
<u>mode</u> <u>vector</u> = [default upb]<u>field</u>; |
||
<u>mode</u> <u> |
<u>mode</u> <u>matrix</u> = [default upb, default upb]<u>field</u>; |
||
# crude exception handling # |
# crude exception handling # |
||
Line 113: | Line 113: | ||
# define an operator to slice array into quarters # |
# define an operator to slice array into quarters # |
||
<u>op</u> <u>top</u> = (<u> |
<u>op</u> <u>top</u> = (<u>matrix</u> m)<u>int</u>: ( ⌊m + ⌈m ) %2, |
||
<u>bot</u> = (<u> |
<u>bot</u> = (<u>matrix</u> m)<u>int</u>: <u>top</u> m + 1, |
||
<u>left</u> = (<u> |
<u>left</u> = (<u>matrix</u> m)<u>int</u>: ( 2⌊m + 2⌈m ) %2, |
||
<u>right</u> = (<u> |
<u>right</u> = (<u>matrix</u> m)<u>int</u>: <u>left</u> m + 1, |
||
<u>left</u> = (<u> |
<u>left</u> = (<u>vector</u> v)<u>int</u>: ( ⌊v + ⌈v ) %2, |
||
<u>right</u> = (<u> |
<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; # Operator priority - same as LWB & UPB # |
<u>prio</u> <u>top</u> = 8, <u>bot</u> = 8, <u>left</u> = 8, <u>right</u> = 8; # Operator priority - same as LWB & UPB # |
||
<u>op</u> × = (<u> |
<u>op</u> × = (<u>vector</u> a, b)<u>field</u>: ( # dot product # |
||
<u>if</u> |
<u>if</u> (⌊a, ⌈a) ≠ (⌊b, ⌈b) <u>then</u> raise index error <u>fi</u>; |
||
<u>if</u> ⌊a = ⌈a <u>then</u> |
<u>if</u> ⌊a = ⌈a <u>then</u> |
||
a[⌈a] × b[⌈b] |
a[⌈a] × b[⌈b] |
||
<u>else</u> |
<u>else</u> |
||
<u>field</u> begin, end; |
<u>field</u> begin, end; |
||
[]<u>proc</u> <u>void</u> |
[]<u>proc</u> <u>void</u> schedule=( |
||
<u>void</u>: begin:=a[:<u>left</u> a] × b[:<u>left</u> b], |
<u>void</u>: begin:=a[:<u>left</u> a] × b[:<u>left</u> b], |
||
<u>void</u>: end :=a[<u>right</u> a:] × b[<u>right</u> b:] |
<u>void</u>: end :=a[<u>right</u> a:] × b[<u>right</u> b:] |
||
); |
); |
||
<u>if</u> <u>level</u> idle cpus = 0 <u>then</u> # use current CPU # |
<u>if</u> <u>level</u> idle cpus = 0 <u>then</u> # use current CPU # |
||
<u>for</u> thread <u>to</u> |
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> schedule[thread] <u>od</u> |
||
<u>else</u> |
<u>else</u> |
||
# PAR # ( # run vector in parallel # |
|||
schedule[1], # assume parent CPU # |
|||
( ↓idle cpus; |
( ↓idle cpus; schedule[2]; ↑idle cpus) |
||
) |
) |
||
<u>fi</u>; |
<u>fi</u>; |
||
Line 143: | Line 143: | ||
); |
); |
||
<u>op</u> × = (<u> |
<u>op</u> × = (<u>matrix</u> a, b)<u>matrix</u>: # matrix multiply # |
||
<u>if</u> ⌊a = ⌈a |
<u>if</u> (⌊a, 2⌊b) = (⌈a, 2⌈b) <u>then</u> |
||
a[⌊a, ] × b[, |
a[⌊a, ] × b[, 2⌈b] # dot product # |
||
<u>else</u> |
<u>else</u> |
||
[⌈a, |
[⌈a, 2⌈b] <u>field</u> out; |
||
<u>if</u> |
<u>if</u> (2⌊a, 2⌈a) ≠ (⌊b, ⌈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>struct</u>(<u>bool</u> required, <u>proc</u> <u>void</u> thread) schedule = ( |
||
( <u>true</u>, # calculate top left corner # |
( <u>true</u>, # calculate top left corner # |
||
Line 154: | Line 154: | ||
( ⌊a ≠ ⌈a, # calculate bottom left corner # |
( ⌊a ≠ ⌈a, # calculate bottom left corner # |
||
<u>void</u>: out[<u>bot</u> a:, :<u>left</u> b] := a[<u>bot</u> a:, ] × b[, :<u>left</u> b]), |
<u>void</u>: out[<u>bot</u> a:, :<u>left</u> b] := a[<u>bot</u> a:, ] × b[, :<u>left</u> b]), |
||
( |
( 2⌊b ≠ 2⌈b, # calculate top right corner # |
||
<u>void</u>: out[:<u>top</u> a, <u>right</u> b:] := a[:<u>top</u> a, ] × b[, <u>right</u> b:]), |
<u>void</u>: out[:<u>top</u> a, <u>right</u> b:] := a[:<u>top</u> a, ] × b[, <u>right</u> b:]), |
||
( |
( (⌊a, 2⌊b) ≠ (⌈a, 2⌈b) , # calculate bottom right corner # |
||
<u>void</u>: out[<u>bot</u> a:, <u>right</u> b:] := a[<u>bot</u> a:, ] × b[, <u>right</u> b:]) |
<u>void</u>: out[<u>bot</u> a:, <u>right</u> b:] := a[<u>bot</u> a:, ] × b[, <u>right</u> b:]) |
||
); |
); |
||
<u>if</u> <u>level</u> idle cpus = 0 <u>then</u> # use current CPU # |
<u>if</u> <u>level</u> idle cpus = 0 <u>then</u> # use current CPU # |
||
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> ( |
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> (required→schedule[thread] | thread→schedule[thread] ) <u>od</u> |
||
<u>else</u> |
<u>else</u> |
||
# 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 # |
# 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) |
||
) |
) |
||
<u>fi</u>; |
<u>fi</u>; |
||
Line 174: | Line 174: | ||
<u>format</u> real fmt = $g(-6,2)$; # width of 6, with no '+' sign, 2 decimals # |
<u>format</u> real fmt = $g(-6,2)$; # width of 6, with no '+' sign, 2 decimals # |
||
<u>proc</u> real matrix printf= (<u>format</u> real fmt, <u> |
<u>proc</u> real matrix printf= (<u>format</u> real fmt, <u>matrix</u> m)<u>void</u>:( |
||
<u>format</u> |
<u>format</u> vector fmt = $"("n(2⌈m-1)(f(real fmt)",")f(real fmt)")"$; |
||
<u>format</u> matrix fmt = $x"("n(⌈m-1)(f( |
<u>format</u> matrix fmt = $x"("n(⌈m-1)(f(vector fmt)","lxx)f(vector fmt)");"$; |
||
# finally print the result # |
# finally print the result # |
||
printf((matrix fmt,m)) |
printf((matrix fmt,m)) |
||
Line 182: | Line 182: | ||
# Some sample matrices to test # |
# Some sample matrices to test # |
||
<u> |
<u>matrix</u> a=((1, 1, 1, 1), # matrix A # |
||
(2, 4, 8, 16), |
(2, 4, 8, 16), |
||
(3, 9, 27, 81), |
(3, 9, 27, 81), |
||
(4, 16, 64, 256)); |
(4, 16, 64, 256)); |
||
<u> |
<u>matrix</u> b=(( 4 , -3 , 4/3, -1/4 ), # matrix B # |
||
(-13/3, 19/4, -7/3, 11/24), |
(-13/3, 19/4, -7/3, 11/24), |
||
( 3/2, -2 , 7/6, -1/4 ), |
( 3/2, -2 , 7/6, -1/4 ), |
||
( -1/6, 1/4, -1/6, 1/24)); |
( -1/6, 1/4, -1/6, 1/24)); |
||
<u> |
<u>matrix</u> c = a × b; # actual multiplication example of A x B # |
||
print((" A x B =",new line)); |
print((" A x B =",new line)); |
||
Line 199: | Line 199: | ||
exception index error: |
exception index error: |
||
putf(stand error, $x"Exception: index error."l$) |
putf(stand error, $x"Exception: index error."l$) |
||
⚫ | |||
⚫ | |||
=={{header|BASIC}}== |
=={{header|BASIC}}== |