Source code in the xpl0 programming language

proc Mat4x1Mul(M, V);   \Multiply matrix M times column vector V
real M,     \4x4 matrix  [M] * [V] -> [V]
     V;     \column vector
real W(4);  \working copy of column vector
int  R;     \row
[for R:= 0 to 4-1 do
    W(R):= M(R,0)*V(0) + M(R,1)*V(1) + M(R,2)*V(2) + M(R,3)*V(3);
for R:= 0 to 4-1 do V(R):= W(R);
];

proc Mat4x4Mul(M, N);   \Multiply matrix M times matrix N
real M, N;   \4x4 matrices       [M] * [N] -> [N]
real W(4,4); \working copy of matrix N
int  C;      \column
[for C:= 0 to 4-1 do
       [W(0,C):= M(0,0)*N(0,C) + M(0,1)*N(1,C) + M(0,2)*N(2,C) + M(0,3)*N(3,C);
        W(1,C):= M(1,0)*N(0,C) + M(1,1)*N(1,C) + M(1,2)*N(2,C) + M(1,3)*N(3,C);
        W(2,C):= M(2,0)*N(0,C) + M(2,1)*N(1,C) + M(2,2)*N(2,C) + M(2,3)*N(3,C);
        W(3,C):= M(3,0)*N(0,C) + M(3,1)*N(1,C) + M(3,2)*N(2,C) + M(3,3)*N(3,C);
        ];
for C:= 0 to 4-1 do
       [N(0,C):= W(0,C);
        N(1,C):= W(1,C);
        N(2,C):= W(2,C);
        N(3,C):= W(3,C);
        ];
];