def METHOD_NAME(a: te.Tensor, b: te.Tensor) -> te.Tensor:
    a_shape = list(a.shape)
    b_shape = list(b.shape)
    a_prepended = False
    b_appended = False
    if len(a_shape) == 1:
        a_prepended = True
        a_shape.insert(0, 1)
    if len(b_shape) == 1:
        b_appended = True
        b_shape.append(1)
    is_a_larger = len(a_shape) > len(b_shape)
    offset = len(a_shape) - len(b_shape) if is_a_larger else len(b_shape) - len(a_shape)
    a_relax = relax.Var("a", relax.TensorStructInfo(a.shape))
    b_relax = relax.Var("b", relax.TensorStructInfo(b.shape))
    f_infer_sinfo = call.op.get_attr("FInferStructInfo")
    output_shape = f_infer_sinfo(relax.op.matmul(a_relax, b_relax), bb).shape
    def matmul_compute(*idx_spatial):
        k = te.reduce_axis((0, a_shape[-1]), name="k")
        def multiply_compute(idx_reduce):
            a_indices = []
            b_indices = []
            for i in range(offset):
                if is_a_larger:
                    a_indices.append(idx_spatial[i])
                else:
                    b_indices.append(idx_spatial[i])
            for i in range(offset, len(output_shape) - (2 - a_prepended - b_appended)):
                a_dim = a_shape[i if is_a_larger else i - offset]
                b_dim = b_shape[i if not is_a_larger else i - offset]
                dim_equal = a_dim == b_dim
                if not isinstance(dim_equal, tir.IntImm) or dim_equal == 0:
                    a_dim_is_one = isinstance(a_dim, tir.IntImm) and a_dim == 1
                    b_dim_is_one = isinstance(b_dim, tir.IntImm) and b_dim == 1
                    a_indices.append(0 if a_dim_is_one else idx_spatial[i])
                    b_indices.append(0 if b_dim_is_one else idx_spatial[i])
                else:
                    a_indices.append(idx_spatial[i])
                    b_indices.append(idx_spatial[i])
            if not a_prepended:
                a_indices.append(idx_spatial[-2 + b_appended])
            a_indices.append(idx_reduce)
            b_indices.append(idx_reduce)
            if not b_appended:
                b_indices.append(idx_spatial[-1])
            dtype = call.attrs.out_dtype
            if dtype != "":
                return a(*a_indices).astype(dtype) * b(*b_indices).astype(dtype)
            return a(*a_indices) * b(*b_indices)
        return te.sum(multiply_compute(k), axis=k)
    return te.compute(
        output_shape,
        lambda *idx: matmul_compute(*idx),  # pylint: disable=unnecessary-lambda
        name="matmul",
    )