docs: def _exec_vector(self, a, bd, mask):
          """Solves the kriging system as a vectorized operation. This method
          can take a lot of memory for large grids and/or large datasets."""
  
          npt = bd.shape[0]
          n = self.X_ADJUSTED.shape[0]
          zero_index = None
          zero_value = False
  
          a_inv = scipy.linalg.inv(a)
  
          if np.any(np.absolute(bd) <= self.eps):
              zero_value = True
              zero_index = np.where(np.absolute(bd) <= self.eps)
  
          b = np.zeros((npt, n+1, 1))
          b[:, :n, 0] = - self.variogram_function(self.variogram_model_parameters, bd)
          if zero_value:
              b[zero_index[0], zero_index[1], 0] = 0.0
          b[:, n, 0] = 1.0
  
          if (~mask).any():
              mask_b = np.repeat(mask[:, np.newaxis, np.newaxis], n+1, axis=1)
              b = np.ma.array(b, mask=mask_b)
  
          x = np.dot(a_inv, b.reshape((npt, n+1)).T).reshape((1, n+1, npt)).T
          zvalues = np.sum(x[:, :n, 0] * self.Z, axis=1)
          sigmasq = np.sum(x[:, :, 0] * -b[:, :, 0], axis=1)
  
          return zvalues, sigmasq