Below are functions and example programs I have written using the GNU Scientific Library (GSL). There are more waiting in queue until I have time to put them here.

gsl-matrix-list-btac.c.
Synopsis: Print a matrix list to a binary file, then read it and print it to stdout in reverse order.
Input: [rows] [columns] [number of matrices] [filename]

gsl-matrix-gram-schmidt.c
Synopsis: Randomly fills a square matrix, orthonormalizes the columns, and checks the results.
Input: [rows] [columns]
About: Uses gsl_vector_normalize and gsl_matrix_gram-schmidt.

gsl-matrix-list-tac.c (tac is the opposite of cat)
Synopsis: Read a matrix list in from a file, reverse it, and output it to a file.
Input: [rows] [columns] [input filename] [output filename]
About: Uses gsl_matrix_fscanf to read and gsl_matrix_fprintf to write.

gsl-matrix-list-fprintf.c
Synopsis: Print a list of random matrices to a file.
Input: [rows] [columns] [number of matrices] [output filename]
for uniform distribution of random numbers in [0,1].

To compile and link

> gcc -Wall -I/path/to/include -c [name].c
> gcc -L/path/to/lib [name].o -o [name] -lgsl -lgslcblas -lm

Running ./[name] will give the proper usage.

Print a list of 10 random 2-by-3 matrices to in.dat

./gsl-matrix-list-fprintf 2 3 10 in.dat

./gsl-matrix-list-tac 2 3 in.dat out.dat

void gsl_matrix_tac(FILE *in, FILE *out, int rows, int columns)
Synopsis: This function recursively reads a matrix list from file and outputs its reverse. Not efficient for large lists. See gsl-matrix-list-btac.c below.
About: Uses gsl_matrix_fscan.

int gsl_matrix_fscan(FILE *in, gsl_matrix *A)
Synopsis: This function scans a matrix from FILE *in unless EOF is reached.

void gsl_matrix_gram_schmidt(gsl_vector *v, gsl_matrix *A, int flag)
Synopsis: This function Gram-Schmidt orthonormalizes a square matrix A. The flag = 0 applies the GS procedure to the rows of A and flag = 1 applies it to the columns of A. The norms of A are stored in the vector v.
About: Uses gsl_vector_normalize. One could alternatively use the QR decomposition with the GSL functions: gsl_linalg_QR_decomp, gsl_linalg_QR_solve, gsl_linalg_QR_unpack (see here).

void gsl_vector_normalize(gsl_vector *a)
Synopsis: Normalizes a vector.

GSL manual
srand48 Seed uniformly distributed double-precision pseudo-random number generator.
drand48 Generate uniformly distributed pseudo-random, non-negative, double-precision, floating-point values, uniformly distributed over the interval [0,1].