Recall complexity to solve $N\times N$ linear system from the Poisson equation:
Given a domain $\Omega$. We represent it is as a union of two domain $\Omega_1$ and $\Omega_2$ with intersection $\Gamma$.
Then, the problem
$$ \Delta u = f, \quad u_{\partial \Omega} = g $$is equivalent to solving two equations
1) $ \Delta u_1 = f_1, \quad u_1|_{\partial \Omega_1} = g|_{ \partial \Omega_1} $
2) $ \Delta u_2 = f_2, \quad u_2|_{\partial \Omega_2} = g|_{ \partial \Omega_2} $
with additional constraints
$$u_1 \mid_{\Gamma} = u_2 \mid_\Gamma,$$and also continuity of the solution on the intersection $\Gamma$:
$$ \left.\frac{\partial u_1}{\partial n} \right|_{\Gamma} = \left.\frac{\partial u_2}{\partial n} \right|_{\Gamma}.$$Can be formally written as an equation on the interface with unknown $\lambda = u \mid_{\Gamma}$: $$ S \lambda = \chi, $$ which is called Poincare-Steklov equation.
The simplest iterative method for the domain decomposition is the so-called Neumann-Dirichlet iterative method.
We successively solve the subproblems: $$-\Delta u^{(k+1)}_1 = f, \quad \mbox{in } \Omega_1, \quad u^{(k+1)}_1 = 0, \quad \mbox{in } \partial \Omega_1 \setminus \Gamma, \quad u^{(k+1)}_1 = \lambda^k \quad \mbox{in } \Gamma.$$
$$ -\Delta u^{(k+1)}_2 = f, \quad \mbox{in } \Omega_2, \quad u^{(k+1)}_1 = 0, \quad \mbox{in } \partial \Omega_2 \setminus \Gamma, \quad \frac{\partial u^{(k+1)}_2}{\partial n_2} = -\frac{\partial u^{(k+1)}_1}{\partial n_1}, \quad \mbox{in } \Gamma. $$And then $$ \lambda^{k+1} = \theta u_2^{(k+1)}\mid_{\Gamma} + (1-\theta) \lambda_k $$
Domain decomposition methods are understood in the best way when treated algebraically.
In the matrix form, the system has the form
$$ \begin{bmatrix} A_{11} & A_{1 \Gamma} & 0 \\ A_{\Gamma 1} & A_{\Gamma} & A_{\Gamma 2} \\ 0 & A_{2 \Gamma} & A_{22} \end{bmatrix}\begin{bmatrix} u_1 \\ u_{\Gamma} \\ u_2 \end{bmatrix} = \begin{bmatrix} f_1 \\ f_{\Gamma} \\ f_2 \end{bmatrix}. $$In FEM, the interface matrix $A_{\Gamma}$ is formed from the contributions from two domains:
$$A_{\Gamma} = A^{(1)}_{\Gamma} + A^{(2)}_{\Gamma}.$$If we eliminate $u_1$ and $u_2$ using standard block Gaussian elimination, we have the following equation on the interface:
$$\left(A_{\Gamma} - A_{\Gamma 1} A_{11}^{-1} A_{1 \Gamma} - A_{\Gamma 2} A_{22}^{-1} A_{2 \Gamma}\right) u_{\Gamma} = f_{\Gamma} - A_{\Gamma} A^{-1}_{11} f_1 - A_{\Gamma 2} A^{-1}_{22} f_2. $$The Schur complement in the continuous form is the Poincare-Steklov operator.
The Neumann-Dirichlet iteration is the simple iteration (Richardson method) applied to the Schur complement:
$$\lambda^{k+1} = \lambda^k + \theta S^{-1}_2 (- S \lambda^k + \chi).$$In order to multiply $S^{-1}_2$ by vector, it is sufficient to solve linear system
$$ \begin{bmatrix} A^{(2)}_{\Gamma} & A_{\Gamma 2} \\ A_{2 \Gamma} & A_{22} \end{bmatrix} \begin{bmatrix} \eta_{\Gamma} \\ \eta_2 \end{bmatrix} = \begin{bmatrix} \chi_{\Gamma} \\ 0 \end{bmatrix}. $$And that is a Laplace equation with Neumann boundary condition!
In case of multiple domains, similar construct can be written down, especially for the case of strips.
The Schur complement will be the sum of $S_i$, where $S_i$ is the interface matrix of the form
$$S_i = N_i (A^{(i)}_{\Gamma} - A_{\Gamma i} A_{ii}^{-1} A_{i \Gamma}) N^{\top}_i.$$Consider the matrix $B_u$ of the form
$$B_u = \begin{bmatrix} A_{11} & A_{1 \Gamma} & 0 \\ A_{\Gamma 1} & A^{(1)}_{\Gamma} & 0 \\ 0 & A_{2 \Gamma} & A_{22} \end{bmatrix} $$This is a non-symmetric matrix, however on the specific subspace
$$U = \{ A_{2 \Gamma} x_{\Gamma} + A_{22} x_2 = 0 \}$$it is a symmetric-positive definite and spectrally equivalent to $A$.
Therefore, if at each iteration step the error belongs to the subspace $U$ (which is the case), then $B_u$ can be an efficient preconditioner.
The symmetric variant has the form
$$B_s = \begin{bmatrix} A_{11} & A_{1 \Gamma} & 0 \\ A_{\Gamma 1} & A^{(1)}_{\Gamma}+ A_{\Gamma 2} A^{-1}_{22} A_{2 \Gamma} & A_{\Gamma 2} \\ 0 & A_{2 \Gamma} & A_{22} \end{bmatrix} $$The solution of a linear system with the matrix $B_s$ requires the solution of linear systems with matrices $A_{22}$, $B_1$, $A_{22}$.
Moreover, $B_s$ is spectrally equivalent to the matrix $A$, if $A = A^{\top}$.
i.e.
$$C_{21} A \leq B_s \leq A.$$The main disadvantage of $B_s$ is that the subdomain problems have to be solved exactly.
Neumann-Neumann method can be generalized to the multiple-domain case.
The preconditioner has the form
$$B_s = \sum_{i=1}^m N_i D_i S^{-1}_i D_i N^{\top}_i,$$where $D_i$ is a diagonal matrix that contains the inverse numbers of the domains, containing the particular node, $N_i$ is the local assembly matrix
(i.e. $\sum_{i=1}^m N_i D_i N^{\top}_i = I.$).
The convergence estimate by Mandel given
$$\mathrm{cond}(B_S S) \leq C\left(1 + \log \frac{d}{h} \right)\max \left\{\frac{C}{d^2}; \log \frac{d}{h} \right\}.$$In practice, the dependence on the number of subdomains is smaller, and can be removed by solving an additional problem on the grid that is created by the subdomain subdivision.
Let $R_i$ be the projection matrix onto the indices, corresponding to the domain (remember the subspace correction algorithm).
Then,
$$A_i = R_i A R^{\top}_i,$$and one step can be written as
$$A_i z_i = R_i (f - Au), \quad u := u + R_i z_i, \quad i = 1, \ldots, m.$$Then the transition matrix has the form
$$e^{n+1} = (I - P_m)\ldots (I - P_1) e^n,$$where $$P_i = R_i A^{-1}_i R_i^{\top}.$$ The matrix $P_i$ is a projector (let us show that!).