The following multiplies two dense matrices.

function matrix_multiply(A,B) =
  {{sum({x*y: x in rowA; y in columnB})
    : columnB in transpose(B)}
   : rowA in A} $

-----------------

 Parallel step: x*y for (x in rowA; y in columnB)

-----------------


The following inverts a dense matrix. 
It uses Gauss-Jordan elimination without pivoting.

function Gauss_Jordan(A,i) =
if (i == #A) 
then 
  A
else 
    let
	(irow,Ap) = head_rest(A);
	val = irow[i];
	irow = {v/val : v in irow};
	Ap = {let scale = jrow[i]
	      in {v - scale*x : x in irow; v in jrow}
	      : jrow in Ap}
in Gauss_Jordan(Ap++[irow],i+1) $

function matrix_inverse(A) =
  let
    n = #A;

    % Pad the matrix with the identity matrix (i.e. A ++ I) %
    Ap = {row ++ rep(dist(0.,n),1.0,i):
	 row in A; i in [0:n]};

    % Run Gauss-Jordan elimination on padded matrix %
    Ap = Gauss_Jordan(Ap,0);

  % Drop the identity matrix at the front %
  in {drop(row,n) : row in Ap} $


-----------------

 The following solves a linear system of the form Ax = b. 
It uses Gauss-Jordan elimination without pivoting 
(the same routine as used in Matrix Inversion above).

function Gauss_Jordan(A,i) =
if (i == #A) then A
else 
    let
	(irow,Ap) = head_rest(A);
	val = irow[i];
	irow = {v/val : v in irow};
	Ap = {let scale = jrow[i]
	      in {v - scale*x : x in irow; v in jrow}
	      : jrow in Ap}
in Gauss_Jordan(Ap++[irow],i+1) $


function linear_solve(A,b) =
  let
    n = #A;

    % Pad the matrix with b %
    Ap = {row ++ [v]: row in A; v in b};

    % Run Gauss-Jordan elimination on padded matrix %
    Ap = Gauss_Jordan(Ap,0);

  % Extract last colum %
  in {row[n] : row in Ap} $