def METHOD_NAME(input_dim: int, output_dim: int, N: int, Ntest: int) -> None:
    covar_shape = (output_dim, Ntest, Ntest)
    X, Y = rng.randn(N, input_dim), rng.randn(N, output_dim)
    Xtest = rng.randn(Ntest, input_dim)
    kernel = Matern32()
    model_gp = gpflow.models.GPR((X, Y), kernel=kernel)
    mu1, var = model_gp.predict_f(Xtest, full_cov=False)
    mu2, covar = model_gp.predict_f(Xtest, full_cov=True)
    assert np.allclose(mu1, mu2, atol=1.0e-10)
    assert covar.shape == covar_shape
    assert var.shape == (Ntest, output_dim)
    for i in range(output_dim):
        assert np.allclose(var[:, i], np.diag(covar[i, :, :]))