diff --git a/examples/exact_values_example_1D.py b/examples/exact_values_example_1D.py
new file mode 100644
index 0000000..5ceeadc
--- /dev/null
+++ b/examples/exact_values_example_1D.py
@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+"""
+Exact Values
+=================
+
+PyKrige demonstration and usage
+as a non-exact interpolator in 1D. 
+"""
+
+from pykrige.ok import OrdinaryKriging
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+np.random.seed(42)
+
+x = np.linspace(0, 12.5, 50)
+xpred = np.linspace(0, 12.5, 393)
+y = np.sin(x) * np.exp(-0.25 * x) + np.random.normal(-0.25, 0.25, 50)
+
+# compare OrdinaryKriging as an exact and non exact interpolator
+uk = OrdinaryKriging(
+    x, np.zeros(x.shape), y, variogram_model="linear", exact_values=False
+)
+uk_exact = OrdinaryKriging(x, np.zeros(x.shape), y, variogram_model="linear")
+
+y_pred, y_std = uk.execute("grid", xpred, np.array([0.0]), backend="loop")
+y_pred_exact, y_std_exact = uk_exact.execute(
+    "grid", xpred, np.array([0.0]), backend="loop"
+)
+
+
+y_pred = np.squeeze(y_pred)
+y_std = np.squeeze(y_std)
+
+y_pred_exact = np.squeeze(y_pred_exact)
+y_std_exact = np.squeeze(y_std_exact)
+
+
+fig, ax = plt.subplots(1, 1, figsize=(10, 4))
+
+ax.scatter(x, y, label="Input Data")
+ax.plot(xpred, y_pred_exact, label="Exact Prediction")
+ax.plot(xpred, y_pred, label="Non Exact Prediction")
+
+ax.fill_between(
+    xpred,
+    y_pred - 3 * y_std,
+    y_pred + 3 * y_std,
+    alpha=0.3,
+    label="Confidence interval",
+)
+ax.legend(loc=9)
+ax.set_ylim(-1.8, 1.3)
+ax.legend(loc=9)
+plt.xlabel("X")
+plt.ylabel("Field")
+plt.show()
diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx
index be15dbe..056c544 100644
--- a/pykrige/lib/cok.pyx
+++ b/pykrige/lib/cok.pyx
@@ -48,7 +48,8 @@ cpdef _c_exec_loop(double [:, ::1] a_all,
             b[k] *= -1
         b[n] = 1.0
 
-        check_b_vect(n, bd, b, eps)
+        if not pars['exact_values']:
+            check_b_vect(n, bd, b, eps)
 
 
         # Do the BLAS matrix-vector multiplication call
diff --git a/pykrige/ok.py b/pykrige/ok.py
index 6935ac3..b34e4c5 100644
--- a/pykrige/ok.py
+++ b/pykrige/ok.py
@@ -15,6 +15,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
+.. [2] N. Cressie, Statistics for spatial data, 
+   (Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
 """
@@ -140,11 +142,21 @@ class OrdinaryKriging:
         coordinates, x is interpreted as longitude and y as latitude
         coordinates, both given in degree. Longitudes are expected in
         [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'.
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     eps = 1.0e-10  # Cutoff for comparison to zero
@@ -173,11 +185,17 @@ def __init__(
         enable_plotting=False,
         enable_statistics=False,
         coordinates_type="euclidean",
+        exact_values=True,
     ):
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -585,11 +603,11 @@ def _get_kriging_matrix(self, n):
             )
         a = np.zeros((n + 1, n + 1))
         a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d)
+
         np.fill_diagonal(a, 0.0)
         a[n, :] = 1.0
         a[:, n] = 1.0
         a[n, n] = 0.0
-
         return a
 
     def _exec_vector(self, a, bd, mask):
@@ -609,7 +627,7 @@ def _exec_vector(self, a, bd, mask):
 
         b = np.zeros((npt, n + 1, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
         b[:, n, 0] = 1.0
 
@@ -647,7 +665,7 @@ def _exec_loop(self, a, bd_all, mask):
 
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
             x = np.dot(a_inv, b)
@@ -683,7 +701,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx):
                 zero_value = False
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 
@@ -705,25 +723,6 @@ def execute(
     ):
         """Calculates a kriged grid and the associated variance.
 
-        This is now the method that performs the main kriging calculation.
-        Note that currently measurements (i.e., z values) are considered
-        'exact'. This means that, when a specified coordinate for interpolation
-        is exactly the same as one of the data points, the variogram evaluated
-        at the point is forced to be zero. Also, the diagonal of the kriging
-        matrix is also always forced to be zero. In forcing the variogram
-        evaluated at data points to be zero, we are effectively saying that
-        there is no variance at that point (no uncertainty,
-        so the value is 'exact').
-
-        In the future, the code may include an extra 'exact_values' boolean
-        flag that can be adjusted to specify whether to treat the measurements
-        as 'exact'. Setting the flag to false would indicate that the
-        variogram should not be forced to be zero at zero distance
-        (i.e., when evaluated at data points). Instead, the uncertainty in
-        the point will be equal to the nugget. This would mean that the
-        diagonal of the kriging matrix would be set to
-        the nugget instead of to zero.
-
         Parameters
         ----------
         style : str
@@ -876,6 +875,7 @@ def execute(
                     "eps",
                     "variogram_model_parameters",
                     "variogram_function",
+                    "exact_values",
                 ]
             }
 
diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py
index 679d2c9..396f0a3 100644
--- a/pykrige/ok3d.py
+++ b/pykrige/ok3d.py
@@ -14,6 +14,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
+.. [2] N. Cressie, Statistics for spatial data, 
+(Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
 """
@@ -151,11 +153,21 @@ class OrdinaryKriging3D:
         Default is False (off).
     enable_plotting : bool, optional
         Enables plotting to display variogram. Default is False (off).
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     eps = 1.0e-10  # Cutoff for comparison to zero
@@ -186,11 +198,17 @@ def __init__(
         anisotropy_angle_z=0.0,
         verbose=False,
         enable_plotting=False,
+        exact_values=True,
     ):
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -594,7 +612,7 @@ def _exec_vector(self, a, bd, mask):
 
         b = np.zeros((npt, n + 1, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
         b[:, n, 0] = 1.0
 
@@ -632,7 +650,7 @@ def _exec_loop(self, a, bd_all, mask):
 
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 
@@ -669,7 +687,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx):
                 zero_index = None
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 
diff --git a/pykrige/uk.py b/pykrige/uk.py
index 1e42e00..7858e17 100644
--- a/pykrige/uk.py
+++ b/pykrige/uk.py
@@ -15,7 +15,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
-
+.. [2] N. Cressie, Statistics for spatial data, 
+   (Wiley Series in Probability and Statistics, 1993) 137 p.
 Copyright (c) 2015-2020, PyKrige Developers
 """
 import numpy as np
@@ -172,11 +173,21 @@ class UniversalKriging:
         Default is False (off).
     enable_plotting : boolean, optional
         Enables plotting to display variogram. Default is False (off).
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     UNBIAS = True  # This can be changed to remove the unbiasedness condition
@@ -212,6 +223,7 @@ def __init__(
         functional_drift=None,
         verbose=False,
         enable_plotting=False,
+        exact_values=True,
     ):
 
         # Deal with mutable default argument
@@ -225,6 +237,11 @@ def __init__(
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -820,6 +837,7 @@ def _get_kriging_matrix(self, n, n_withdrifts):
         else:
             a = np.zeros((n_withdrifts, n_withdrifts))
         a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d)
+
         np.fill_diagonal(a, 0.0)
 
         i = n
@@ -888,7 +906,7 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids)
         else:
             b = np.zeros((npt, n_withdrifts, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
 
         i = n
@@ -979,7 +997,7 @@ def _exec_loop(self, a, bd_all, xy, xy_orig, mask, n_withdrifts, spec_drift_grid
             else:
                 b = np.zeros((n_withdrifts, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
 
             i = n
@@ -1039,24 +1057,6 @@ def execute(
         """Calculates a kriged grid and the associated variance.
         Includes drift terms.
 
-        This is now the method that performs the main kriging calculation.
-        Note that currently measurements (i.e., z values) are considered
-        'exact'. This means that, when a specified coordinate for interpolation
-        is exactly the same as one of the data points, the variogram evaluated
-        at the point is forced to be zero. Also, the diagonal of the kriging
-        matrix is also always forced to be zero. In forcing the variogram
-        evaluated at data points to be zero, we are effectively saying that
-        there is no variance at that point (no uncertainty,
-        so the value is 'exact').
-
-        In the future, the code may include an extra 'exact_values' boolean
-        flag that can be adjusted to specify whether to treat the measurements
-        as 'exact'. Setting the flag to false would indicate that the variogram
-        should not be forced to be zero at zero distance (i.e., when evaluated
-        at data points). Instead, the uncertainty in the point will be equal to
-        the nugget. This would mean that the diagonal of the kriging matrix
-        would be set to the nugget instead of to zero.
-
         Parameters
         ----------
         style : str
diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py
index c94664a..b78d6ca 100644
--- a/pykrige/uk3d.py
+++ b/pykrige/uk3d.py
@@ -14,6 +14,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
+.. [2] N. Cressie, Statistics for spatial data, 
+   (Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
 """
@@ -166,11 +168,21 @@ class UniversalKriging3D:
         Default is False (off).
     enable_plotting : boolean, optional
         Enables plotting to display variogram. Default is False (off).
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     UNBIAS = True  # This can be changed to remove the unbiasedness condition
@@ -206,6 +218,7 @@ def __init__(
         functional_drift=None,
         verbose=False,
         enable_plotting=False,
+        exact_values=True,
     ):
 
         # Deal with mutable default argument
@@ -219,6 +232,11 @@ def __init__(
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -712,7 +730,7 @@ def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids):
         else:
             b = np.zeros((npt, n_withdrifts, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
 
         i = n
@@ -788,7 +806,7 @@ def _exec_loop(self, a, bd_all, xyz, mask, n_withdrifts, spec_drift_grids):
             else:
                 b = np.zeros((n_withdrifts, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
 
             i = n
diff --git a/tests/test_core.py b/tests/test_core.py
index d17038b..5a13b04 100755
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -436,6 +436,61 @@ def test_core_krige_3d():
     assert ss == approx(0.0, rel=1e-3)
 
 
+def test_non_exact():
+    # custom data for this test
+    data = np.array(
+        [[0.0, 0.0, 0.47], [1.5, 1.5, 0.56], [3, 3, 0.74], [4.5, 4.5, 1.47],]
+    )
+
+    # construct grid points so diagonal
+    # is identical to input points
+    gridx = np.arange(0.0, 4.51, 1.5)
+    gridy = np.arange(0.0, 4.51, 1.5)
+
+    ok = OrdinaryKriging(
+        data[:, 0],
+        data[:, 1],
+        data[:, 2],
+        variogram_model="exponential",
+        variogram_parameters=[500.0, 3000.0, 5.0],
+    )
+    z, ss = ok.execute("grid", gridx, gridy, backend="vectorized")
+
+    ok_non_exact = OrdinaryKriging(
+        data[:, 0],
+        data[:, 1],
+        data[:, 2],
+        variogram_model="exponential",
+        variogram_parameters=[500.0, 3000.0, 5.0],
+        exact_values=False,
+    )
+    z_non_exact, ss_non_exact = ok_non_exact.execute(
+        "grid", gridx, gridy, backend="vectorized"
+    )
+
+    in_values = np.diag(z)
+
+    # test that krig field
+    # at input location are identical
+    # to the inputs themselves  with
+    # exact_values == True
+    assert_allclose(in_values, data[:, 2])
+
+    # test that krig field
+    # at input location are different
+    # than the inputs themselves
+    # with exact_values == False
+    assert ~np.allclose(in_values, data[:, 2])
+
+    # test that off diagonal values are the same
+    # by filling with dummy value and comparing
+    # each entry in array
+    np.fill_diagonal(z, 0.0)
+    np.fill_diagonal(z_non_exact, 0.0)
+
+    assert_allclose(z, z_non_exact)
+
+
 def test_ok(validation_ref):
 
     # Test to compare OK results to those obtained using KT3D_H2O.
@@ -519,6 +574,13 @@ def test_ok_execute(sample_data_2d):
 
     ok = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2])
 
+    with pytest.raises(ValueError):
+        OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values="blurg")
+
+    ok_non_exact = OrdinaryKriging(
+        data[:, 0], data[:, 1], data[:, 2], exact_values=False
+    )
+
     with pytest.raises(ValueError):
         ok.execute("blurg", gridx, gridy)
 
@@ -538,6 +600,18 @@ def test_ok_execute(sample_data_2d):
     assert np.amax(ss) != np.amin(ss)
     assert not np.ma.is_masked(z)
 
+    z1, ss1 = ok_non_exact.execute("grid", gridx, gridy, backend="loop")
+    assert_allclose(z1, z)
+    assert_allclose(ss1, ss)
+
+    z, ss = ok_non_exact.execute("grid", gridx, gridy, backend="loop")
+    shape = (gridy.size, gridx.size)
+    assert z.shape == shape
+    assert ss.shape == shape
+    assert np.amax(z) != np.amin(z)
+    assert np.amax(ss) != np.amin(ss)
+    assert not np.ma.is_masked(z)
+
     with pytest.raises(IOError):
         ok.execute("masked", gridx, gridy, backend="vectorized")
     mask = np.array([True, False])
@@ -570,6 +644,14 @@ def test_ok_execute(sample_data_2d):
     assert z[0, 0] is np.ma.masked
     assert ss[0, 0] is np.ma.masked
 
+    z, ss = ok_non_exact.execute(
+        "masked", gridx, gridy, mask=mask_ref.T, backend="loop"
+    )
+    assert np.ma.is_masked(z)
+    assert np.ma.is_masked(ss)
+    assert z[0, 0] is np.ma.masked
+    assert ss[0, 0] is np.ma.masked
+
     with pytest.raises(ValueError):
         ok.execute(
             "points",
@@ -594,11 +676,20 @@ def test_cython_ok(sample_data_2d):
     data, (gridx, gridy, _), mask_ref = sample_data_2d
 
     ok = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2])
+    ok_non_exact = OrdinaryKriging(
+        data[:, 0], data[:, 1], data[:, 2], exact_values=False
+    )
+
     z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop")
     z2, ss2 = ok.execute("grid", gridx, gridy, backend="C")
     assert_allclose(z1, z2)
     assert_allclose(ss1, ss2)
 
+    z1, ss1 = ok_non_exact.execute("grid", gridx, gridy, backend="loop")
+    z2, ss2 = ok_non_exact.execute("grid", gridx, gridy, backend="C")
+    assert_allclose(z1, z2)
+    assert_allclose(ss1, ss2)
+
     closest_points = 4
 
     z1, ss1 = ok.execute(
@@ -610,6 +701,15 @@ def test_cython_ok(sample_data_2d):
     assert_allclose(z1, z2)
     assert_allclose(ss1, ss2)
 
+    z1, ss1 = ok_non_exact.execute(
+        "grid", gridx, gridy, backend="loop", n_closest_points=closest_points
+    )
+    z2, ss2 = ok_non_exact.execute(
+        "grid", gridx, gridy, backend="C", n_closest_points=closest_points
+    )
+    assert_allclose(z1, z2)
+    assert_allclose(ss1, ss2)
+
 
 def test_uk(validation_ref):
 
@@ -818,6 +918,24 @@ def test_uk_execute(sample_data_2d):
         drift_terms=["regional_linear"],
     )
 
+    with pytest.raises(ValueError):
+        UniversalKriging(
+            data[:, 0],
+            data[:, 1],
+            data[:, 2],
+            variogram_model="linear",
+            drift_terms=["regional_linear"],
+            exact_values="blurg",
+        )
+
+    uk_non_exact = UniversalKriging(
+        data[:, 0],
+        data[:, 1],
+        data[:, 2],
+        variogram_model="linear",
+        drift_terms=["regional_linear"],
+    )
+
     with pytest.raises(ValueError):
         uk.execute("blurg", gridx, gridy)
     with pytest.raises(ValueError):
@@ -831,6 +949,18 @@ def test_uk_execute(sample_data_2d):
     assert np.amax(ss) != np.amin(ss)
     assert not np.ma.is_masked(z)
 
+    z1, ss1 = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized")
+    assert_allclose(z1, z)
+    assert_allclose(ss1, ss)
+
+    z, ss = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized")
+    shape = (gridy.size, gridx.size)
+    assert z.shape == shape
+    assert ss.shape == shape
+    assert np.amax(z) != np.amin(z)
+    assert np.amax(ss) != np.amin(ss)
+    assert not np.ma.is_masked(z)
+
     z, ss = uk.execute("grid", gridx, gridy, backend="loop")
     shape = (gridy.size, gridx.size)
     assert z.shape == shape
@@ -871,6 +1001,14 @@ def test_uk_execute(sample_data_2d):
     assert z[0, 0] is np.ma.masked
     assert ss[0, 0] is np.ma.masked
 
+    z, ss = uk_non_exact.execute(
+        "masked", gridx, gridy, mask=mask_ref.T, backend="loop"
+    )
+    assert np.ma.is_masked(z)
+    assert np.ma.is_masked(ss)
+    assert z[0, 0] is np.ma.masked
+    assert ss[0, 0] is np.ma.masked
+
     with pytest.raises(ValueError):
         uk.execute(
             "points",
@@ -914,15 +1052,32 @@ def test_ok_uk_produce_same_result(validation_ref):
         verbose=False,
         enable_plotting=False,
     )
+    uk_non_exact = UniversalKriging(
+        data[:, 0],
+        data[:, 1],
+        data[:, 2],
+        variogram_model="linear",
+        verbose=False,
+        enable_plotting=False,
+        exact_values=False,
+    )
     z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="vectorized")
     assert_allclose(z_ok, z_uk)
     assert_allclose(ss_ok, ss_uk)
 
+    z_uk, ss_uk = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized")
+    assert_allclose(z_ok, z_uk)
+    assert_allclose(ss_ok, ss_uk)
+
     z_ok, ss_ok = ok.execute("grid", gridx, gridy, backend="loop")
     z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="loop")
     assert_allclose(z_ok, z_uk)
     assert_allclose(ss_ok, ss_uk)
 
+    z_uk, ss_uk = uk_non_exact.execute("grid", gridx, gridy, backend="loop")
+    assert_allclose(z_ok, z_uk)
+    assert_allclose(ss_ok, ss_uk)
+
 
 def test_ok_backends_produce_same_result(validation_ref):
 
@@ -1742,6 +1897,28 @@ def test_ok3d(validation_ref):
         variogram_model="exponential",
         variogram_parameters=[500.0, 3000.0, 0.0],
     )
+
+    with pytest.raises(ValueError):
+        OrdinaryKriging3D(
+            data[:, 0],
+            data[:, 1],
+            np.zeros(data[:, 1].shape),
+            data[:, 2],
+            variogram_model="exponential",
+            variogram_parameters=[500.0, 3000.0, 0.0],
+            exact_values="blurg",
+        )
+
+    ok3d_non_exact = OrdinaryKriging3D(
+        data[:, 0],
+        data[:, 1],
+        np.zeros(data[:, 1].shape),
+        data[:, 2],
+        variogram_model="exponential",
+        variogram_parameters=[500.0, 3000.0, 0.0],
+        exact_values=False,
+    )
+
     k, ss = k3d.execute(
         "grid", gridx_ref, gridy_ref, np.array([0.0]), backend="vectorized"
     )
@@ -1972,6 +2149,17 @@ def test_ok3d_backends_produce_same_result(sample_data_3d):
     k3d = OrdinaryKriging3D(
         data[:, 0], data[:, 1], data[:, 2], data[:, 3], variogram_model="linear"
     )
+
+    ok3d_non_exact = OrdinaryKriging3D(
+        data[:, 0],
+        data[:, 1],
+        np.zeros(data[:, 1].shape),
+        data[:, 2],
+        variogram_model="exponential",
+        variogram_parameters=[500.0, 3000.0, 0.0],
+        exact_values=False,
+    )
+
     k_k3d_v, ss_k3d_v = k3d.execute(
         "grid", gridx_ref, gridy_ref, gridz_ref, backend="vectorized"
     )
@@ -1981,6 +2169,15 @@ def test_ok3d_backends_produce_same_result(sample_data_3d):
     assert_allclose(k_k3d_v, k_k3d_l, rtol=1e-05, atol=1e-8)
     assert_allclose(ss_k3d_v, ss_k3d_l, rtol=1e-05, atol=1e-8)
 
+    k, ss = ok3d_non_exact.execute(
+        "grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="loop"
+    )
+    k1, ss1 = ok3d_non_exact.execute(
+        "grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="vectorized"
+    )
+    assert_allclose(k1, k)
+    assert_allclose(ss1, ss)
+
 
 def test_ok3d_execute(sample_data_3d):