  Linear systems with a tensor product structure arise naturally when
considering the discretization of Laplace type differential equations or, more
generally, multidimensional operators with separable coefficients. In this
work, we focus on the numerical solution of linear systems of the form $$
\left(I\otimes \dots\otimes I \otimes A_1+\dots + A_d\otimes I \otimes\dots
\otimes I\right)x=b,$$ where the matrices $A_t\in\mathbb R^{n\times n}$ are
symmetric positive definite and belong to the class of hierarchically
semiseparable matrices.
  We propose and analyze a nested divide-and-conquer scheme, based on the
technology of low-rank updates, that attains the quasi-optimal computational
cost $\mathcal O(n^d (\log(n) + \log(\kappa)^2 + \log(\kappa)
\log(\epsilon^{-1})))$ where $\kappa$ is the condition number of the linear
system, and $\epsilon$ the target accuracy. Our theoretical analysis highlights
the role of inexactness in the nested calls of our algorithm and provides worst
case estimates for the amplification of the residual norm. The performances are
validated on 2D and 3D case studies.
