Matrix multiplication: Difference between revisions
(→Go: New section) |
(→{{header|Prolog}}: Added PureBasic) |
||
Line 985: | Line 985: | ||
mmult(M1, M2, M3) :- transpose(M2,MT), maplist(mm_helper(MT), M1, M3). |
mmult(M1, M2, M3) :- transpose(M2,MT), maplist(mm_helper(MT), M1, M3). |
||
mm_helper(M2, I1, M3) :- maplist(dot(I1), M2, M3).</lang> |
mm_helper(M2, I1, M3) :- maplist(dot(I1), M2, M3).</lang> |
||
=={{header|PureBasic}}== |
|||
Matrices represented as integer arrays with rows in the first dimension and columns in the second. |
|||
<lang PureBasic>Procedure multiplyMatrix(Array a(2), Array b(2), Array prd(2)) |
|||
Protected ar = ArraySize(a()) ;#rows for matrix a |
|||
Protected ac = ArraySize(a(), 2) ;#cols for matrix a |
|||
Protected br = ArraySize(b()) ;#rows for matrix b |
|||
Protected bc = ArraySize(b(), 2) ;#cols for matrix b |
|||
If ac = br |
|||
Dim prd(ar, bc) |
|||
Protected i, j, k |
|||
For i = 0 To ar |
|||
For j = 0 To bc |
|||
For k = 0 To br ;ac |
|||
prd(i, j) = prd(i, j) + (a(i, k) * b(k, j)) |
|||
Next |
|||
Next |
|||
Next |
|||
ProcedureReturn #True ;multiplication performed, product in prd() |
|||
Else |
|||
ProcedureReturn #False ;multiplication not performed, dimensions invalid |
|||
EndIf |
|||
EndProcedure</lang> |
|||
Additional code to demonstrate use. |
|||
<lang PureBasic>DataSection |
|||
Data.i 2,3 ;matrix a (#rows, #cols) |
|||
Data.i 1,2,3, 4,5,6 ;elements by row |
|||
Data.i 3,1 ;matrix b (#rows, #cols) |
|||
Data.i 1, 5, 9 ;elements by row |
|||
EndDataSection |
|||
Procedure displayMatrix(Array a(2), text.s) |
|||
Protected i, j |
|||
Protected columns = ArraySize(a(), 2), Rows = ArraySize(a(), 1) |
|||
PrintN(text + ": (" + Str(Rows + 1) + ", " + Str(columns + 1) + ")") |
|||
For i = 0 To Rows |
|||
For j = 0 To columns |
|||
Print(LSet(Str(a(i, j)), 4, " ")) |
|||
Next |
|||
PrintN("") |
|||
Next |
|||
PrintN("") |
|||
EndProcedure |
|||
Procedure loadMatrix(Array a(2)) |
|||
Protected Rows, columns, i, j |
|||
Read.i Rows |
|||
Read.i columns |
|||
Dim a(Rows - 1, columns - 1) |
|||
For i = 0 To Rows - 1 |
|||
For j = 0 To columns - 1 |
|||
Read.i a(i, j) |
|||
Next |
|||
Next |
|||
EndProcedure |
|||
Dim a(0,0) |
|||
Dim b(0,0) |
|||
Dim c(0,0) |
|||
If OpenConsole() |
|||
loadMatrix(a()) |
|||
displayMatrix(a(), "matrix a") |
|||
loadMatrix(b()) |
|||
displayMatrix(b(), "matrix b") |
|||
If multiplyMatrix(a(), b(), c()) |
|||
displayMatrix(c(), "product of a * b") |
|||
Else |
|||
PrintN("product of a * b is undefined") |
|||
EndIf |
|||
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit") |
|||
Input() |
|||
CloseConsole() |
|||
EndIf</lang> |
|||
Sample output: |
|||
<pre>matrix a: (2, 3) |
|||
1 2 3 |
|||
4 5 6 |
|||
matrix b: (3, 1) |
|||
1 |
|||
5 |
|||
9 |
|||
product of a * b: (2, 1) |
|||
38 |
|||
83</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 16:58, 28 May 2010
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 dimensions, so long as the number of columns of the first matrix is equal to the number of rows of the second matrix.
Ada
Ada has matrix multiplication predefined for any floating-point or complex type. The implementation is provided by the standard library packages Ada.Numerics.Generic_Real_Arrays and Ada.Numerics.Generic_Complex_Arrays correspondingly. The following example illustrates use of real matrix multiplication for the type Float: <lang ada>with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
procedure Matrix_Product is
procedure Put (X : Real_Matrix) is type Fixed is delta 0.01 range -100.0..100.0; begin for I in X'Range (1) loop for J in X'Range (2) loop Put (Fixed'Image (Fixed (X (I, J)))); end loop; New_Line; end loop; end Put; A : constant Real_Matrix := ( ( 1.0, 1.0, 1.0, 1.0), ( 2.0, 4.0, 8.0, 16.0), ( 3.0, 9.0, 27.0, 81.0), ( 4.0, 16.0, 64.0, 256.0) ); B : constant Real_Matrix := ( ( 4.0, -3.0, 4.0/3.0, -1.0/4.0 ), (-13.0/3.0, 19.0/4.0, -7.0/3.0, 11.0/24.0), ( 3.0/2.0, -2.0, 7.0/6.0, -1.0/4.0 ), ( -1.0/6.0, 1.0/4.0, -1.0/6.0, 1.0/24.0) );
begin
Put (A * B);
end Matrix_Product;</lang> Sample output:
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
The following code illustrates how matrix multiplication could be implemented from scratch: <lang ada>package Matrix_Ops is
type Matrix is array (Natural range <>, Natural range <>) of Float; function "*" (Left, Right : Matrix) return Matrix;
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 Constraint_Error; 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;</lang>
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 );
<lang algol68># 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 );</lang> # 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));
Parallel processing
Alternatively - for multicore CPUs - use the following reinvention of Strassen's O(n^log2(7)) recursive matrix multiplication algorithm:
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$)
APL
Matrix multiply in APL is just +.×. For example:
<lang apl> x ← +.×
A ← ↑A*¨⊂A←⍳4 ⍝ Same A as in other examples (1 1 1 1⍪ 2 4 8 16⍪ 3 9 27 81,[0.5] 4 16 64 256) B ← ⌹A ⍝ Matrix inverse of A 'F6.2' ⎕FMT A x 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</lang>
AutoHotkey
ahk discussion <lang autohotkey>Matrix("b"," ; rows separated by "," , 1 2 ; entries separated by space or tab , 2 3 , 3 0") MsgBox % "B`n`n" MatrixPrint(b) Matrix("c"," , 1 2 3 , 3 2 1") MsgBox % "C`n`n" MatrixPrint(c)
MatrixMul("a",b,c) MsgBox % "B * C`n`n" MatrixPrint(a)
MsgBox % MatrixMul("x",b,b)
Matrix(_a,_v) { ; Matrix structure: m_0_0 = #rows, m_0_1 = #columns, m_i_j = element[i,j], i,j > 0
Local _i, _j = 0 Loop Parse, _v, `, If (A_LoopField != "") { _i := 0, _j ++ Loop Parse, A_LoopField, %A_Space%%A_Tab% If (A_LoopField != "") _i++, %_a%_%_i%_%_j% := A_LoopField } %_a% := _a, %_a%_0_0 := _j, %_a%_0_1 := _i
} MatrixPrint(_a) {
Local _i = 0, _t Loop % %_a%_0_0 { _i++ Loop % %_a%_0_1 _t .= %_a%_%A_Index%_%_i% "`t" _t .= "`n" } Return _t
} MatrixMul(_a,_b,_c) {
Local _i = 0, _j, _k, _s If (%_b%_0_0 != %_c%_0_1) Return "ERROR: inner dimensions " %_b%_0_0 " != " %_c%_0_1 %_a% := _a, %_a%_0_0 := %_b%_0_0, %_a%_0_1 := %_c%_0_1 Loop % %_c%_0_1 { _i++, _j := 0 Loop % %_b%_0_0 { _j++, _k := _s := 0 Loop % %_b%_0_1 _k++, _s += %_b%_%_k%_%_j% * %_c%_%_i%_%_k% %_a%_%_i%_%_j% := _s } }
}</lang>
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
<lang 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 const 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);
}</lang> 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}};
C++
<lang cpp>#include <iostream>
- include <blitz/tinymat.h>
int main() {
using namespace blitz;
TinyMatrix<double,3,3> A, B, C;
A = 1, 2, 3, 4, 5, 6, 7, 8, 9;
B = 1, 0, 0, 0, 1, 0, 0, 0, 1;
C = product(A, B);
std::cout << C << std::endl;
}</lang> Output:
(3,3): [ 1 2 3 ] [ 4 5 6 ] [ 7 8 9 ]
Common Lisp
<lang 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))</lang>
D
<lang d>import std.stdio: writefln; import std.string: format, join;
T[][] matrixMul(T)(T[][] m1, T[][] m2) {
bool isRectangular(T[][] matrix) { foreach (row; matrix) if (row.length != matrix[0].length) return false; return true; }
T[][] result; if (isRectangular(m1) && isRectangular(m2) && m1[0].length == m2.length) { result = new T[][](m1.length, m2[0].length);
foreach (i, m1_row_i; m1) for (int j; j < m2[0].length; j++) { T aux = 0; foreach (k, m2_row_k; m2) aux += m1_row_i[k] * m2_row_k[j]; result[i][j] = aux; } } else throw new Exception("matrixMul Error"); return result;
}
string prettyPrint(T)(T[][] matrix) {
string[] result; foreach (row; matrix) result ~= format(row); return "[" ~ result.join(",\n ") ~ "]";
}
void main() {
float[][] a = [[1, 2], [3, 4], [3, 6]]; float[][] b = [[-3, -8, 3,], [-2, 1, 4]];
writefln("A = \n", prettyPrint(a)); writefln("\nB = \n", prettyPrint(b)); writefln("\nA * B = \n", prettyPrint(matrixMul(a, b)));
}</lang>
ELLA
Sample originally from ftp://ftp.dra.hmg.gb/pub/ella (a now dead link) - Public release.
Code for matrix multiplication hardware design verification: <lang ella>MAC ZIP = ([INT n]TYPE t: vector1 vector2) -> [n][2]t:
[INT k = 1..n](vector1[k], vector2[k]).
MAC TRANSPOSE = ([INT n][INT m]TYPE t: matrix) -> [m][n]t:
[INT i = 1..m] [INT j = 1..n] matrix[j][i].
MAC INNER_PRODUCT{FN * = [2]TYPE t -> TYPE s, FN + = [2]s -> s}
= ([INT n][2]t: vector) -> s: IF n = 1 THEN *vector[1] ELSE *vector[1] + INNER_PRODUCT {*,+} vector[2..n] FI.
MAC MATRIX_MULT {FN * = [2]TYPE t->TYPE s, FN + = [2]s->s} = ([INT n][INT m]t: matrix1, [m][INT p]t: matrix2) -> [n][p]s: BEGIN
LET transposed_matrix2 = TRANSPOSE matrix2.
OUTPUT [INT i = 1..n][INT j = 1..p]
INNER_PRODUCT{*,+}ZIP(matrix1[i],transposed_matrix2[j])
END.
TYPE element = NEW elt/(1..20),
product = NEW prd/(1..1200).
FN PLUS = (product: integer1 integer2) -> product:
ARITH integer1 + integer2.
FN MULT = (element: integer1 integer2) -> product:
ARITH integer1 * integer2.
FN MULT_234 = ([2][3]element:matrix1, [3][4]element:matrix2) ->
[2][4]product: MATRIX_MULT{MULT,PLUS}(matrix1, matrix2).
FN TEST = () -> [2][4]product: ( LET m1 = ((elt/2, elt/1, elt/1),
(elt/3, elt/6, elt/9)), m2 = ((elt/6, elt/1, elt/3, elt/4), (elt/9, elt/2, elt/8, elt/3), (elt/6, elt/4, elt/1, elt/2)). OUTPUT MULT_234 (m1, m2)
).
COM test: just displaysignal MOC</lang>
Factor
The built-in word m.
multiplies matrices:
( scratchpad ) USE: math.matrices { { 1 2 } { 3 4 } } { { -3 -8 3 } { -2 1 4 } } m. . { { -7 -6 11 } { -17 -20 25 } }
Forth
<lang forth>include fsl-util.f
3 3 float matrix A{{ ATemplate:3 3fread 1e 2e 3e 4e 5e 6e 7e 8e 9e 3 3 float matrix B{{ BTemplate:3 3fread 3e 3e 3e 2e 2e 2e 1e 1e 1e 3 3 float matrix C{{ \ result A{{ B{{ C{{ mat* C{{ }}print</lang>
Fortran
In ISO Fortran 90 or later, use the SIZE and MATMUL intrinsic functions: <lang fortran>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 )
print *, 'A' do i = 1, n
print *, a(i,:)
end do
print *, print *, 'B' do i = 1, m
print *, b(i,:)
end do
print *, print *, 'C = AB' do i = 1, n
print *, c(i,:)
end do</lang>
Go
<lang go>package main
import (
"fmt"
)
type Value float64 type Matrix [][]Value
func Multiply(m1, m2 Matrix) (m3 Matrix, ok bool) {
rows, cols, extra := len(m1), len(m2[0]), len(m2) if len(m1[0]) != extra { return nil, false } m3 = make(Matrix, rows) for i := 0; i < rows; i++ { m3[i] = make([]Value,cols) for j := 0; j < cols; j++ { for k := 0; k < extra; k++ { m3[i][j] += m1[i][k] * m2[k][j] } } } return m3, true
}
func (m Matrix) String() string {
rows := len(m) cols := len(m[0]) out := "[" for r := 0; r < rows; r++ { if r > 0 { out += ",\n " } out += "[ " for c := 0; c < cols; c++ { if c > 0 { out += ", " } out += fmt.Sprintf("%7.3f", m[r][c]) } out += " ]" } out += "]" return out
}
func main() {
A := Matrix{[]Value{1, 1, 1, 1}, []Value{2, 4, 8, 16}, []Value{3, 9, 27, 81}, []Value{4, 16, 64, 256}} B := Matrix{[]Value{ 4.0 , -3.0 , 4.0/3, -1.0/4 }, []Value{-13.0/3, 19.0/4, -7.0/3, 11.0/24}, []Value{ 3.0/2, -2.0 , 7.0/6, -1.0/4 }, []Value{ -1.0/6, 1.0/4, -1.0/6, 1.0/24}} P,ok := Multiply(A,B) if !ok { panic("Invalid dimensions") } fmt.Printf("Matrix A:\n%s\n\n", A) fmt.Printf("Matrix B:\n%s\n\n", B) fmt.Printf("Product of A and B:\n%s\n\n", P)
}</lang>
Output:
Matrix A: [[ 1.000, 1.000, 1.000, 1.000 ], [ 2.000, 4.000, 8.000, 16.000 ], [ 3.000, 9.000, 27.000, 81.000 ], [ 4.000, 16.000, 64.000, 256.000 ]] Matrix B: [[ 4.000, -3.000, 1.333, -0.250 ], [ -4.333, 4.750, -2.333, 0.458 ], [ 1.500, -2.000, 1.167, -0.250 ], [ -0.167, 0.250, -0.167, 0.042 ]] Product of A and B: [[ 1.000, 0.000, -0.000, -0.000 ], [ 0.000, 1.000, -0.000, -0.000 ], [ 0.000, 0.000, 1.000, -0.000 ], [ 0.000, 0.000, 0.000, 1.000 ]]
Haskell
A somewhat inefficient version with lists (transpose is expensive):
<lang haskell>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]]</lang>
A more efficient version, based on arrays:
<lang haskell>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]</lang>
HicEst
<lang hicest>REAL :: m=4, n=2, p=3, a(m,n), b(n,p), res(m,p)
a = $ ! initialize to 1, 2, ..., m*n b = $ ! initialize to 1, 2, ..., n*p
res = 0 DO i = 1, m
DO j = 1, p DO k = 1, n res(i,j) = res(i,j) + a(i,k) * b(k,j) ENDDO ENDDO
ENDDO
DLG(DefWidth=4, Text=a, Text=b,Y=0, Text=res,Y=0)</lang> <lang hicest>a b res 1 2 1 2 3 9 12 15 3 4 4 5 6 19 26 33 5 6 29 40 51 7 8 39 54 69 </lang>
IDL
<lang idl>result = arr1 # arr2</lang>
J
Matrix multiply in J is just +/ .*
. For example:
<lang j> mp =: +/ .* NB. Matrix product
A =: ^/~>:i. 4 NB. Same A as in other examples (1 1 1 1, 2 4 8 16, 3 9 27 81,:4 16 64 256) B =: %.A NB. Matrix inverse of A '6.2' 8!:2 A mp 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</lang> The notation is for a generalized inner product so that <lang j>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</lang> etc.
The general inner product extends to multidimensional arrays, requiring only that x and y be conformable (trailing dimension of array x equals the leading dimension of array y). For example, the matrix multiplication of two dimensional arrays requires x to have the same numbers of rows as y has columns, as you would expect.
Java
<lang java>public static double[][] mult(double a[][], double b[][]){//a[m][n], b[n][p]
if(a.length == 0) return new double[0][0]; 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;
}</lang>
JavaScript
for the print()
function
Extends Matrix Transpose#JavaScript <lang javascript>// returns a new matrix Matrix.prototype.mult = function(other) {
if (this.width != other.height) { throw "error: incompatible sizes"; }
var result = []; for (var i = 0; i < this.height; i++) { result[i] = []; for (var j = 0; j < other.width; j++) { var sum = 0; for (var k = 0; k < this.width; k++) { sum += this.mtx[i][k] * other.mtx[k][j]; } result[i][j] = sum; } } return new Matrix(result);
}
var a = new Matrix([[1,2],[3,4]]) var b = new Matrix([[-3,-8,3],[-2,1,4]]); print(a.mult(b));</lang> output
-7,-6,11 -17,-20,25
Mathematica
<lang mathematica>M1 = {{1, 2},
{3, 4}, {5, 6}, {7, 8}}
M2 = {{1, 2, 3},
{4, 5, 6}}
M = M1.M2</lang>
Or without the variables:
<lang mathematica>{{1, 2}, {3, 4}, {5, 6}, {7, 8}}.{{1, 2, 3}, {4, 5, 6}}</lang>
The result is: <lang mathematica>{{9, 12, 15}, {19, 26, 33}, {29, 40, 51}, {39, 54, 69}}</lang>
MATLAB
<lang Matlab>function [output] = matrixmultiplication(matrixA, matrixB)
output = matrixA*matrixB;</lang>
Nial
<lang nial>|A := 4 4 reshape 1 1 1 1 2 4 8 16 3 9 27 81 4 16 64 256 =1 1 1 1 =2 4 8 16 =3 9 27 81 =4 16 64 256 |B := inverse A
|A innerproduct B =1. 0. 8.3e-17 -2.9e-16 =1.3e-15 1. -4.4e-16 -3.3e-16 =0. 0. 1. 4.4e-16 =0. 0. 0. 1.</lang>
OCaml
This version works on arrays of arrays of ints: <lang ocaml>let matrix_multiply x y =
let x0 = Array.length x and y0 = Array.length y in let y1 = if y0 = 0 then 0 else Array.length y.(0) in let z = Array.make_matrix x0 y1 0 in for i = 0 to x0-1 do for j = 0 to y1-1 do for k = 0 to y0-1 do z.(i).(j) <- z.(i).(j) + x.(i).(k) * y.(k).(j) done done done; z</lang>
# matrix_multiply [|[|1;2|];[|3;4|]|] [|[|-3;-8;3|];[|-2;1;4|]|];; - : int array array = [|[|-7; -6; 11|]; [|-17; -20; 25|]|]
This version works on lists of lists of ints: <lang ocaml>(* equivalent to (apply map ...) *) let rec mapn f lists =
assert (lists <> []); if List.mem [] lists then [] else f (List.map List.hd lists) :: mapn f (List.map List.tl lists)
let matrix_multiply m1 m2 =
List.map (fun row -> mapn (fun column -> List.fold_left (+) 0 (List.map2 ( * ) row column)) m2) m1</lang>
# matrix_multiply [[1;2];[3;4]] [[-3;-8;3];[-2;1;4]];; - : int list list = [[-7; -6; 11]; [-17; -20; 25]]
Octave
<lang octave>a = zeros(4); % prepare the matrix % 1 1 1 1 % 2 4 8 16 % 3 9 27 81 % 4 16 64 256 for i = 1:4
for j = 1:4 a(i, j) = i^j; endfor
endfor b = inverse(a); a * b</lang>
Perl
For most applications involving extensive matrix arithmetic, using the CPAN module called "PDL" (that stands for "Perl Data Language") would probably be the easiest and most efficient approach. That said, here's an implementation of matrix multiplication in plain Perl.
<lang perl>sub mmult
{our @a; local *a = shift; our @b; local *b = shift; my @p = []; my $rows = @a; my $cols = @{$b[0]}; my $n = @b - 1; for (my $r = 0 ; $r < $rows ; ++$r) {for (my $c = 0 ; $c < $cols ; ++$c) {$p[$r][$c] += $a[$r][$_] * $b[$_][$c] foreach 0 .. $n;}} return [@p];}</lang>
This function takes two references to arrays of arrays and returns the product as a reference to a new anonymous array of arrays.
PicoLisp
<lang PicoLisp>(de matMul (Mat1 Mat2)
(mapcar '((Row) (apply mapcar Mat2 '(@ (apply + (mapcar * Row (rest)))) ) ) Mat1 ) )
(matMul
'((1 2 3) (4 5 6)) '((6 -1) (3 2) (0 -3)) )</lang>
Output:
-> ((12 -6) (39 -12))
PL/I
<lang PL/I> /* Matrix multiplication of A by B, yielding C */ MMULT: procedure (a, b, c);
declare (a, b, c)(*,*) float controlled; declare (i, j, m, n, p) fixed binary;
if hbound(a,2) ^= hbound(b,1) then do; put skip list ('Matrices are incompatible for matrix multiplication'); signal error; end;
m = hbound(a, 1); p = hbound(b, 2); if allocation(c) > 0 then free c;
allocate c(m,p);
do i = 1 to m; do j = 1 to p; c(i,j) = sum(a(i,*) * b(*,j) ); end; end;
end MMULT; </lang>
Pop11
<lang pop11>define matmul(a, b) -> c;
lvars ba = boundslist(a), bb = boundslist(b); lvars i, i0 = ba(1), i1 = ba(2); lvars j, j0 = bb(1), j1 = bb(2); lvars k, k0 = bb(3), k1 = bb(4); if length(ba) /= 4 then throw([need_2d_array ^a]) endif; if length(bb) /= 4 then throw([need_2d_array ^b]) endif; if ba(3) /= j0 or ba(4) /= j1 then throw([dimensions_do_not_match ^a ^b]); endif; newarray([^i0 ^i1 ^k0 ^k1], 0) -> c; for i from i0 to i1 do for k from k0 to k1 do for j from j0 to j1 do c(i, k) + a(i, j)*b(j, k) -> c(i, k); endfor; endfor; endfor;
enddefine;</lang>
Prolog
<lang prolog>% SWI-Prolog has transpose/2 in its clpfd library
- - use_module(library(clpfd)).
% N is the dot product of lists V1 and V2. dot(V1, V2, N) :- maplist(product,V1,V2,P), sumlist(P,N). product(N1,N2,N3) :- N3 is N1*N2.
% Matrix multiplication with matrices represented % as lists of lists. M3 is the product of M1 and M2 mmult(M1, M2, M3) :- transpose(M2,MT), maplist(mm_helper(MT), M1, M3). mm_helper(M2, I1, M3) :- maplist(dot(I1), M2, M3).</lang>
PureBasic
Matrices represented as integer arrays with rows in the first dimension and columns in the second. <lang PureBasic>Procedure multiplyMatrix(Array a(2), Array b(2), Array prd(2))
Protected ar = ArraySize(a()) ;#rows for matrix a Protected ac = ArraySize(a(), 2) ;#cols for matrix a Protected br = ArraySize(b()) ;#rows for matrix b Protected bc = ArraySize(b(), 2) ;#cols for matrix b If ac = br Dim prd(ar, bc) Protected i, j, k For i = 0 To ar For j = 0 To bc For k = 0 To br ;ac prd(i, j) = prd(i, j) + (a(i, k) * b(k, j)) Next Next Next ProcedureReturn #True ;multiplication performed, product in prd() Else ProcedureReturn #False ;multiplication not performed, dimensions invalid EndIf
EndProcedure</lang> Additional code to demonstrate use. <lang PureBasic>DataSection
Data.i 2,3 ;matrix a (#rows, #cols) Data.i 1,2,3, 4,5,6 ;elements by row Data.i 3,1 ;matrix b (#rows, #cols) Data.i 1, 5, 9 ;elements by row
EndDataSection
Procedure displayMatrix(Array a(2), text.s)
Protected i, j Protected columns = ArraySize(a(), 2), Rows = ArraySize(a(), 1) PrintN(text + ": (" + Str(Rows + 1) + ", " + Str(columns + 1) + ")") For i = 0 To Rows For j = 0 To columns Print(LSet(Str(a(i, j)), 4, " ")) Next PrintN("") Next PrintN("")
EndProcedure
Procedure loadMatrix(Array a(2))
Protected Rows, columns, i, j Read.i Rows Read.i columns Dim a(Rows - 1, columns - 1) For i = 0 To Rows - 1 For j = 0 To columns - 1 Read.i a(i, j) Next Next
EndProcedure
Dim a(0,0) Dim b(0,0) Dim c(0,0)
If OpenConsole()
loadMatrix(a()) displayMatrix(a(), "matrix a") loadMatrix(b()) displayMatrix(b(), "matrix b") If multiplyMatrix(a(), b(), c()) displayMatrix(c(), "product of a * b") Else PrintN("product of a * b is undefined") EndIf Print(#CRLF$ + #CRLF$ + "Press ENTER to exit") Input() CloseConsole()
EndIf</lang> Sample output:
matrix a: (2, 3) 1 2 3 4 5 6 matrix b: (3, 1) 1 5 9 product of a * b: (2, 1) 38 83
Python
<lang python>a=((1, 1, 1, 1), # matrix A #
(2, 4, 8, 16), (3, 9, 27, 81), (4, 16, 64, 256))
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.))
def MatrixMul( mtx_a, mtx_b):
tpos_b = zip( *mtx_b) rtn = [[ sum( ea*eb for ea,eb in zip(a,b)) for b in tpos_b] for a in mtx_a] return rtn
v = MatrixMul( a, b )
print 'v = (' for r in v:
print '[', for val in r: print '%8.2f '%val, print ']'
print ')'
u = MatrixMul(b,a)
print 'u = ' for r in u:
print '[', for val in r: print '%8.2f '%val, print ']'
print ')'</lang>
Another one,
<lang python>from operator import mul
def matrixMul(m1, m2):
return map( lambda row: map( lambda *column: sum(map(mul, row, column)), *m2), m1)</lang>
R
<lang r>a %*% b</lang>
Ruby
Using
<lang ruby>require 'matrix'
Matrix[[1, 2],
[3, 4]] * Matrix[[-3, -8, 3], [-2, 1, 4]]</lang>
Output:
Matrix[[-7, -6, 11], [-17, -20, 25]]
Version for lists:
<lang ruby>def matrix_mult(a, b)
a.map do |ar| b.transpose.map do |bc| ar.zip(bc).map {|x,y| x*y}.inject {|z,w| z+w} end end
end</lang>
Scala
Assuming an array of arrays representation:
<lang scala>def mult[A](a: Array[Array[A]], b: Array[Array[A]])(implicit n: Numeric[A]) = {
import n._ for (row <- a) yield for(col <- b.transpose) yield row zip col map Function.tupled(_*_) reduceLeft (_+_)
}</lang>
For any subclass of Seq
(which does not include Java-specific arrays):
<lang scala>def mult[A, CC[X] <: Seq[X], DD[Y] <: Seq[Y]](a: CC[DD[A]], b: CC[DD[A]]) (implicit n: Numeric[A]): CC[DD[A]] = {
import n._ for (row <- a) yield for(col <- b.transpose) yield row zip col map Function.tupled(_*_) reduceLeft (_+_)
}</lang>
Examples:
scala> Array(Array(1, 2), Array(3, 4)) res0: Array[Array[Int]] = Array(Array(1, 2), Array(3, 4)) scala> Array(Array(-3, -8, 3), Array(-2, 1, 4)) res1: Array[Array[Int]] = Array(Array(-3, -8, 3), Array(-2, 1, 4)) scala> mult(res0, res1) res2: Array[scala.collection.mutable.GenericArray[Int]] = Array(GenericArray(-7, -6, 11), GenericArray(-17, -20, 25)) scala> res0.map(_.toList).toList res5: List[List[Int]] = List(List(1, 2), List(3, 4)) scala> res1.map(_.toList).toList res6: List[List[Int]] = List(List(-3, -8, 3), List(-2, 1, 4)) scala> mult(res5, res6) res7: Seq[Seq[Int]] = List(List(-7, -6, 11), List(-17, -20, 25))
A fully generic multiplication that returns the same collection as received is possible, but much more verbose.
Scheme
This version works on lists of lists: <lang scheme>(define (matrix-multiply matrix1 matrix2)
(map (lambda (row) (apply map (lambda column (apply + (map * row column))) matrix2)) matrix1))</lang>
> (matrix-multiply '((1 2) (3 4)) '((-3 -8 3) (-2 1 4))) ((-7 -6 11) (-17 -20 25))
Seed7
<lang 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; var float: accumulator is 0.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 accumulator := 0.0; for k range 1 to length(left) do accumulator +:= left[i][k] * right[k][j]; end for; result[i][j] := accumulator; end for; end for; end if; end func;</lang>
Original source: [1]
SQL
<lang 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;</lang>
Tcl
<lang tcl>package require Tcl 8.5 namespace path ::tcl::mathop proc matrix_multiply {a b} {
lassign [size $a] a_rows a_cols lassign [size $b] b_rows b_cols if {$a_cols != $b_rows} { error "incompatible sizes: a($a_rows, $a_cols), b($b_rows, $b_cols)" } set temp [lrepeat $a_rows [lrepeat $b_cols 0]] for {set i 0} {$i < $a_rows} {incr i} { for {set j 0} {$j < $b_cols} {incr j} { set sum 0 for {set k 0} {$k < $a_cols} {incr k} { set sum [+ $sum [* [lindex $a $i $k] [lindex $b $k $j]]] } lset temp $i $j $sum } } return $temp
}</lang>
Using the print_matrix
procedure defined in Matrix Transpose#Tcl
% print_matrix [matrix_multiply {{1 2} {3 4}} {{-3 -8 3} {-2 1 4}}] -7 -6 11 -17 -20 25
TI-83 BASIC
Store your matrices in [A] and [B]. <lang ti83b>Disp [A]*[B]</lang> An error will show if the matrices have invalid dimensions for multiplication.
TI-89 BASIC
<lang ti89b>[1,2; 3,4; 5,6; 7,8] → m1 [1,2,3; 4,5,6] → m2 m1 * m2</lang>
Or without the variables:
<lang ti89b>[1,2; 3,4; 5,6; 7,8] * [1,2,3; 4,5,6]</lang>
The result (without prettyprinting) is:
<lang ti89b>[[9,12,15][19,26,33][29,40,51][39,54,69]]</lang>
Ursala
There is a library function for matrix multiplication of IEEE double precision floating point numbers. This example shows how to define and use a matrix multiplication function over any chosen field given only the relevant product and sum functions, in this case for the built in rational number type.
<lang Ursala>#import rat
a =
<
<1/1, 1/1, 1/1, 1/1>, <2/1, 4/1, 8/1, 16/1>, <3/1, 9/1, 27/1, 81/1>, <4/1, 16/1, 64/1, 256/1>>
b =
<
< 4/1, -3/1, 4/3, -1/4>, <-13/3, 19/4, -7/3, 11/24>, < 3/2, -2/1, 7/6, -1/4>, < -1/6, 1/4, -1/6, 1/24>>
mmult = *rK7lD *rlD sum:-0.+ product*p
- cast %qLL
test = mmult(a,b)</lang> output:
< <1/1,0/1,0/1,0/1>, <0/1,1/1,0/1,0/1>, <0/1,0/1,1/1,0/1>, <0/1,0/1,0/1,1/1>>
ZPL
<lang ZPL> program matmultSUMMA;
prototype GetSingleDim(infile:file):integer; prototype GetInnerDim(infile1:file; infile2:file):integer;
config var
Afilename: string = ""; Bfilename: string = "";
Afile: file = open(Afilename,file_read); Bfile: file = open(Bfilename,file_read);
default_size:integer = 4; m:integer = GetSingleDim(Afile); n:integer = GetInnerDim(Afile,Bfile); p:integer = GetSingleDim(Bfile);
iters: integer = 1;
printinput: boolean = false; verbose: boolean = true; dotiming: boolean = false;
region
RA = [1..m,1..n]; RB = [1..n,1..p]; RC = [1..m,1..p]; FCol = [1..m,*]; FRow = [*,1..p];
var
A : [RA] double; B : [RB] double; C : [RC] double; Aflood : [FCol] double; Bflood : [FRow] double;
procedure ReadA(); var step:double; [RA] begin
if (Afile != znull) then read(Afile,A); else step := 1.0/(m*n); A := ((Index1-1)*n + Index2)*step + 1.0; end; end;
procedure ReadB();
var step:double;
[RB] begin
if (Bfile != znull) then read(Bfile,B); else step := 1.0/(n*p); B := ((Index1-1)*p + Index2)*step + 1.0; end; end;
procedure matmultSUMMA();
var
i: integer; it: integer; runtime: double;
[RC] begin
ReadA(); ReadB();
if (printinput) then [RA] writeln("A is:\n",A); [RB] writeln("B is:\n",B); end;
ResetTimer();
for it := 1 to iters do C := 0.0; -- zero C for i := 1 to n do [FCol] Aflood := >>[,i] A; -- flood A col [FRow] Bflood := >>[i,] B; -- flood B row
C += (Aflood * Bflood); -- multiply end; end;
runtime := CheckTimer();
if (verbose) then writeln("C is:\n",C); end;
if (dotiming) then writeln("total runtime = %12.6f":runtime); writeln("actual runtime = %12.6f":runtime/iters); end; end;
procedure GetSingleDim(infile:file):integer;
var dim:integer;
begin
if (infile != znull) then read(infile,dim); else dim := default_size; end; return dim;
end;
procedure GetInnerDim(infile1:file; infile2:file):integer;
var
col:integer; row:integer; retval:integer;
begin
retval := -1; if (infile1 != znull) then read(infile1,col); retval := col; end; if (infile2 != znull) then read(infile2,row); if (retval = -1) then retval := row; else if (row != col) then halt("ERROR: Inner dimensions don't match"); end; end; end; if (retval = -1) then retval := default_size; end; return retval;
end; </lang>