docs: def _exec_loop(self, a, bd_all, xyz, mask, n_withdrifts, spec_drift_grids):
          """Solves the kriging system by looping over all specified points.
          Less memory-intensive, but involves a Python-level loop."""
  
          npt = bd_all.shape[0]
          n = self.X_ADJUSTED.shape[0]
          kvalues = np.zeros(npt)
          sigmasq = np.zeros(npt)
  
          a_inv = scipy.linalg.inv(a)
  
          for j in np.nonzero(~mask)[0]:   # Note that this is the same thing as range(npt) if mask is not defined,
              bd = bd_all[j]               # otherwise it takes the non-masked elements.
              if np.any(np.absolute(bd) <= self.eps):
                  zero_value = True
                  zero_index = np.where(np.absolute(bd) <= self.eps)
              else:
                  zero_value = False
                  zero_index = None
  
              if self.UNBIAS:
                  b = np.zeros((n_withdrifts+1, 1))
              else:
                  b = np.zeros((n_withdrifts, 1))
              b[:n, 0] = - self.variogram_function(self.variogram_model_parameters, bd)
              if zero_value:
                  b[zero_index[0], 0] = 0.0
  
              i = n
              if self.regional_linear_drift:
                  b[i, 0] = xyz[j, 2]
                  i += 1
                  b[i, 0] = xyz[j, 1]
                  i += 1
                  b[i, 0] = xyz[j, 0]
                  i += 1
              if self.specified_drift:
                  for spec_vals in spec_drift_grids:
                      b[i, 0] = spec_vals.flatten()[i]
                      i += 1
              if self.functional_drift:
                  for func in self.functional_drift_terms:
                      b[i, 0] = func(xyz[j, 2], xyz[j, 1], xyz[j, 0])
                      i += 1
              if i != n_withdrifts:
                  warnings.warn("Error in setting up kriging system. "
                                "Kriging may fail.", RuntimeWarning)
              if self.UNBIAS:
                  b[n_withdrifts, 0] = 1.0
  
              x = np.dot(a_inv, b)
              kvalues[j] = np.sum(x[:n, 0] * self.VALUES)
              sigmasq[j] = np.sum(x[:, 0] * -b[:, 0])
  
          return kvalues, sigmasq