Matrix multiplication: Difference between revisions
m (→[[ALGOL 68]]) |
m (→[[ALGOL 68]]: turn parallel on) |
||
Line 134: | Line 134: | ||
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> schedule[thread] <u>od</u> |
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> schedule[thread] <u>od</u> |
||
<u>else</u> |
<u>else</u> |
||
<u>par</u>( # run vector in parallel # |
|||
schedule[1], # assume parent CPU # |
schedule[1], # assume parent CPU # |
||
( ↓idle cpus; schedule[2]; ↑idle cpus) |
( ↓idle cpus; schedule[2]; ↑idle cpus) |
||
Line 162: | Line 162: | ||
<u>for</u> thread <u>to</u> ⌈schedule <u>do</u> (required→schedule[thread] | thread→schedule[thread] ) <u>od</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> |
||
<u>par</u>( # 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), |
Revision as of 00:38, 22 April 2008
You are encouraged to solve this task according to the task description, using any language you may know.
Multiply two Matrices together, they can be of any dimension as long as the number of columns of the first matrix is equal to the number of rows of the second matrix
Ada
package Matrix_Ops is type Matrix is array(Natural range <>, Natural range <>) of Float; function "*" (Left, Right : Matrix) return Matrix; Dimension_Violation : exception; end Matrix_Ops;
package body Matrix_Ops is --------- -- "*" -- --------- function "*" (Left, Right : Matrix) return Matrix is Temp : Matrix(Left'Range(1), Right'Range(2)) := (Others =>(Others => 0.0)); begin if Left'Length(2) /= Right'Length(1) then raise Dimension_Violation; end if; for I in Left'range(1) loop for J in Right'range(2) loop for K in Left'range(2) loop Temp(I,J) := Temp(I,J) + Left(I, K)*Right(K, J); end loop; end loop; end loop; return Temp; end "*"; end Matrix_Ops;
ALGOL 68
An example of user defined Vector and Matrix Multiplication Operators:
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 # 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$)
Output:
Product of a and b: (( 1.00, -0.00, -0.00, -0.00), ( -0.00, 1.00, -0.00, -0.00), ( -0.00, -0.00, 1.00, -0.00), ( -0.00, -0.00, -0.00, 1.00));
Alternatively, for multicore CPUs, the following recursive algorithm can be used:
int default upb := 3; mode field = long real; 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; 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: ( ⌊m + ⌈m ) %2, bot = (matrix m)int: top m + 1, left = (matrix m)int: ( 2⌊m + 2⌈m ) %2, right = (matrix m)int: left m + 1, left = (vector v)int: ( ⌊v + ⌈v ) %2, right = (vector v)int: left v + 1; prio top = 8, bot = 8, left = 8, right = 8; # Operator priority - same as LWB & UPB # op × = (vector a, b)field: ( # dot product # if (⌊a, ⌈a) ≠ (⌊b, ⌈b) then raise index error fi; if ⌊a = ⌈a then a[⌈a] × b[⌈b] else field begin, end; []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 ⌈schedule do schedule[thread] od else par( # run vector in parallel # schedule[1], # assume parent CPU # ( ↓idle cpus; schedule[2]; ↑idle cpus) ) fi; begin+end fi ); op × = (matrix a, b)matrix: # matrix multiply # if (⌊a, 2⌊b) = (⌈a, 2⌈b) then a[⌊a, ] × b[, 2⌈b] # dot product # else [⌈a, 2⌈b] field out; if (2⌊a, 2⌈a) ≠ (⌊b, ⌈b) then raise index error fi; []struct(bool required, proc void thread) schedule = ( ( true, # calculate top left corner # void: out[:top a, :left b] := a[:top a, ] × b[, :left b]), ( ⌊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 # void: out[:top a, right b:] := a[:top a, ] × b[, right b:]), ( (⌊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 else 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 # ( required→schedule[3] | ↓idle cpus; thread→schedule[3]; ↑idle cpus), ( required→schedule[2] | ↓idle cpus; thread→schedule[2]; ↑idle cpus) ) fi; out fi; 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⌈m-1)(f(real fmt)",")f(real fmt)")"$; format matrix fmt = $x"("n(⌈m-1)(f(vector fmt)","lxx)f(vector fmt)");"$; # finally print the result # printf((matrix fmt,m)) ); # Some sample matrices to 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 c = a × b; # actual multiplication example of A x B # print((" A x B =",new line)); real matrix printf(real fmt, c) ∎ 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.
BASIC
Assume the matrices to be multiplied are a and b
IF (LEN(a,2) = LEN(b)) 'if valid dims n = LEN(a,2) m = LEN(a) p = LEN(b,2) DIM ans(0 TO m - 1, 0 TO p - 1) FOR i = 0 TO m - 1 FOR j = 0 TO p - 1 FOR k = 0 TO n - 1 ans(i, j) = ans(i, j) + (a(i, k) * b(k, j)) NEXT k, j, i 'print answer FOR i = 0 TO m - 1 FOR j = 0 TO p - 1 PRINT ans(i, j); NEXT j PRINT NEXT i ELSE PRINT "invalid dimensions" END IF
C
#include <stdio.h> #define dim 4 /* fixed length square matrices */ const int SLICE=0; /* coder hints */ typedef double field_t; /* field_t type is long float */ typedef field_t vec_t[dim]; typedef field_t *ref_vec_t; /* address of first element */ typedef vec_t matrix_t[dim]; typedef field_t *ref_matrix_t; /* address of first element */ typedef char* format; /* define the vector/matrix_t operators */ field_t v_times_v (vec_t a, vec_t b, int step_b){ /* basically the dot product if step_b==1*/ field_t result=0; for( int i=0; i<sizeof a; i++ ) result+= a[i]*b[i*step_b]; return result; } ref_vec_t v_eq_v_times_m(vec_t result, vec_t a, matrix_t b){ for( int j=0; j<sizeof b; j++ ) result[j]=v_times_v(a,&b[SLICE][j],sizeof b[SLICE] / sizeof (field_t)); return &result[SLICE]; } ref_matrix_t m_eq_m_times_m (matrix_t result, matrix_t a, matrix_t b){ for( int k=0; k<sizeof result; k++ ) v_eq_v_times_m(&result[k][SLICE],&a[k][SLICE],b); return &result[SLICE][SLICE]; } /* Some sample matrices to test */ matrix_t a={{1, 1, 1, 1}, /* matrix_t A */ {2, 4, 8, 16}, {3, 9, 27, 81}, {4, 16, 64, 256}}; matrix_t b={{ 4.0 , -3.0 , 4.0/3, -1.0/4 }, /* matrix_t B */ {-13.0/3, 19.0/4, -7.0/3, 11.0/24}, { 3.0/2, -2.0 , 7.0/6, -1.0/4 }, { -1.0/6, 1.0/4, -1.0/6, 1.0/24}}; int main(){ matrix_t prod; m_eq_m_times_m(prod,a,b); /* actual multiplication example of A x B */ #define field_fmt "%6.2f" /* width of 6, with no '+' sign, 2 decimals */ #define vec_fmt "{"field_fmt","field_fmt","field_fmt","field_fmt"}" #define matrix_fmt " {"vec_fmt",\n "vec_fmt",\n "vec_fmt",\n "vec_fmt"};" format result_fmt = " Product of a and b: \n"matrix_fmt"\n"; /* finally print the result */ vprintf(result_fmt,(void*)&prod); }
Output:
Product of a and b: {{ 1.00, 0.00, -0.00, -0.00}, { 0.00, 1.00, -0.00, -0.00}, { 0.00, 0.00, 1.00, -0.00}, { 0.00, 0.00, 0.00, 1.00}};
Common Lisp
(defun matrix-multiply (a b) (flet ((col (mat i) (mapcar #'(lambda (row) (elt row i)) mat)) (row (mat i) (elt mat i))) (loop for row from 0 below (length a) collect (loop for col from 0 below (length (row b 0)) collect (apply #'+ (mapcar #'* (row a row) (col b col))))))) ;; example use: (matrix-multiply '((1 2) (3 4)) '((-3 -8 3) (-2 1 4)))
(defun matrix-multiply (matrix1 matrix2) (mapcar (lambda (row) (apply #'mapcar (lambda (&rest column) (apply #'+ (mapcar #'* row column))) matrix2)) matrix1))
D
module mxmul ; import std.stdio ; import std.string ; bool isRectangular(T)(T[][] a){ if(a.length < 1 || a[0].length < 1) return false ; for(int i = 1 ; i < a.length; i++) if(a[i].length != a[i-1].length) return false ; return true ; } T[][] mmul(T=real)(T[][] lhs, T[][] rhs) { T[][] res ; if(isRectangular(lhs) && isRectangular(rhs) && lhs[0].length == rhs.length){ res = new T[][](lhs.length, rhs[0].length) ; for(int i = 0 ; i < lhs.length ; i++) for(int j = 0 ; j < rhs[0].length ; j++) { res[i][j] = 0 ; for(int k = 0 ; k < rhs.length ; k++) res[i][j] += lhs[i][k] * rhs[k][j] ; } } else throw new Exception("Mul. Error") ; return res ; } string toString(T=real)(T[][] a){ // for pretty print string[] s; foreach(e ; a) s ~= format("%8s", e)[1..$-1] ; return "\n<" ~ join(s,"\n ") ~ ">" ; } void main() { float[][] m = [[0.5,0,0,1],[0.0,0.5,0,0],[0.0,0,0.5,-1]] ; float[][] n = [[0.5,0,0],[0.0,0.5,0],[0.0,0,0.5],[1.0,0,-1]] ; writefln(" m = ", m.toString()) ; writefln(" n (m's transpose) = ", n.toString()) ; writefln(" m * n = ", m.mmul(n).toString()) ; writefln(" n * m = ", n.mmul(m).toString()) ; writefln("(n * m) * n = ", n.mmul(m).mmul(n).toString()) ; }
Fortran
In ISO Fortran 90 or later, use the SIZE and MATMUL intrinsic functions:
REAL, DIMENSION(N,M) :: A = RESHAPE( (/ (I, I=1, N*M) /), (/ N, M /) ) REAL, DIMENSION(M,K) :: B = RESHAPE( (/ (I, I=1, M*K) /), (/ M, K /) ) REAL, DIMENSION(SIZE(A,1), SIZE(B,2)) :: C ! C is an array whose first dimension (row) size is the same ! as A's first dimension size, and whose second dimension ! (column) size is the same as B's second dimension size. C = MATMUL( A, B )
Haskell
A somewhat inefficient version with lists (transpose is expensive):
import Data.List mmult :: Num a => [[a]] -> [[a]] -> [[a]] mmult a b = [ [ sum $ zipWith (*) ar bc | bc <- (transpose b) ] | ar <- a ] -- Example use: test = [[1, 2], [3, 4]] `mmult` [[-3, -8, 3], [-2, 1, 4]]
A more efficient version, based on arrays:
import Data.Array mmult :: (Ix i, Num a) => Array (i,i) a -> Array (i,i) a -> Array (i,i) a mmult x y | x1 /= y0 || x1' /= y0' = error "range mismatch" | otherwise = array ((x0,y1),(x0',y1')) l where ((x0,x1),(x0',x1')) = bounds x ((y0,y1),(y0',y1')) = bounds y ir = range (x0,x0') jr = range (y1,y1') kr = range (x1,x1') l = [((i,j), sum [x!(i,k) * y!(k,j) | k <- kr]) | i <- ir, j <- jr]
IDL
result = arr1 # arr2
J
x +/ .* y
where x and y are conformable arrays (trailing dimension of array x equals the leading dimension of array y). The notation is for a generalized inner product so that
x ~:/ .*. y NB. boolean inner product (~: is "not equal" (exclusive or) and *. is "and") x *./ .= y NB. which rows of x are the same as vector y? x + / .= y NB. number of places where each row of x equals vector y
etc.
Java
public static double[][] mult(double a[][], double b[][]){//a[m][n], b[n][p] if(a[0].length != b.length) return null; //invalid dims int n = a[0].length; int m = a.length; int p = b[0].length; double ans[][] = new double[m][p]; for(int i = 0;i < m;i++){ for(int j = 0;j < p;j++){ for(int k = 0;k < n;k++){ ans[i][j] += a[i][k] * b[k][j]; } } } return ans; }
Seed7
const type: matrix is array array float; const func matrix: (in matrix: left) * (in matrix: right) is func result var matrix: result is matrix.value; local var integer: i is 0; var integer: j is 0; var integer: k is 0; begin if length(left[1]) <> length(right) then raise RANGE_ERROR; else result := length(left) times length(right[1]) times 0.0; for i range 1 to length(left) do for j range 1 to length(right) do for k range 1 to length(left) do result[i][j] +:= left[i][k] * right[k][j]; end for; end for; end for; end if; end func;
SQL
CREATE TABLE a (x integer, y integer, e real); CREATE TABLE b (x integer, y integer, e real); -- test data -- A is a 2x2 matrix INSERT INTO a VALUES(0,0,1); INSERT INTO a VALUES(1,0,2); INSERT INTO a VALUES(0,1,3); INSERT INTO a VALUES(1,1,4); -- B is a 2x3 matrix INSERT INTO b VALUES(0,0,-3); INSERT INTO b VALUES(1,0,-8); INSERT INTO b VALUES(2,0,3); INSERT INTO b VALUES(0,1,-2); INSERT INTO b VALUES(1,1, 1); INSERT INTO b VALUES(2,1,4); -- C is 2x2 * 2x3 so will be a 2x3 matrix SELECT rhs.x, lhs.y, (SELECT sum(a.e*b.e) FROM a, b WHERE a.y = lhs.y AND b.x = rhs.x AND a.x = b.y) INTO TABLE c FROM a AS lhs, b AS rhs WHERE lhs.x = 0 AND rhs.y = 0;
TI-83 BASIC
Store your matrices in [A] and [B].
Disp [A]*[B]
An error will show if the matrices have invalid dimensions for multiplication.