Here are a few routines (and a header file) culminating in a routine to produce random orthogonal matrices.

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