define matmul(a, b) -> c;
    lvars ba = boundslist(a), bb = boundslist(b);
    lvars i, i0 = ba(1), i1 = ba(2);
    lvars j, j0 = bb(1), j1 = bb(2);
    lvars k, k0 = bb(3), k1 = bb(4);
    if length(ba) /= 4 then
        throw([need_2d_array ^a])
    endif;
    if length(bb) /= 4 then
        throw([need_2d_array ^b])
    endif;
    if ba(3) /= j0 or ba(4) /= j1 then
        throw([dimensions_do_not_match ^a ^b]);
    endif;
    newarray([^i0 ^i1 ^k0 ^k1], 0) -> c;
    for i from i0 to i1 do
        for k from k0 to k1 do
            for j from j0 to j1 do
                c(i, k) + a(i, j)*b(j, k) -> c(i, k);
            endfor;
        endfor;
    endfor;
enddefine;