inv : float array array -> float array array Inverts the matrix or raises Linear.Singular for singular matrices.

mul : float array array -> float array array -> float array array multiplies two matrices.

determ : float array array -> float returns the determinant.

trans : 'a array array -> 'a array array returns the transpose or raises Linear.ragged if vectors are of varying length.

ip : float array -> float array -> float returns inner product (sum of termwise products).

mtv : float array array -> float array -> float array Matrix times vector.

vneg : float array -> float array negates a vector.

iden : int -> float array array creates n by n identity matrix.

sm : float -> float array -> float array scalar times vector yielding vector.

vadd : float array -> float array -> float array adds two vectors.

gs : (float array -> float array -> float) -> float array array -> float array array computes Gram-Schmidt of 2nd argument. For an n by m matrix C, for any j from 1 to n, the first j vectors of C, and the first j vectors Y = gs g C span the same space. g is a positive definite bilinear form and for 0≤i<n and 0≤i<n
g C_i Y_j = δ_ij
g encapsulates a metric tensor. The conventional Gram-Schmidt is gs ip.

eigen : float array array -> (float * float array) array computes eigenvalues and vectors of real symmetric matrices. The algorithm is Jacobi’s.

module Fac = struct
open Linear
let oth eig = aini (len eig)
  (fun j -> match eig.(j) with (evl, evc) -> evc)
let fac symm = let eig = eigen symm and l = len symm in
  (oth eig, aini l (fun i -> aini l (fun j ->
    if i = j then match eig.(j) with (a, b) -> a else 0.)))
end;;