Matrix multiplication

Matrix multiplication
You are encouraged to solve this task according to the task description, using any language you may know.

Multiply two Matrices together, they can be of any dimension as long as the number of columns of the first matrix is equal to the number of rows of the second matrix

Ada

package Matrix_Ops is
   type Matrix is array(Natural range <>, Natural range <>) of Float;
   function "*" (Left, Right : Matrix) return Matrix;
   Dimension_Violation : exception;
end Matrix_Ops;

package body Matrix_Ops is

   ---------
   -- "*" --
   ---------

   function "*" (Left, Right : Matrix) return Matrix is
      Temp : Matrix(Left'Range(1), Right'Range(2)) := (Others =>(Others => 0.0));
   begin
      if Left'Length(2) /= Right'Length(1) then
         raise Dimension_Violation;
      end if;
     
      for I in Left'range(1) loop
         for J in Right'range(2) loop
            for K in Left'range(2) loop
               Temp(I,J) := Temp(I,J) + Left(I, K)*Right(K, J);
            end loop;
         end loop;
      end loop;
      return Temp;
   end "*";

end Matrix_Ops;

ALGOL 68

An example of user defined Vector and Matrix Multiplication Operators:

MODE FIELD = LONG REAL; # field type is LONG REAL #

# crude exception handling #
PROC VOID raise index error := VOID: GOTO exception index error;

# define the vector/matrix operators #
OP * = (FLEX []FIELD 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 * = ([]FIELD a, [,]FIELD b)[]FIELD: ( # overload vec 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
  );

OP * = ([,]FIELD a, b)[,]FIELD: ( # 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 #
[,]FIELD a=((1,  1,  1,   1), # matrix A #
            (2,  4,  8,  16),
            (3,  9, 27,  81),
            (4, 16, 64, 256));

[,]FIELD 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));

[,]FIELD prod = a * b; # actual multiplication example of A x B #

FORMAT field fmt = $g(-6,2)$; # width of 6, with no '+' sign, 2 decimals #
FORMAT vec fmt = $"("n(2 UPB prod-1)(f(field fmt)",")f(field fmt)")"$;
FORMAT matrix fmt = $x"("n(UPB prod-1)(f(vec fmt)","lxx)f(vec fmt)");"$;
 
FORMAT result fmt = $"Product of a and b: "lf(matrix fmt)l$;

# finally print the result #
printf((result fmt,prod))
EXIT 

exception index error: 
  putf(stand error, $x"Error: Array bound 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));

Common Lisp

(defun matrix-multiply (a b)
  (flet ((col (mat i) (mapcar #'(lambda (row) (elt row i)) mat))
         (row (mat i) (elt mat i)))
    (loop for row from 0 below (length a)
          collect (loop for col from 0 below (length (row b 0))
                        collect (apply #'+ (mapcar #'* (row a row) (col b col)))))))

;; example use:
(matrix-multiply '((1 2) (3 4)) '((-3 -8 3) (-2 1 4)))

(defun matrix-multiply (matrix1 matrix2)
 (mapcar
  (lambda (row)
   (apply #'mapcar
    (lambda (&rest column)
     (apply #'+ (mapcar #'* row column))) matrix2)) matrix1))

Haskell

A somewhat inefficient version with lists (transpose is expensive):

 import Data.List

 mmult :: Num a => [[a]] -> [[a]] -> [[a]] 
 mmult a b = [ [ sum $ zipWith (*) ar bc | bc <- (transpose b) ] | ar <- a ]
 
 -- Example use:
 test = [[1, 2],
         [3, 4]] `mmult` [[-3, -8, 3],
                          [-2,  1, 4]]

A more efficient version, based on arrays:

 import Data.Array
 
 mmult :: (Ix i, Num a) => Array (i,i) a -> Array (i,i) a -> Array (i,i) a 
 mmult x y 
   | x1 /= y0 || x1' /= y0'  = error "range mismatch"
   | otherwise               = array ((x0,y1),(x0',y1')) l
   where
     ((x0,x1),(x0',x1')) = bounds x
     ((y0,y1),(y0',y1')) = bounds y
     ir = range (x0,x0')
     jr = range (y1,y1')
     kr = range (x1,x1')
     l  = [((i,j), sum [x!(i,k) * y!(k,j) | k <- kr]) | i <- ir, j <- jr]

J

x +/ .* y

where x and y are conformable arrays (trailing dimension of array x equals the leading dimension of array y). The notation is for a generalized inner product so that

x ~:/ .*. y   NB. boolean inner product (~: is "not equal" (exclusive or) and *. is "and")
x *./ .=  y   NB. which rows of x are the same as vector y?
x + / .=  y   NB. number of places where each row of x equals vector y

etc.

Java

public static double[][] mult(double a[][], double b[][]){//a[m][n], b[n][p]
   if(a[0].length != b.length) return null; //invalid dims

   int n = a[0].length;
   int m = a.length;
   int p = b[0].length;

   double ans[][] = new double[m][p];

   for(int i = 0;i < m;i++){
      for(int j = 0;j < p;j++){
         for(int k = 0;k < n;k++){
            ans[i][j] += a[i][k] * b[k][j];
         }
      }
   }
   return ans;
}

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;