In: Computer Science
Implement the Matrix Multiplication in Algol-68. I recommend using Algol 68 Genie for your interpreter/compiler.
An example of user defined Vector and Matrix Multiplication Operators using algol-68genie: MODE FIELD = LONG REAL; # field type is LONG REAL # INT default upb:=3; MODE VECTOR = [default upb]FIELD; MODE MATRIX = [default upb,default upb]FIELD; # crude exception handling # PROC VOID raise index error := VOID: GOTO exception index error; # define the vector/matrix operators # 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; FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD; result ); 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; FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD; result ); # this is the task portion # 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; FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD; result ); # Some sample matrices to test # test:( 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 # (-13/3, 19/4, -7/3, 11/24), ( 3/2, -2 , 7/6, -1/4 ), ( -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 # PROC real matrix printf= (FORMAT real fmt, MATRIX m)VOID:( FORMAT vector fmt = $"("n(2 UPB m-1)(f(real fmt)",")f(real fmt)")"$; FORMAT matrix fmt = $x"("n(UPB m-1)(f(vector fmt)","lxx)f(vector fmt)");"$; # finally print the result # printf((matrix fmt,m)) ); # finally print the result # print(("Product of a and b: ",new line)); real matrix printf(real fmt, prod) EXIT exception index error: putf(stand error, $x"Exception: index error."l$) )
Strassen matrix multiplication:
( # .nf Strassen's O(n^log2(7)) matrix multiplication algorithm L.Allison wacsvax UWA 26 June 1984 Dept. Computer Science, Monash University, Australia # OP * = ([,] INT a,b) [,]INT: IF INT all=1 UPB a; all = 1 LWB a THEN # 1x1 # a[1,1]*b[1,1] ELSE # nxn # INT half = all % 2; # i.e. all div 2 # [,]INT a11=a[1:half,1:half], # i.e. top left quarter of a # a12=a[1:half,half+1:all], # top right # a21=a[half+1:all,1:half], # etc. # a22=a[half+1:all,half+1:all], b11=b[1:half,1:half], b12=b[1:half,half+1:all], b21=b[half+1:all,1:half], b22=b[half+1:all,half+1:all]; LOC[1:all,1:all]INT c; # c will become a*b # REF[,]INT c11=c[1:half,1:half], c12=c[1:half,half+1:all], c21=c[half+1:all,1:half], c22=c[half+1:all,half+1:all]; [,]INT m1=(a12-a22)*(b21+b22), # NB. only 7 [,]*[,] # m2=(a11+a22)*(b11+b22), m3=(a11-a21)*(b11+b12), m4=(a11+a12)*b22, m5=a11*(b12-b22), m6=a22*(b21-b11), m7=(a21+a22)*b11; c11:=m1+m2-m4+m6; c12:=m4+m5; c21:=m6+m7; c22:=m2-m3+m5-m7; c # i.e. return c # FI; OP + = ([,]INT a,b)[,]INT: ( [1:1 UPB a, 1:2 UPB a]INT c; FOR i TO 1 UPB a DO FOR j TO 2 UPB a DO c[i,j]:=a[i,j]+b[i,j] OD OD; c ); OP - = ([,]INT a,b)[,]INT: ( [1:1 UPB a, 1:2 UPB a]INT c; FOR i TO 1 UPB a DO FOR j TO 2 UPB a DO c[i,j]:=a[i,j]-b[i,j] OD OD; c ); PROC show = ([,]INT a)VOID: ( print(newline); FOR i TO 1 UPB a DO FOR j TO 2 UPB a DO print( whole(a[i,j],6) ) OD; print(newline) OD; print(newline) ); [,]INT a=((1,2,3,4),(5,6,7,8),(9,10,11,12),(13,14,15,16)), # a little test # b=((17,18,19,20),(21,22,23,24),(25,26,27,28),(29,30,31,32)); show(a); show(b); show(a*b) )