Perform simple mathematical operations on two matrices
From CodeCodex
| Image:Lightbulb.png Related content: |
| <DPL>
category=Math shownamespace=false </DPL> |
Contents |
[edit] Implementations
[edit] Algol 68
Basic matrix operations.
MODE MATELT = LONG REAL; # any suitable real type #
# matrices should be [,] MATELT. lower bounds must be 1. #
PROC ident = (INT rank) [,] MATELT :
# creates an identity square matrix of the specified rank. #
BEGIN
[1 : rank, 1 : rank] MATELT m;
FOR i TO rank DO
FOR j TO rank DO
m[i, j] := (i = j | 1 | 0)
OD
OD;
m
END;
OP + = ([,] MATELT a, b) [,] MATELT :
# sum of two matrices, which must have identical dimensions. #
BEGIN
INT
nr rows = 1 UPB a # = 1 UPB b #,
nr cols = 2 UPB a # = 2 UPB b #;
[1 : nr rows, 1 : nr cols] MATELT c;
FOR i TO nr rows DO
FOR j TO nr cols DO
c[i, j] := a[i, j] + b[i, j]
OD
OD;
c
END;
OP * = ([,] MATELT a, b) [,] MATELT :
# product of two matrices, which must conform -- second dimension
of a must equal first dimension of b. #
BEGIN
INT
nr rows = 1 UPB a,
nr cols = 2 UPB b,
nr conform = 1 UPB b # = 2 UPB a #;
[1 : nr rows, 1 : nr cols] MATELT c;
FOR i TO nr rows DO
FOR j TO nr cols DO
MATELT sum := 0;
FOR k TO nr conform DO
sum +:= a[i, k] * b[k, j]
OD;
c[i, j] := sum
OD
OD;
c
END;
OP * = ([,] MATELT a, MATELT b) [,] MATELT :
# scalar matrix product. #
BEGIN
INT nr rows = 1 UPB a, nr cols = 2 UPB a;
[1 : nr rows, 1 : nrcols] MATELT c;
FOR i TO nr rows DO
FOR j TO nr cols DO
c[i, j] := a[i, j] * b
OD
OD;
c
END;
OP * = (MATELT a, [,] MATELT b) [,] MATELT :
# scalar matrix product. #
b * a;
PROC detx = (PROC (INT, INT) MATELT m, INT rank) MATELT :
# internal routine for determinant and inverse computations.
m is callback to do actual accessing of matrix elements. #
CASE rank IN
m(1, 1),
m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)
OUT
MATELT result := 0;
INT sign := 1;
FOR i TO rank DO
MATELT term =
sign
*
m(i, 1)
*
detx
(
(INT i2, j2) MATELT :
m
(
(i2 < i | i2 | i2 + 1),
j2 + 1
),
rank - 1
);
result +:= term;
sign := - sign
OD;
result
ESAC;
PROC det = ([,] MATELT m) MATELT :
# computes determinant of matrix m, which must be square. #
# Method used is "expansion by minors"
<http://mathworld.wolfram.com/DeterminantExpansionbyMinors.html>
which is not the most efficient, but does illustrate how compactly
things can be expressed using callbacks. #
detx
(
(INT i, j) MATELT : m[i, j],
1 UPB m
);
PROC inv = ([,] MATELT m) [,] MATELT :
# computes inverse of matrix m, which must be square.
This is done using minors
<http://mathworld.wolfram.com/MatrixInverse.html>. #
BEGIN
INT rank = 1 UPB m # = 2 UPB m #;
[1 : rank, 1 : rank] MATELT inv m;
CASE rank IN
inv m[1, 1] := 1.0 / m[1, 1]
OUT
MATELT det m = det(m);
INT sign := 1;
FOR i TO rank DO
FOR j TO rank DO
inv m[j, i] :=
sign
*
detx
(
(INT i2, j2) MATELT :
m
[
(i2 < i | i2 | i2 + 1),
(j2 < j | j2 | j2 + 1)
],
rank - 1
)
/
det m;
sign := - sign
OD;
IF NOT ODD rank THEN
sign := - sign
FI
OD
ESAC;
inv m
END;
PROC make matrix =
(
INT nr rows,
INT nr cols,
PROC (INT, INT) MATELT val
) [,] MATELT :
# returns a matrix with elements computed by val. #
BEGIN
[1 : nr rows, 1 : nrcols] MATELT m;
FOR i TO nr rows DO
FOR j TO nr cols DO
m[i, j] := val(i, j)
OD
OD;
m
END;
PROC put matrix = (REF FILE out, [,] MATELT m) VOID :
# prints a formatted representation of matrix m to out. #
BEGIN
FOR i TO 1 UPB m DO
put(out, "[");
FOR j TO 2 UPB m DO
IF j /= 1 THEN
put(out, " ")
FI;
put(out, m[i, j])
OD;
put(out, ("]", new line))
OD
END;
CO some tests of above routines CO
[,] MATELT m = make matrix(4, 4, (INT i, j) MATELT : 1.0 / (2 * i + j));
# something not too far removed from a Hilbert matrix,
quite badly conditioned #
[,] MATELT inv m = inv(m);
print(("m:", new line));
put matrix(standout, m);
print(("m + 1:", new line));
put matrix(standout, m + ident(4));
print(("inv(m):", new line));
put matrix(standout, inv m);
print(("m * inv(m):", new line)); # should be close to identity matrix #
put matrix(standout, m * inv m);
print(("det(m) = ", det(m), new line))
[edit] Java
| This code requires JMathTools, an external library, to run. |
import static org.math.array.LinearAlgebra.*;
public class LinearAlgebraExample {
public static void main(String[] args) {
// random 4 x 3 matrix
double[][] A = random(4, 3);
// random 4 x 3 matrix
double[][] B = random(4, 3);
// random 4 x 4 matrix + Id
double[][] C = plus(random(4, 4),identity(4));
// linear algebra
double[][] D = plus(A, B);
double[][] E = plus(A, 1);
double[][] F = minus(A, B);
double[][] G = minus(2, A);
double[][] H = times(A, transpose(B));
double[][] I = times(A, 10);
double[][] J = divide(C, B);
double[][] K = divide(A, 10);
double[][] L = inverse(C);
// print matrices in command line
System.out.println("A = \n" + tostring(A));
System.out.println("B = \n" + tostring(B));
System.out.println("C = \n" + tostring(C));
System.out.println("D = A + B = \n" + tostring(D));
System.out.println("E = A + 1 = \n" + tostring(E));
System.out.println("F = A - B = \n" + tostring(F));
System.out.println("G = 2 - A = \n" + tostring(G));
System.out.println("H = A * t(B) = \n" + tostring(H));
System.out.println("I = A * 10 = \n" + tostring(I));
System.out.println("J = C / B = \n" + tostring(J));
System.out.println("K = A / 10 = \n" + tostring(K));
System.out.println("L = C^-1 = \n" + tostring(L));
}
}
[edit] Perl
TIMTOWTDI. See Math::Cephes::Matrix or Math::MatrixReal or PDL::MatrixOps or Math::Matrix.
[edit] Ruby
require 'matrix' puts a = Matrix[[11, 12, 13], [21, 22, 23], [31, 32, 33]] puts b = Matrix.identity(3) puts c = Matrix[[3.0, 2.0, 1.0], [1.0, 3.0, 1.0], [2.0, 2.0, 1.0]] puts d = a + b puts (e = a + 1) rescue p $! puts f = a - b puts (g = 2 - a) rescue p $! puts h = a * b.transpose puts i = a * 10 puts j = c / b puts k = a / 10 puts l = c.inverse
