Random Normal Deviate: rnd() returns a random normal deviate. This is a normal (Gaussian) distribution with zero mean and standard deviation of 1. The pdf (probability density function) is e−x2/2/√(2π).
Matrix Multiply: mm(m, n, j, A, B, P) places the product of the m by n matrix A, and the n by j matrix B, in the space P for an m by j matrix.
Gram-Schmidt: For 0 < m ≤ n, gs(m, n, A) transforms the m by n matrix A, into an orthogonal matrix. For any k ≤ m the first k rows of the resulting matrix span the same subspace as the first k rows of the original matrix. For 0 < i, j ≤ m the inner product of row i and row j is δij.
Random Orthogonal Matrix: For 0 < j ≤ k, grom(j, k, A) places a random orthogonal j by k matrix in the space A. The probability distribution is invariant under the Haar measure of the Lie group O(k).
A demo.
clang demo.c lob.c ./a.out