Perform simple mathematical operations on two matrices

From CodeCodex

Related content:

Implementations[edit]

Algol 68[edit]

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))

Java[edit]

40px

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));
 
   }
}

Perl[edit]

TIMTOWTDI. See Math::Cephes::Matrix or Math::MatrixReal or PDL::MatrixOps or Math::Matrix.

Ruby[edit]

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