Jump to content

Matrix multiplication: Difference between revisions

Matrix multiplication (view source)

Revision as of 00:29, 22 April 2008

78 bytes added ,  16 years ago
m
→‎[[ALGOL 68]]
m (→‎[[ALGOL 68]]: use the bold form of the if ~ then ~ else ~ fi syntax. - easier to read blocks.)
Line 37:
MODE FIELD = LONG REAL; # field type is LONG REAL #
INT default upb:=3;
MODE VECVECTOR = [default upb]FIELD;
MODE MATMATRIX = [default upb,default upb]FIELD;
# crude exception handling #
Line 44:
# define the vector/matrix operators #
OP * = (VECVECTOR a,b)FIELD: ( # basically the dot product #
FIELD result:=0;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
Line 51:
);
OP * = (VECVECTOR a, MATMATRIX b)VECVECTOR: ( # overload vecvector times matrix #
[2 LWB b:2 UPB b]FIELD result;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
Line 59:
<pre style="background-color:#ffe">
# this is the task portion #
OP * = (MATMATRIX a, b)MATMATRIX: ( # overload matrix times matrix #
[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;
Line 67:
</pre>
# Some sample matrices to test #
MATMATRIX a=((1, 1, 1, 1), # matrix A #
(2, 4, 8, 16),
(3, 9, 27, 81),
(4, 16, 64, 256));
MATMATRIX 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));
MATMATRIX prod = a * b; # actual multiplication example of A x B #
FORMAT real fmt = $g(-6,2)$; # width of 6, with no '+' sign, 2 decimals #
PROC real matrix printf= (FORMAT real fmt, MATMATRIX m)VOID:(
FORMAT vecvector fmt = $"("n(2 UPB m-1)(f(real fmt)",")f(real fmt)")"$;
FORMAT matrix fmt = $x"("n(UPB m-1)(f(vecvector fmt)","lxx)f(vecvector fmt)");"$;
# finally print the result #
printf((matrix fmt,m))
Line 102:
 
Alternatively, for multicore CPUs, the following recursive algorithm can be used:
<u>int</u> default upb := 13;
<u>mode</u> <u>field</u> = <u>long</u> <u>real</u>;
<u>mode</u> <u>vecvector</u> = [default upb]<u>field</u>;
<u>mode</u> <u>matmatrix</u> = [default upb, default upb]<u>field</u>;
# crude exception handling #
Line 113:
# define an operator to slice array into quarters #
<u>op</u> <u>top</u> = (<u>matmatrix</u> m)<u>int</u>: ( ⌊m + ⌈m ) %2,
<u>bot</u> = (<u>matmatrix</u> m)<u>int</u>: <u>top</u> m + 1,
<u>left</u> = (<u>matmatrix</u> m)<u>int</u>: ( 2 ⌊m2⌊m + 2 ⌈m2⌈m ) %2,
<u>right</u> = (<u>matmatrix</u> m)<u>int</u>: <u>left</u> m + 1,
<u>left</u> = (<u>vecvector</u> v)<u>int</u>: ( ⌊v + ⌈v ) %2,
<u>right</u> = (<u>vecvector</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>op</u> × = (<u>vecvector</u> a, b)<u>field</u>: ( # dot product #
<u>if</u> ⌊a≠⌊b(⌊a, ∨⌈a) ⌈a≠⌈b≠ (⌊b, ⌈b) <u>then</u> raise index error <u>fi</u>;
<u>if</u> ⌊a = ⌈a <u>then</u>
a[⌈a] × b[⌈b]
<u>else</u>
<u>field</u> begin, end;
[]<u>proc</u> <u>void</u> thread schedule=(
<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>if</u> <u>level</u> idle cpus = 0 <u>then</u> # use current CPU #
<u>for</u> thread <u>to</u> ⌈thread schedule⌈schedule <u>do</u> thread schedule[thread] <u>od</u>
<u>else</u>
<u>par</u># PAR # ( # run vector in parallel #
thread schedule[1], # assume parent CPU #
( ↓idle cpus; thread schedule[2]; ↑idle cpus)
)
<u>fi</u>;
Line 143:
);
<u>op</u> × = (<u>matmatrix</u> a, b)<u>matmatrix</u>: # matrix multiply #
<u>if</u> (⌊a, 2⌊b) = (⌈a, ∧ 2 ⌊b = 2 ⌈b2⌈b) <u>then</u>
a[⌊a, ] × b[, 2 ⌈b2⌈b] # dot product #
<u>else</u>
[⌈a, 2 ⌈b2⌈b] <u>field</u> out;
<u>if</u> 2(2⌊a, ⌊a≠⌊b2⌈a) ∨≠ 2(⌊b, ⌈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>true</u>, # calculate top left corner #
Line 154:
( ⌊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]),
( 2 ⌊b2⌊b ≠ 2 ⌈b2⌈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:]),
( 2(⌊a, ⌊b2⌊b) ≠ 2(⌈a, ⌈b2⌈b) ∧⌊a ≠ ⌈a, # 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>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> (required →schedulerequired→schedule[thread] | thread →schedulethread→schedule[thread] ) <u>od</u>
<u>else</u>
<u>par</u># PAR # ( # run vector in parallel #
thread →schedulethread→schedule[1], # thread is always required, and assume parent CPU #
( required →schedulerequired→schedule[4] | ↓idle cpus; thread →schedulethread→schedule[4]; ↑idle cpus),
# try to do opposite corners of matrix in parallel if CPUs are limited #
( required →schedulerequired→schedule[3] | ↓idle cpus; thread →schedulethread→schedule[3]; ↑idle cpus),
( required →schedulerequired→schedule[2] | ↓idle cpus; thread →schedulethread→schedule[2]; ↑idle cpus)
)
<u>fi</u>;
Line 174:
<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>matmatrix</u> m)<u>void</u>:(
<u>format</u> vecvector fmt = $"("n(2 ⌈m2⌈m-1)(f(real fmt)",")f(real fmt)")"$;
<u>format</u> matrix fmt = $x"("n(⌈m-1)(f(vecvector fmt)","lxx)f(vecvector fmt)");"$;
# finally print the result #
printf((matrix fmt,m))
Line 182:
# Some sample matrices to test #
<u>matmatrix</u> a=((1, 1, 1, 1), # matrix A #
(2, 4, 8, 16),
(3, 9, 27, 81),
(4, 16, 64, 256));
<u>matmatrix</u> 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));
<u>matmatrix</u> c = a × b; # actual multiplication example of A x B #
print((" A x B =",new line));
Line 199:
exception index error:
putf(stand error, $x"Exception: index error."l$)
Note: In the ALGOL 68 Revised Report the keywords, types and operators were permitted in a different typeface. Hence the underlingunderline of these identifiers above. Effectively this puts these identifiers into a different name space.▼
 
▲Note: In the ALGOL 68 Revised Report the keywords, types and operators were permitted in a different typeface. Hence the underling of these identifiers above. Effectively this puts these identifiers into a different name space.
 
=={{header|BASIC}}==
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/48912"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information