Matrix multiplication: Difference between revisions
m →{{header|REXX}}: added comment, changed sep character, changed indentation of DO loops. -- ~~~~ |
m →{{header|REXX}}: fixed up a comment. -- ~~~~ |
||
Line 2,328: | Line 2,328: | ||
end |
end |
||
Brows=r-1; Bcols=c-1 |
Brows=r-1; Bcols=c-1 |
||
⚫ | |||
⚫ | |||
do i =1 for Arows /*multiply matrix A & B ──► C */ |
do i =1 for Arows /*multiply matrix A & B ──► C */ |
||
do j =1 for Bcols |
do j =1 for Bcols |
Revision as of 00:03, 28 July 2012
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: <lang algol68>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$)
)</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));
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
BBC BASIC
BBC BASIC has built-in matrix multiplication: <lang bbcbasic> DIM matrix1(3,1), matrix2(1,2), product(3,2)
matrix1() = 1, 2, \ \ 3, 4, \ \ 5, 6, \ \ 7, 8 matrix2() = 1, 2, 3, \ \ 4, 5, 6 product() = matrix1() . matrix2()
FOR row% = 0 TO DIM(product(),1) FOR col% = 0 TO DIM(product(),2) PRINT product(row%,col%),; NEXT PRINT NEXT
</lang> Output:
9 12 15 19 26 33 29 40 51 39 54 69
C
For performance critical work involving matrices, especially large or sparse ones, always consider using an established library such as BLAS first. <lang c>#include <stdio.h>
- include <stdlib.h>
/* Make the data structure self-contained. Element at row i and col j
is x[i * w + j]. More often than not, though, you might want to represent a matrix some other way */
typedef struct { int h, w; double *x;} matrix_t, *matrix;
inline double dot(double *a, double *b, int len, int step) { double r = 0; while (len--) { r += *a++ * *b; b += step; } return r; }
matrix mat_new(int h, int w) { matrix r = malloc(sizeof(matrix) + sizeof(double) * w * h); r->h = h, r->w = w; r->x = (double*)(r + 1); return r; }
matrix mat_mul(matrix a, matrix b) { matrix r; double *p, *pa; int i, j; if (a->w != b->h) return 0;
r = mat_new(a->h, b->w); p = r->x; for (pa = a->x, i = 0; i < a->h; i++, pa += a->w) for (j = 0; j < b->w; j++) *p++ = dot(pa, b->x + j, a->w, b->w); return r; }
void mat_show(matrix a) { int i, j; double *p = a->x; for (i = 0; i < a->h; i++, putchar('\n')) for (j = 0; j < a->w; j++) printf("\t%7.3f", *p++); putchar('\n'); }
int main() { double da[] = { 1, 1, 1, 1, 2, 4, 8, 16, 3, 9, 27, 81, 4,16, 64, 256 }; double db[] = { 4.0, -3.0, 4.0/3, -13.0/3, 19.0/4, -7.0/3, 3.0/2, -2.0, 7.0/6, -1.0/6, 1.0/4, -1.0/6};
matrix_t a = { 4, 4, da }, b = { 4, 3, db }; matrix c = mat_mul(&a, &b);
/* mat_show(&a), mat_show(&b); */ mat_show(c); /* free(c) */ return 0; }</lang>
C#
This code should work with any version of the .NET Framework and C# language
<lang csharp>public class Matrix { int n; int m; double[,] a;
public Matrix(int n, int m) { if (n <= 0 || m <= 0) throw new ArgumentException("Matrix dimensions must be positive"); this.n = n; this.m = m; a = new double[n, m]; }
//indices start from one public double this[int i, int j] { get { return a[i - 1, j - 1]; } set { a[i - 1, j - 1] = value; } }
public int N { get { return n; } } public int M { get { return m; } }
public static Matrix operator*(Matrix _a, Matrix b) { int n = _a.N; int m = b.M; int l = _a.M; if (l != b.N) throw new ArgumentException("Illegal matrix dimensions for multiplication. _a.M must be equal b.N"); Matrix result = new Matrix(_a.N, b.M); for(int i = 0; i < n; i++) for (int j = 0; j < m; j++) { double sum = 0.0; for (int k = 0; k < l; k++) sum += _a.a[i, k]*b.a[k, j]; result.a[i, j] = sum; } return result; } }</lang>
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 ]
Generic solution
main.cpp <lang cpp>
- include <iostream>
- include "matrix.h"
- if !defined(ARRAY_SIZE)
#define ARRAY_SIZE(x) (sizeof((x)) / sizeof((x)[0]))
- endif
int main() {
int am[2][3] = { {1,2,3}, {4,5,6}, }; int bm[3][2] = { {1,2}, {3,4}, {5,6} };
Matrix<int> a(ARRAY_SIZE(am), ARRAY_SIZE(am[0]), am[0], ARRAY_SIZE(am)*ARRAY_SIZE(am[0])); Matrix<int> b(ARRAY_SIZE(bm), ARRAY_SIZE(bm[0]), bm[0], ARRAY_SIZE(bm)*ARRAY_SIZE(bm[0])); Matrix<int> c;
try { c = a * b; for (unsigned int i = 0; i < c.rowNum(); i++) { for (unsigned int j = 0; j < c.colNum(); j++) { std::cout << c[i][j] << " "; } std::cout << std::endl; } } catch (MatrixException& e) { std::cerr << e.message() << std::endl; return e.errorCode(); }
} /* main() */ </lang>
matrix.h <lang cpp>
- ifndef _MATRIX_H
- define _MATRIX_H
- include <sstream>
- include <string>
- include <vector>
- define MATRIX_ERROR_CODE_COUNT 5
- define MATRIX_ERR_UNDEFINED "1 Undefined exception!"
- define MATRIX_ERR_WRONG_ROW_INDEX "2 The row index is out of range."
- define MATRIX_ERR_MUL_ROW_AND_COL_NOT_EQUAL "3 The row number of second matrix must be equal with the column number of first matrix!"
- define MATRIX_ERR_MUL_ROW_AND_COL_BE_GREATER_THAN_ZERO "4 The number of rows and columns must be greater than zero!"
- define MATRIX_ERR_TOO_FEW_DATA "5 Too few data in matrix."
class MatrixException { private:
std::string message_; int errorCode_;
public:
MatrixException(std::string message = MATRIX_ERR_UNDEFINED);
inline std::string message() { return message_; };
inline int errorCode() { return errorCode_; };
};
MatrixException::MatrixException(std::string message) {
errorCode_ = MATRIX_ERROR_CODE_COUNT + 1; std::stringstream ss(message); ss >> errorCode_; if (errorCode_ < 1) { errorCode_ = MATRIX_ERROR_CODE_COUNT + 1; } std::string::size_type pos = message.find(' '); if (errorCode_ <= MATRIX_ERROR_CODE_COUNT && pos != std::string::npos) { message_ = message.substr(pos + 1); } else { message_ = message + " (This an unknown and unsupported exception!)"; }
}
/**
* Generic class for matrices. */
template <class T> class Matrix { private:
std::vector<T> v; // the data of matrix unsigned int m; // the number of rows unsigned int n; // the number of columns
protected:
virtual void clear() { v.clear(); m = n = 0; }
public:
Matrix() { clear(); } Matrix(unsigned int, unsigned int, T* = 0, unsigned int = 0); Matrix(unsigned int, unsigned int, const std::vector<T>&);
virtual ~Matrix() { clear(); } Matrix& operator=(const Matrix&); std::vector<T> operator[](unsigned int) const; Matrix operator*(const Matrix&);
inline unsigned int rowNum() const { return m; }
inline unsigned int colNum() const { return n; }
inline unsigned int size() const { return v.size(); }
inline void add(const T& t) { v.push_back(t); }
};
template <class T> Matrix<T>::Matrix(unsigned int row, unsigned int col, T* data, unsigned int dataLength) {
clear(); if (row > 0 && col > 0) { m = row; n = col; unsigned int mxn = m * n; if (dataLength && data) { for (unsigned int i = 0; i < dataLength && i < mxn; i++) { v.push_back(data[i]); } } }
}
template <class T> Matrix<T>::Matrix(unsigned int row, unsigned int col, const std::vector<T>& data) {
clear(); if (row > 0 && col > 0) { m = row; n = col; unsigned int mxn = m * n; if (data.size() > 0) { for (unsigned int i = 0; i < mxn && i < data.size(); i++) { v.push_back(data[i]); } } }
}
template<class T> Matrix<T>& Matrix<T>::operator=(const Matrix<T>& other) {
clear(); if (other.m > 0 && other.n > 0) { m = other.m; n = other.n; unsigned int mxn = m * n; for (unsigned int i = 0; i < mxn && i < other.size(); i++) { v.push_back(other.v[i]); } } return *this;
}
template<class T> std::vector<T> Matrix<T>::operator[](unsigned int index) const {
std::vector<T> result; if (index >= m) { throw MatrixException(MATRIX_ERR_WRONG_ROW_INDEX); } else if ((index + 1) * n > size()) { throw MatrixException(MATRIX_ERR_TOO_FEW_DATA); } else { unsigned int begin = index * n; unsigned int end = begin + n; for (unsigned int i = begin; i < end; i++) { result.push_back(v[i]); } } return result;
}
template<class T> Matrix<T> Matrix<T>::operator*(const Matrix<T>& other) {
Matrix result(m, other.n); if (n != other.m) { throw MatrixException(MATRIX_ERR_MUL_ROW_AND_COL_NOT_EQUAL); } else if (m <= 0 || n <= 0 || other.n <= 0) { throw MatrixException(MATRIX_ERR_MUL_ROW_AND_COL_BE_GREATER_THAN_ZERO); } else if (m * n > size() || other.m * other.n > other.size()) { throw MatrixException(MATRIX_ERR_TOO_FEW_DATA); } else { for (unsigned int i = 0; i < m; i++) { for (unsigned int j = 0; j < other.n; j++) { T temp = v[i * n] * other.v[j]; for (unsigned int k = 1; k < n; k++) { temp += v[i * n + k] * other.v[k * other.n + j]; } result.v.push_back(temp); } } } return result;
}
- endif /* _MATRIX_H */
</lang>
Output:
22 28 49 64
Clojure
<lang lisp> (defn transpose
[s] (apply map vector s))
(defn nested-for
[f x y] (map (fn [a] (map (fn [b] (f a b)) y)) x))
(defn matrix-mult
[a b] (nested-for (fn [x y] (reduce + (map * x y))) a (transpose b)))
(def ma [[1 1 1 1] [2 4 8 16] [3 9 27 81] [4 16 64 256]]) (def mb [[4 -3 4/3 -1/4] [-13/3 19/4 -7/3 11/24] [3/2 -2 7/6 -1/4] [-1/6 1/4 -1/6 1/24]])</lang>
Output:
=> (matrix-mult ma mb) ((1 0 0 0) (0 1 0 0) (0 0 1 0) (0 0 0 1))
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)))</lang>
<lang lisp>(defun matrix-multiply (matrix1 matrix2)
(mapcar (lambda (row) (apply #'mapcar (lambda (&rest column) (apply #'+ (mapcar #'* row column))) matrix2)) matrix1))</lang>
The following version uses 2D arrays as inputs.
<lang lisp>(defun mmul (A B)
(let* ((m (car (array-dimensions A))) (n (cadr (array-dimensions A))) (l (cadr (array-dimensions B))) (C (make-array `(,m ,l) :initial-element 0))) (loop for i from 0 to (- m 1) do (loop for k from 0 to (- l 1) do (setf (aref C i k) (loop for j from 0 to (- n 1) sum (* (aref A i j) (aref B j k)))))) C))</lang>
Example use:
<lang lisp>(mmul #2a((1 2) (3 4)) #2a((-3 -8 3) (-2 1 4)))
- 2A((-7 -6 11) (-17 -20 25))
</lang>
Another version:
<lang lisp>(defun mmult (a b)
(loop with m = (array-dimension a 0) with n = (array-dimension a 1) with l = (array-dimension b 1) with c = (make-array (list m l) :initial-element 0) for i below m do (loop for k below l do (setf (aref c i k) (loop for j below n sum (* (aref a i j) (aref b j k))))) finally (return c)))</lang>
D
<lang d>import std.stdio, std.string, std.conv, std.numeric,
std.array, std.algorithm;
bool isRectangular(T)(in T[][] M) pure /*nothrow*/ {
//return all!(row => row.length == M[0].length)(M); return !canFind!(row => row.length != M[0].length)(M);
}
T[][] matrixMul(T)(in T[][] A, in T[][] B) pure nothrow in {
assert(isRectangular(A) && isRectangular(B) && !A.empty && !B.empty && A[0].length == B.length);
} body {
auto result = new T[][](A.length, B[0].length); auto aux = new T[B.length]; foreach (j; 0 .. B[0].length) { foreach (k, row; B) aux[k] = row[j]; foreach (i, ai; A) result[i][j] = dotProduct(ai, aux); } return result;
}
void main() {
immutable a = [[1, 2], [3, 4], [3, 6]]; immutable b = [[-3, -8, 3,], [-2, 1, 4]];
enum form = "[%([%(%d, %)],\n %)]]"; writefln("A = \n" ~ form ~ "\n", a); writefln("B = \n" ~ form ~ "\n", b); writefln("A * B = \n" ~ form, matrixMul(a, b));
}</lang>
- Output:
A = [[1, 2], [3, 4], [3, 6]] B = [[-3, -8, 3], [-2, 1, 4]] A * B = [[-7, -6, 11], [-17, -20, 25], [-21, -18, 33]]
Stronger Statically Typed Version
All array sizes are verified at compile-time (and no matrix is copied). Same output. <lang d>import std.stdio, std.string, std.conv, std.numeric,
std.array, std.algorithm, std.traits;
template TMMul(M1, M2) { // helper
alias Unqual!(typeof(M1[0][0]))[M2[0].length][M1.length] TMMul;
}
void matrixMul(T, T2, size_t k, size_t m, size_t n)
(in ref T[m][k] A, in ref T[n][m] B, /*out*/ ref T2[n][k] result) pure nothrow
if (is(T2 == Unqual!T)) {
T2[m] aux; foreach (j; 0 .. n) { foreach (k, row; B) aux[k] = row[j]; foreach (i, ai; A) result[i][j] = dotProduct(ai, aux); }
}
void main() {
immutable int[2][3] a = [[1, 2], [3, 4], [3, 6]]; immutable int[3][2] b = [[-3, -8, 3,], [-2, 1, 4]];
enum form = "[%([%(%d, %)],\n %)]]"; writefln("A = \n" ~ form ~ "\n", a); writefln("B = \n" ~ form ~ "\n", b); TMMul!(typeof(a), typeof(b)) result = void; matrixMul(a, b, result); writefln("A * B = \n" ~ form, result);
}</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>
Euphoria
<lang euphoria>function matrix_mul(sequence a, sequence b)
sequence c if length(a[1]) != length(b) then return 0 else c = repeat(repeat(0,length(b[1])),length(a)) for i = 1 to length(a) do for j = 1 to length(b[1]) do for k = 1 to length(a[1]) do c[i][j] += a[i][k]*b[k][j] end for end for end for return c end if
end function</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 } }
Fantom
Using a list of lists representation. The multiplication is done using three nested loops.
<lang fantom> class Main {
// multiply two matrices (with no error checking) public static Int[][] multiply (Int[][] m1, Int[][] m2) { Int[][] result := [,] m1.each |Int[] row1| { Int[] row := [,] m2[0].size.times |Int colNumber| { Int value := 0 m2.each |Int[] row2, Int index| { value += row1[index] * row2[colNumber] } row.add (value) } result.add (row) } return result }
public static Void main () { m1 := [[1,2,3],[4,5,6]] m2 := [[1,2],[3,4],[5,6]]
echo ("${m1} times ${m2} = ${multiply(m1,m2)}") }
} </lang>
Output:
[[1, 2, 3], [4, 5, 6]] times [[1, 2], [3, 4], [5, 6]] = [[22, 28], [49, 64]]
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 MATMUL intrinsic function to perform Matrix Multiply; use RESHAPE and SIZE intrinsic functions to form the matrices themselves: <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>
Frink
<lang gap>matprod[a is array, b is array] := {
c = makeArray[[length[a], length[b@0]], 0]
a_row = length[a]-1 a_col = length[a@0]-1 b_col = length[b]-1
for row = 0 to a_row { for col = 0 to b_col { for inc = 0 to a_col { c@row@col = c@row@col + (a@row@inc * b@inc@col) } } } return c
}</lang>
GAP
<lang gap># Built-in A := [[1, 2], [3, 4], [5, 6], [7, 8]]; B := [[1, 2, 3], [4, 5, 6]];
PrintArray(A);
- [ [ 1, 2 ],
- [ 3, 4 ],
- [ 5, 6 ],
- [ 7, 8 ] ]
PrintArray(B);
- [ [ 1, 2, 3 ],
- [ 4, 5, 6 ] ]
PrintArray(A * B);
- [ [ 9, 12, 15 ],
- [ 19, 26, 33 ],
- [ 29, 40, 51 ],
- [ 39, 54, 69 ] ]</lang>
Go
2D representation
<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 ]]
Flat representation
<lang go>package main
import "fmt"
type matrix struct {
ele []float64 stride int
}
func matrixFromRows(rows [][]float64) *matrix {
if len(rows) == 0 { return &matrix{nil, 0} } m := &matrix{make([]float64, len(rows)*len(rows[0])), len(rows[0])} for rx, row := range rows { copy(m.ele[rx*m.stride:(rx+1)*m.stride], row) } return m
}
func (m *matrix) print(heading string) {
if heading > "" { fmt.Print("\n", heading, "\n") } for e := 0; e < len(m.ele); e += m.stride { fmt.Printf("%6.3f ", m.ele[e:e+m.stride]) fmt.Println() }
}
func (m1 *matrix) multiply(m2 *matrix) (m3 *matrix, ok bool) {
if m1.stride*m2.stride != len(m2.ele) { return nil, false } m3 = &matrix{make([]float64, (len(m1.ele)/m1.stride)*m2.stride), m2.stride} for m1c0, m3x := 0, 0; m1c0 < len(m1.ele); m1c0 += m1.stride { for m2r0 := 0; m2r0 < m2.stride; m2r0++ { for m1x, m2x := m1c0, m2r0; m2x < len(m2.ele); m2x += m2.stride { m3.ele[m3x] += m1.ele[m1x] * m2.ele[m2x] m1x++ } m3x++ } } return m3, true
}
func main() {
a := matrixFromRows([][]float64{ {1, 1, 1, 1}, {2, 4, 8, 16}, {3, 9, 27, 81}, {4, 16, 64, 256}, }) b := matrixFromRows([][]float64{ { 4, -3, 4. / 3, -1. / 4, }, { -13. / 3, 19. / 4, -7. / 3, 11. / 24, }, { 3. / 2, -2, 7. / 6, -1. / 4, }, { -1. / 6, 1. / 4, -1. / 6, 1. / 24, }, }) p, ok := a.multiply(b) a.print("Matrix A:") b.print("Matrix B:") if !ok { fmt.Println("not conformable for matrix multiplication") return } p.print("Product of A and B:")
}</lang> Output is similar to 2D version.
Groovy
Without Indexed Loops
Uses transposition to avoid indirect element access via ranges of indexes. "assertConformability()" asserts that a & b are both rectangular lists of lists, and that row-length (number of columns) of a is equal to the column-length (number of rows) of b. <lang groovy>def assertConformable = { a, b ->
assert a instanceof List assert b instanceof List assert a.every { it instanceof List && it.size() == b.size() } assert b.every { it instanceof List && it.size() == b[0].size() }
}
def matmulWOIL = { a, b ->
assertConformable(a, b) def bt = b.transpose() a.collect { ai -> bt.collect { btj -> [ai, btj].transpose().collect { it[0] * it[1] }.sum() } }
}</lang>
Without Transposition
Uses ranges of indexes, the way that matrix multiplication is typically defined. Not as elegant, but it avoids expensive transpositions. Reuses "assertConformable()" from above. <lang groovy>def matmulWOT = { a, b ->
assertConformable(a, b) (0..<a.size()).collect { i -> (0..<b[0].size()).collect { j -> (0..<b.size()).collect { k -> a[i][k] * b[k][j] }.sum() } }
}</lang>
Test: <lang groovy>def m4by2 = [ [ 1, 2 ],
[ 3, 4 ], [ 5, 6 ], [ 7, 8 ] ]
def m2by3 = [ [ 1, 2, 3 ],
[ 4, 5, 6 ] ]
matmulWOIL(m4by2, m2by3).each { println it } println() matmulWOT(m4by2, m2by3).each { println it }</lang>
Output:
[9, 12, 15] [19, 26, 33] [29, 40, 51] [39, 54, 69] [9, 12, 15] [19, 26, 33] [29, 40, 51] [39, 54, 69]
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>
Icon and Unicon
Using the provided matrix library:
<lang icon> link matrix
procedure main ()
m1 := [[1,2,3], [4,5,6]] m2 := [[1,2],[3,4],[5,6]] m3 := mult_matrix (m1, m2) write ("Multiply:") write_matrix ("", m1) # first argument is filename, or "" for stdout write ("by:") write_matrix ("", m2) write ("Result: ") write_matrix ("", m3)
end </lang>
And a hand-crafted multiply procedure:
<lang icon> procedure multiply_matrix (m1, m2)
result := [] # to hold the final matrix every row1 := !m1 do { # loop through each row in the first matrix row := [] every colIndex := 1 to *m1 do { # and each column index of the result value := 0 every rowIndex := 1 to *m2 do { value +:= row1[rowIndex] * m2[rowIndex][colIndex] } put (row, value) } put (result, row) # add each row as it is complete } return result
end </lang>
Output:
Multiply: 1 2 3 4 5 6 by: 1 2 3 4 5 6 Result: 22 28 49 64
IDL
<lang idl>result = arr1 # arr2</lang>
J
Matrix multiply in J is +/ .*
. 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 a value in row x equals the corresponding value in 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
K
<lang k> (1 2;3 4)_mul (5 6;7 8) (19 22
43 50)</lang>
Liberty BASIC
There is no native matrix capability. A set of functions is available at http://www.diga.me.uk/RCMatrixFuncs.bas implementing matrices of arbitrary dimension in a string format. <lang lb> MatrixA$ ="4, 4, 1, 1, 1, 1, 2, 4, 8, 16, 3, 9, 27, 81, 4, 16, 64, 256" MatrixB$ ="4, 4, 4, -3, 4/3, -1/4 , -13/3, 19/4, -7/3, 11/24, 3/2, -2, 7/6, -1/4, -1/6, 1/4, -1/6, 1/24"
print "Product of two matrices"
call DisplayMatrix MatrixA$
print " *"
call DisplayMatrix MatrixB$
print " ="
MatrixP$ =MatrixMultiply$( MatrixA$, MatrixB$)
call DisplayMatrix MatrixP$
</lang>
Product of two matrices
| 1.00000 1.00000 1.00000 1.00000 |
| 2.00000 4.00000 8.00000 16.00000 |
| 3.00000 9.00000 27.00000 81.00000 |
| 4.00000 16.00000 64.00000 256.00000 |
| 4.00000 -3.00000 1.33333 -0.25000 |
| -4.33333 4.75000 -2.33333 0.45833 |
| 1.50000 -2.00000 1.16667 -0.25000 |
| -0.16667 0.25000 -0.16667 0.04167 |
=
| 1.00000 0.00000 0.00000 0.00000 |
| 0.00000 1.00000 0.00000 0.00000 |
| 0.00000 0.00000 1.00000 0.00000 |
| 0.00000 0.00000 0.00000 1.00000 |
Logo
<lang logo>TO LISTVMD :A :F :C :NV
- PROCEDURE LISTVMD
- A = LIST
- F = ROWS
- C = COLS
- NV = NAME OF MATRIX / VECTOR NEW
- this procedure transform a list in matrix / vector square or rect
(LOCAL "CF "CC "NV "T "W) MAKE "CF 1 MAKE "CC 1 MAKE "NV (MDARRAY (LIST :F :C) 1) MAKE "T :F * :C FOR [Z 1 :T][MAKE "W ITEM :Z :A MDSETITEM (LIST :CF :CC) :NV :W MAKE "CC :CC + 1 IF :CC = :C + 1 [MAKE "CF :CF + 1 MAKE "CC 1]] OUTPUT :NV END
TO XX
- MAIN PROGRAM
- LRCVS 10.04.12
- THIS PROGRAM multiplies two "square" matrices / vector ONLY!!!
- THE RECTANGULAR NOT WORK!!!
CT CS HT
- FIRST DATA MATRIX / VECTOR
MAKE "A [1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49] MAKE "FA 5 ;"ROWS MAKE "CA 5 ;"COLS
- SECOND DATA MATRIX / VECTOR
MAKE "B [2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50] MAKE "FB 5 ;"ROWS MAKE "CB 5 ;"COLS
IF (OR :FA <> :CA :FB <>:CB) [PRINT "Las_matrices/vector_no_son_cuadradas THROW
"TOPLEVEL ]
IFELSE (OR :CA <> :FB :FA <> :CB) [PRINT
"Las_matrices/vector_no_son_compatibles THROW "TOPLEVEL ][MAKE "MA LISTVMD :A
- FA :CA "MA MAKE "MB LISTVMD :B :FB :CB "MB] ;APPLICATION <<< "LISTVMD"
PRINT (LIST "THIS_IS: "ROWS "X "COLS) PRINT [] PRINT (LIST :MA "=_M1 :FA "ROWS "X :CA "COLS) PRINT [] PRINT (LIST :MB "=_M2 :FA "ROWS "X :CA "COLS) PRINT []
MAKE "T :FA * :CB
MAKE "RE (ARRAY :T 1)
MAKE "CO 0
FOR [AF 1 :CA][
FOR [AC 1 :CA][
MAKE "TEMP 0
FOR [I 1 :CA ][
MAKE "TEMP :TEMP + (MDITEM (LIST :I :AF) :MA) * (MDITEM (LIST :AC :I) :MB)]
MAKE "CO :CO + 1
SETITEM :CO :RE :TEMP]]
PRINT []
PRINT (LIST "THIS_IS: :FA "ROWS "X :CB "COLS)
SHOW LISTVMD :RE :FA :CB "TO ;APPLICATION <<< "LISTVMD"
END
- \
M1 * M2 RESULT / SOLUTION 1 3 5 7 9 2 4 6 8 10 830 1880 2930 3980 5030
11 13 15 17 19 12 14 16 18 20 890 2040 3190 4340 5490 21 23 25 27 29 X 22 24 26 28 30 = 950 2200 3450 4700 5950 31 33 35 37 39 32 34 36 38 40 1010 2360 3710 5060 6410 41 43 45 47 49 42 44 46 48 50 1070 2520 3970 5420 6870
- \
NOW IN LOGO!!!!
THIS_IS: ROWS X COLS
{{1 3 5 7 9} {11 13 15 17 19} {21 23 25 27 29} {31 33 35 37 39} {41 43 45 47 49}} =_M1 5 ROWS X 5 COLS
{{2 4 6 8 10} {12 14 16 18 20} {22 24 26 28 30} {32 34 36 38 40} {42 44 46 48 50}} =_M2 5 ROWS X 5 COLS
THIS_IS: 5 ROWS X 5 COLS
{{830 1880 2930 3980 5030} {890 2040 3190 4340 5490} {950 2200 3450 4700 5950}
{1010 2360 3710 5060 6410} {1070 2520 3970 5420 6870}}</lang>
Lua
<lang lua>function MatMul( m1, m2 )
if #m1[1] ~= #m2 then -- inner matrix-dimensions must agree return nil end
local res = {} for i = 1, #m1 do res[i] = {} for j = 1, #m2[1] do res[i][j] = 0 for k = 1, #m2 do res[i][j] = res[i][j] + m1[i][k] * m2[k][j] end end end return res
end
-- Test for MatMul mat1 = { { 1, 2, 3 }, { 4, 5, 6 } } mat2 = { { 1, 2 }, { 3, 4 }, { 5, 6 } } erg = MatMul( mat1, mat2 ) for i = 1, #erg do
for j = 1, #erg[1] do io.write( erg[i][j] ) io.write(" ") end io.write("\n")
end </lang>
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 / Octave
Matlab contains two methods of multiplying matrices: by using the "mtimes(matrix,matrix)" function, or the "*" operator.
<lang MATLAB>>> A = [1 2;3 4]
A =
1 2 3 4
>> B = [5 6;7 8]
B =
5 6 7 8
>> A * B
ans =
19 22 43 50
>> mtimes(A,B)
ans =
19 22 43 50</lang>
Maxima
<lang maxima>a: matrix([1, 2],
[3, 4], [5, 6], [7, 8])$
b: matrix([1, 2, 3],
[4, 5, 6])$
a . b; /* matrix([ 9, 12, 15],
[19, 26, 33], [29, 40, 51], [39, 54, 69]) */</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>
OxygenBasic
When using matrices in Video graphics, speed is important. Here is a matrix multiplier written in OxygenBasics's x86 Assembly code. <lang oxygenbasic>
'Example of matrix layout mapped to an array of 4x4 cells ' ' 0 4 8 C ' 1 5 9 D ' 2 6 A E ' 3 7 B F '
% MatrixType double
sub MatrixMul(MatrixType *A,*B,*C, sys n) '======================================== ' ' #if leftmatch matrixtype single % OneStep 4 % mtype single #endif ' #if leftmatch matrixtype double % OneStep 8 % mtype double #endif
sys pa=@A, pb=@B, pc=@C sys ColStep=OneStep*n
mov ecx,pa mov edx,pb mov eax,pc
mov esi,n ( call column : dec esi : jg repeat ) exit sub
column: '======
mov edi,n ( call cell : dec edi : jg repeat ) add edx,ColStep sub ecx,ColStep ret
cell: ' row A * column B '=======================
'matrix data is stored ascending vertically then horizontally 'thus rows are minor, columns are major ' push ecx push edx push eax mov eax,4 fldz ( fld mtype [ecx] fmul mtype [edx] faddp st1 add ecx,ColStep 'next column of matrix A add edx,OneStep 'next row of matrix B dec eax jnz repeat ) pop eax fstp mtype [eax] 'assign to next row of matrix C ' pop edx pop ecx add eax,OneStep 'next cell in column of matrix C (columns then rows) add ecx,OneStep 'next row of matrix A ret ' end sub
function ShowMatrix(MatrixType*A,sys n) as string '================================================ string cr=chr(13)+chr(10), tab=chr(9) function="MATRIX " n "x" n cr cr sys i,j,m ' for i=1 to n m=0 for j=1 to n function+=str( A[m+i] ) tab m+=n next function+=cr next end function
'TEST '====
% n 4 MatrixType A[n*n],B[n*n],C[n*n]
'reading vertically (minor) then left to right (major)
A <= 4,0,0,1, 0,4,0,0, 0,0,4,0, 0,0,0,4
B <= 2,0,0,2, 0,2,0,0, 0,0,2,0, 0,0,0,2
MatrixMul A,B,C,n
Print ShowMatrix C,n </lang>
PARI/GP
<lang parigp>M*N</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.
Perl 6
There are three ways in which this example differs significantly from the original Perl 5 code. These are not esoteric differences; all three of these features typically find heavy use in Perl 6.
First, we can use a real signature that can bind two arrays as arguments, because the default in Perl 6 is not to flatten arguments unless the signature specifically requests it. We don't need to pass the arrays with backslashes because the binding choice is made lazily by the signature itself at run time; in Perl 5 this choice must be made at compile time. Also, we can bind the arrays to formal parameters that are really lexical variable names; in Perl 5 they can only be bound to global array objects (via a typeglob assignment).
Second, we use the X cross operator in conjunction with a two-parameter closure to avoid writing nested loops. The X cross operator, along with Z, the zip operator, is a member of a class of operators that expect lists on both sides, so we call them "list infix" operators. We tend to define these operators using capital letters so that they stand out visually from the lists on both sides. The cross operator makes every possible combination of the one value from the first list followed by one value from the second. The right side varies most rapidly, just like an inner loop. (The X and Z operators may both also be used as meta-operators, Xop or Zop, distributing some other operator "op" over their generated list. All metaoperators in Perl 6 may be applied to user-defined operators as well.)
Third is the use of prefix ^ to generate a list of numbers in a range. Here it is used on an array to generate all the indexes of the array. We have a way of indicating a range by the infix .. operator, and you can put a ^ on either end to exclude that endpoint. We found ourselves writing 0 ..^ @a so often that we made ^@a a shorthand for that. It's pronounced "upto". The array is evaluated in a numeric context, so it returns the number of elements it contains, which is exactly what you want for the exclusive limit of the range.
<lang perl6>sub mmult(@a,@b) {
my @p; for ^@a X ^@b[0] -> $r, $c { @p[$r][$c] += @a[$r][$_] * @b[$_][$c] for ^@b; } @p;
}
my @a = [1, 1, 1, 1],
[2, 4, 8, 16], [3, 9, 27, 81], [4, 16, 64, 256];
my @b = [ 4 , -3 , 4/3, -1/4 ],
[-13/3, 19/4, -7/3, 11/24], [ 3/2, -2 , 7/6, -1/4 ], [ -1/6, 1/4, -1/6, 1/24];
.say for mmult(@a,@b);</lang>
Output:
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
Note that these are not rounded values, but exact, since all the math was done in rationals. Hence we need not rely on format tricks to hide floating-point inaccuracies.
Just for the fun of it, here's a functional version that uses no temp variables or side effects. Some people will find this more readable and elegant, and others will, well, not.
<lang perl6>sub mmult(@a,@b) {
for ^@a -> $r {[ for ^@b[0] -> $c { [+] (@a[$r][^@b] Z* @b[^@b]»[$c]) } ]}
}</lang>
Here we use Z with an "op" of *, which is a zip with multiply. This, along with the [+] reduction operator, replaces the inner loop. We chose to split the outer X loop back into two loops to make it convenient to collect each subarray value in [...]. It just collects all the returned values from the inner loop and makes an array of them. The outer loop simply returns the list of arrays. We also had to sneak in a » metaoperator (known as "hyper") to do a parallel subscript lookup. Eventually we'll have shaped arrays, and a multidimensional subscript will automatically slice in all its dimensions--but rakudo doesn't do that yet.
PicoLisp
<lang PicoLisp>(de matMul (Mat1 Mat2)
(mapcar '((Row) (apply mapcar Mat2 '(@ (sum * 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>
Another one, use numpy the most popular array package for python
<lang python>
import numpy as np
np.dot(a,b)
- or if a is an array
a.dot(b)</lang>
R
<lang r>a %*% b</lang>
Racket
<lang scheme>(define (m-mult m1 m2)
(for/list ([r m1]) (for/list ([c (apply map list m2)]) (apply + (map * r c)))))</lang>
Rascal
<lang Rascal>public rel[real, real, real] matrixMultiplication(rel[real x, real y, real v] matrix1, rel[real x, real y, real v] matrix2){ if (max(matrix1.x) == max(matrix2.y)){ p = {<x1,y1,x2,y2, v1*v2> | <x1,y1,v1> <- matrix1, <x2,y2,v2> <- matrix2};
result = {}; for (y <- matrix1.y){ for (x <- matrix2.x){ v = (0.0 | it + v | <x1, y1, x2, y2, v> <- p, x==x2 && y==y1, x1==y2 && y2==x1); result += <x,y,v>; } } return result; } else throw "Matrix sizes do not match.";
//a matrix, given by a relation of the x-coordinate, y-coordinate and value. public rel[real x, real y, real v] matrixA = { <0.0,0.0,12.0>, <0.0,1.0, 6.0>, <0.0,2.0,-4.0>, <1.0,0.0,-51.0>, <1.0,1.0,167.0>, <1.0,2.0,24.0>, <2.0,0.0,4.0>, <2.0,1.0,-68.0>, <2.0,2.0,-41.0> };</lang>
REXX
<lang rexx>/*REXX program multiplies two matrixes together, shows matrixes & result*/ x.= x.1='1 2' x.2='3 4' x.3='5 6' x.4='7 8'
do r=1 while x.r\== /*build the "A" matric from X. #s*/ do c=1 while x.r\==; parse var x.r a.r.c x.r; end end
Arows=r-1; Acols=c-1 y.= y.1='1 2 3' y.2='4 5 6'
do r=1 while y.r\== /*build the "B" matric from Y. #s*/ do c=1 while y.r\==; parse var y.r b.r.c y.r; end end
Brows=r-1; Bcols=c-1 c.=0; L=0 /*L is max width of an element. */
do i =1 for Arows /*multiply matrix A & B ──► C */ do j =1 for Bcols do k=1 for Acols c.i.j = c.i.j + a.i.k * b.k.j; L=max(L,length(c.i.j)) end /*k*/ end /*j*/ end /*i*/
call showMat 'A',Arows,Acols call showMat 'B',Brows,Bcols call showMat 'C',Arows,Bcols exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────SHOWMAT subroutine──────────────────*/ showMat: parse arg mat,rows,cols; say say center(mat 'matrix',cols*(L+1)+4,"─")
do r =1 for rows; _= do c=1 for cols; _=_ right(value(mat'.'r'.'c),L); end; say _ end
return</lang> output
─A matrix─ 1 2 3 4 5 6 7 8 ──B matrix─── 1 2 3 4 5 6 ──C matrix─── 9 12 15 19 26 33 29 40 51 39 54 69
Ruby
Using 'matrix' from the standard library: <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>
UNIX Shell
<lang bash>
- !/bin/bash
DELAY=0 # increase this if printing of matrices should be slower
echo "This script takes two matrices, henceforth called A and B, and returns their product, AB.
For the time being, matrices can have integer components only.
"
read -p "Number of rows of matrix A: " arows read -p "Number of columns of matrix A: " acols brows="$acols" echo echo "Number of rows of matrix B: "$brows read -p "Number of columns of matrix B: " bcols
crows="$arows" ccols="$bcols" echo
echo "Number of rows of matrix AB: " $crows echo "Number of columns of matrix AB: " $ccols echo echo
matrixa=( ) matrixb=( )
- input matrix A
maxlengtha=0 for ((row=1; row<=arows; row++)); do
for ((col=1; col<=acols; col++)); do
checkentry="false" while [ "$checkentry" != "true" ]; do read -p "Enter component A[$row, $col]: " number index=$(((row-1)*acols+col)) matrixa[$index]="$number" [ "${matrixa[$index]}" -eq "$number" ] && checkentry="true" echo done entry="${matrixa[$index]}" [ "${#entry}" -gt "$maxlengtha" ] && maxlengtha="${#entry}"
done echo
done
- print matrix A to guard against errors
if [ "$maxlengtha" -le "5" ]; then
width=8
else
width=$((maxlengtha + 3))
fi
echo "This is matrix A:
"
for ((row=1; row<=arows; row++)); do
for ((col=1; col<=acols; col++)); do
index=$(((row-1)*acols+col)) printf "%${width}d" "${matrixa[$index]}" sleep "$DELAY"
done echo; echo # printf %s "\n\n" does not work...
done
echo echo
- input matrix B
maxlengthb=0 for ((row=1; row<=brows; row++)); do
for ((col=1; col<=bcols; col++)); do
checkentry="false" while [ "$checkentry" != "true" ]; do read -p "Enter component B[$row, $col]: " number index=$(((row-1)*bcols+col)) matrixb[$index]="$number" [ "${matrixb[$index]}" -eq "$number" ] && checkentry="true" echo done entry="${matrixb[$index]}" [ "${#entry}" -gt "$maxlengthb" ] && maxlengthb="${#entry}"
done echo
done
- print matrix B to guard against errors
if [ "$maxlengthb" -le "5" ]; then
width=8
else
width=$((maxlengthb + 3))
fi
echo "This is matrix B:
"
for ((row=1; row<=brows; row++)); do
for ((col=1; col<=bcols; col++)); do
index=$(((row-1)*bcols+col)) printf "%${width}d" "${matrixb[$index]}" sleep "$DELAY"
done echo; echo # printf %s "\n\n" does not work...
done
read -p "Hit enter to continue"
- calculate matrix C := AB
maxlengthc=0 time for ((row=1; row<=crows; row++)); do
for ((col=1; col<=ccols; col++)); do
# calculate component C[$row, $col]
runningtotal=0 for ((j=1; j<=acols; j++)); do rowa="$row" cola="$j" indexa=$(((rowa-1)*acols+cola)) rowb="$j" colb="$col" indexb=$(((rowb-1)*bcols+colb))
entry_from_A=${matrixa[$indexa]} entry_from_B=${matrixb[$indexb]}
subtotal=$((entry_from_A * entry_from_B)) ((runningtotal+=subtotal)) done
number="$runningtotal"
# store component in the result array index=$(((row-1)*ccols+col)) matrixc[$index]="$number"
entry="${matrixc[$index]}" [ "${#entry}" -gt "$maxlengthc" ] && maxlengthc="${#entry}"
done
done
echo read -p "Hit enter to continue" echo
- print the matrix C
if [ "$maxlengthc" -le "5" ]; then
width=8
else
width=$((maxlengthc + 3))
fi
echo "The product matrix is:
"
for ((row=1; row<=crows; row++)); do
for ((col=1; col<=ccols; col++)); do
index=$(((row-1)*ccols+col)) printf "%${width}d" "${matrixc[$index]}" sleep "$DELAY"
done echo; echo # printf %s "\n\n" does not work...
done
echo echo </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>
- Programming Tasks
- Matrices
- Ada
- ALGOL 68
- APL
- AutoHotkey
- BASIC
- BBC BASIC
- C
- C sharp
- C++
- Blitz++
- Clojure
- Common Lisp
- D
- ELLA
- Euphoria
- Factor
- Fantom
- Forth
- Forth Scientific Library
- Fortran
- Frink
- GAP
- Go
- Groovy
- Haskell
- HicEst
- Icon
- Unicon
- IDL
- J
- Java
- JavaScript
- K
- Liberty BASIC
- Logo
- Lua
- Mathematica
- MATLAB
- Octave
- Maxima
- Nial
- OCaml
- OxygenBasic
- PARI/GP
- Perl
- Perl 6
- PicoLisp
- PL/I
- Pop11
- Prolog
- PureBasic
- Python
- R
- Racket
- Rascal
- REXX
- Ruby
- Scala
- Scheme
- Seed7
- SQL
- Tcl
- TI-83 BASIC
- TI-89 BASIC
- UNIX Shell
- Ursala
- ZPL