/* Matrix multiplication of A by B, yielding C */
MMULT: procedure (a, b, c);
   declare (a, b, c)(*,*) float controlled;
   declare (i, j, m, n, p) fixed binary;

   if hbound(a,2) ^= hbound(b,1) then
      do;
         put skip list
            ('Matrices are incompatible for matrix multiplication');
         signal error;
      end;

   m = hbound(a, 1); p = hbound(b, 2);
   if allocation(c) > 0 then free c;

   allocate c(m,p);

   do  i = 1 to m;
      do j = 1 to p;
         c(i,j) = sum(a(i,*) * b(*,j) );
      end;
   end;
end MMULT;