Matrix multiplication: Difference between revisions

Line 37:

MODE FIELD = LONG REAL; # field type is LONG REAL #

INT default upb:=3;

MODE ~~VEC~~ = [default upb]FIELD;

MODE VECTOR = [default upb]FIELD;

MODE ~~MAT~~ = [default upb,default upb]FIELD;

MODE MATRIX = [default upb,default upb]FIELD;

# crude exception handling #

Line 44:

# define the vector/matrix operators #

OP * = (~~VEC~~ a,b)FIELD: ( # basically the dot product #

OP * = (VECTOR 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 * = (~~VEC~~ a, ~~MAT~~ b)~~VEC~~: ( # overload ~~vec~~ times matrix #

OP * = (VECTOR a, MATRIX b)VECTOR: ( # overload vector 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:

# this is the task portion #

OP * = (~~MAT~~ a, b)~~MAT~~: ( # overload matrix times matrix #

OP * = (MATRIX a, b)MATRIX: ( # 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 #

~~MAT~~ 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));

~~MAT~~ 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));

~~MAT~~ prod = a * b; # actual multiplication example of A x B #

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 #

PROC real matrix printf= (FORMAT real fmt, ~~MAT~~ m)VOID:(

PROC real matrix printf= (FORMAT real fmt, MATRIX m)VOID:(

FORMAT ~~vec~~ fmt = $"("n(2 UPB m-1)(f(real fmt)",")f(real fmt)")"$;

FORMAT vector fmt = $"("n(2 UPB m-1)(f(real fmt)",")f(real fmt)")"$;

FORMAT matrix fmt = $x"("n(UPB m-1)(f(~~vec~~ fmt)","lxx)f(~~vec~~ fmt)");"$;

FORMAT matrix fmt = $x"("n(UPB m-1)(f(vector fmt)","lxx)f(vector fmt)");"$;

# finally print the result #

printf((matrix fmt,m))

Line 102:

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

int default upb := 1;

int default upb := 3;

mode field = long real;

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

mode vector = [default upb]field;

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

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

# crude exception handling #

Line 113:

# define an operator to slice array into quarters #

op top = (~~mat~~ m)int: ( ⌊m + ⌈m ) %2,

op top = (matrix m)int: ( ⌊m + ⌈m ) %2,

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

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

left = (~~mat~~ m)int: ( ~~2 ⌊m~~ + ~~2 ⌈m~~ ) %2,

left = (matrix m)int: ( 2⌊m + 2⌈m ) %2,

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

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

left = (~~vec~~ v)int: ( ⌊v + ⌈v ) %2,

left = (vector v)int: ( ⌊v + ⌈v ) %2,

right = (~~vec~~ 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 #

op × = (~~vec~~ a, b)field: ( # dot product #

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

if ~~⌊a≠⌊b~~ ∨ ~~⌈a≠⌈b~~ then raise index error fi;

if (⌊a, ⌈a) ≠ (⌊b, ⌈b) then raise index error fi;

if ⌊a = ⌈a then

a[⌈a] × b[⌈b]

else

field begin, end;

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

[]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 #

for thread to ~~⌈thread schedule~~ do ~~thread~~ schedule[thread] od

for thread to ⌈schedule do schedule[thread] od

else

~~par~~ ( # run vector in parallel #

# PAR # ( # run vector in parallel #

~~thread~~ schedule[1], # assume parent CPU #

schedule[1], # assume parent CPU #

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

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

)

fi;

Line 143:

);

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

op × = (matrix a, b)matrix: # matrix multiply #

if ⌊a = ⌈a ~~∧ 2 ⌊b = 2 ⌈b~~ then

if (⌊a, 2⌊b) = (⌈a, 2⌈b) then

a[⌊a, ] × b[, ~~2 ⌈b~~] # dot product #

a[⌊a, ] × b[, 2⌈b] # dot product #

else

[⌈a, ~~2 ⌈b~~] field out;

[⌈a, 2⌈b] field out;

if 2 ~~⌊a≠⌊b~~ ∨ 2 ~~⌈a≠⌈b~~ then raise index error fi;

if (2⌊a, 2⌈a) ≠ (⌊b, ⌈b) then raise index error fi;

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

( true, # calculate top left corner #

Line 154:

( ⌊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 #

( 2⌊b ≠ 2⌈b, # calculate top right corner #

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

( 2 ⌊b ≠ 2 ⌈b ~~∧⌊a ≠ ⌈a~~, # calculate bottom right corner #

( (⌊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 #

for thread to ⌈schedule do (~~required →schedule~~[thread] | ~~thread →schedule~~[thread] ) od

for thread to ⌈schedule do (required→schedule[thread] | thread→schedule[thread] ) od

else

~~par~~ ( # run vector in parallel #

# PAR # ( # run vector in parallel #

~~thread →schedule~~[1], # thread is always required, and assume parent CPU #

thread→schedule[1], # thread is always required, and assume parent CPU #

( ~~required →schedule~~[4] | ↓idle cpus; ~~thread →schedule~~[4]; ↑idle cpus),

( 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[3] | ↓idle cpus; thread→schedule[3]; ↑idle cpus),

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

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

)

fi;

Line 174:

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

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

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

format ~~vec~~ fmt = $"("n(~~2 ⌈m~~-1)(f(real fmt)",")f(real fmt)")"$;

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

format matrix fmt = $x"("n(⌈m-1)(f(~~vec~~ fmt)","lxx)f(~~vec~~ fmt)");"$;

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

# finally print the result #

printf((matrix fmt,m))

Line 182:

# Some sample matrices to test #

~~mat~~ 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));

~~mat~~ 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));

~~mat~~ 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));

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 underline 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}}==