def matmul(C: BufferizedNDArray, A: BufferizedNDArray, B: BufferizedNDArray) -> BufferizedNDArray:
    A_: BufferizedNDArray = A
    A__val_slot: np_buf_t(float64) = unpack(A_.val)
    B_: BufferizedNDArray = B
    B__val_slot: np_buf_t(float64) = unpack(B_.val)
    C_: BufferizedNDArray = C
    C__val_slot: np_buf_t(float64) = unpack(C_.val)
    m: finch.int64 = A_.shape.element_0
    n: finch.int64 = B_.shape.element_1
    p: finch.int64 = A_.shape.element_1
    for i in range(0, length(slot(C__val_slot, np_buf_t(float64)))):
        store(slot(C__val_slot, np_buf_t(float64)), i, 0.0)
    for i in range(0, m):
        i__pos: finch.int64 = add(0, mul(A_.strides.element_0, i))
        i__pos_2: finch.int64 = add(0, mul(C_.strides.element_0, i))
        for k in range(0, p):
            k__pos: finch.int64 = add(i__pos, mul(A_.strides.element_1, k))
            k__pos_2: finch.int64 = add(0, mul(B_.strides.element_0, k))
            for j in range(0, n):
                j__pos: finch.int64 = add(k__pos_2, mul(B_.strides.element_1, j))
                j__pos_2: finch.int64 = add(i__pos_2, mul(C_.strides.element_1, j))
                a_ik: finch.float64 = load(slot(A__val_slot, np_buf_t(float64)), k__pos)
                b_kj: finch.float64 = load(slot(B__val_slot, np_buf_t(float64)), j__pos)
                c_ij: finch.float64 = mul(a_ik, b_kj)
                store(slot(C__val_slot, np_buf_t(float64)), j__pos_2, add(load(slot(C__val_slot, np_buf_t(float64)), j__pos_2), c_ij))
    repack(C__val_slot)
    matmul_return: BufferizedNDArray = C
    return matmul_return
