code: def solvedbi_sm(ah, rho, b, c=None, axis=4):
      r
  
      a = np.conj(ah)
      if c is None:
          c = solvedbi_sm_c(ah, a, rho, axis)
      if have_numexpr:
          cb = inner(c, b, axis=axis)
          return ne.evaluate('(b - (a * cb)) / rho')
      else:
          return (b - (a * inner(c, b, axis=axis))) / rho
docs: r"""
      Solve a diagonal block linear system with a scaled identity term
      using the Sherman-Morrison equation.
  
      The solution is obtained by independently solving a set of linear
      systems of the form (see :cite:`wohlberg-2016-efficient`)
  
      .. math::
        (\rho I + \mathbf{a} \mathbf{a}^H ) \; \mathbf{x} = \mathbf{b} \;\;.
  
      In this equation inner products and matrix products are taken along
      the specified axis of the corresponding multi-dimensional arrays; the
      solutions are independent over the other axes.
  
      Parameters
      ----------
      ah : array_like
        Linear system component :math:`\mathbf{a}^H`
      rho : float
        Linear system parameter :math:`\rho`
      b : array_like
        Linear system component :math:`\mathbf{b}`
      c : array_like, optional (default None)
        Solution component :math:`\mathbf{c}` that may be pre-computed using
        :func:`solvedbi_sm_c` and cached for re-use.
      axis : int, optional (default 4)
        Axis along which to solve the linear system
  
      Returns
      -------
      x : ndarray
        Linear system solution :math:`\mathbf{x}`