docs: def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids):
          """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)
  
          if self.UNBIAS:
              b = np.zeros((npt, n_withdrifts+1, 1))
          else:
              b = np.zeros((npt, n_withdrifts, 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
  
          i = n
          if self.regional_linear_drift:
              b[:, i, 0] = xyz[:, 2]
              i += 1
              b[:, i, 0] = xyz[:, 1]
              i += 1
              b[:, i, 0] = xyz[:, 0]
              i += 1
          if self.specified_drift:
              for spec_vals in spec_drift_grids:
                  b[:, i, 0] = spec_vals.flatten()
                  i += 1
          if self.functional_drift:
              for func in self.functional_drift_terms:
                  b[:, i, 0] = func(xyz[:, 2], xyz[:, 1], xyz[:, 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
  
          if (~mask).any():
              mask_b = np.repeat(mask[:, np.newaxis, np.newaxis],
                                 n_withdrifts+1, axis=1)
              b = np.ma.array(b, mask=mask_b)
  
          if self.UNBIAS:
              x = np.dot(a_inv, b.reshape((npt, n_withdrifts+1)).T).reshape((1, n_withdrifts+1, npt)).T
          else:
              x = np.dot(a_inv, b.reshape((npt, n_withdrifts)).T).reshape((1, n_withdrifts, npt)).T
          kvalues = np.sum(x[:, :n, 0] * self.VALUES, axis=1)
          sigmasq = np.sum(x[:, :, 0] * -b[:, :, 0], axis=1)
  
          return kvalues, sigmasq