Solving a complex, sparse, linear system using the Schur complement

Question

Solution method

I am repetitively solving sparse linear systems (for the need of ARNOLDI iterations) of the type: $$\underbrace{\begin{bmatrix} J_1 & J_2 \\ J_3 & J_4 \end{bmatrix}}_J \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_z= \underbrace{\begin{bmatrix} b_x \\ b_y \end{bmatrix}}_b$$

I build the Schur-complement ($J_4$ is non-singular): $$J_D=J_1-J_2J_4^{-1}J_3$$ Then, I solve for $x$: $$J_Dx=b_x-J_2J_4^{-1}b_y$$ and then for $y$: $$J_4y=b_y-J_3x$$ At this point I have a full solution of $z$.

Intermediate calculations

To calculate $J_2J_4^{-1}J_3$, I first compute $\tilde{J_2}=J_2J_4^{-1}$ as follows: $$\tilde{J_2}^T=(J_2J_4^{-1})^T=(J_4^{-1})^TJ_2^T$$ Thus, I get $\tilde{J_2}$ by solving the system: $$J_4^T\tilde{J_2}^T=J_2^T$$ and tacking the transpose.

Edit: The above is implemented in Matlab and I use the conjugate transpose ' (mathworks.com/help/matlab/ref/ctranspose.html) for both the real and the complex. Also I use dot (mathworks.com/help/matlab/ref/dot.html) for the dot products.

Problem

This way of solving works perfectly when $J$ is real (even if $b$ is complex). That is, solving $Jz=b$ directly or by building the Schur-complement gives the same $z$ as expected. However, when $J_4$ becomes complex, then I get different results... I suspect that I did something really wrong when taking the transpose for $\tilde{J_2}$ or doing the dot product for $J_2J_4^{-1}b_y$. Yet, I fail to detect my error.

Note: The reason I solve using the Schur-complement is that matrices $J_2,J_3,J_4$ have some particular structures that allow for good parallelization of the operations involved. E.g. $J_4$ is block diagonal and non-singular.

Answer 1

Your approach looks good, so you may have a bug somewhere in your implementation. As a validation, I ran the following simple test case involving a randomly generated dense complex matrix $J$ and a randomly generated real right-hand side $b$.

N = 100;
M = 30;
J = exp(2i*pi*rand(N));
b = rand(N,1);
Jd = (J(M+1:end,M+1:end)' \ J(1:M,M+1:end)')';
x = b(1:M) - Jd*b(M+1:end);
Jd = J(1:M,1:M) - Jd*J(M+1:end,1:M);
x = Jd \ x;
y = b(M+1:end) - J(M+1:end,1:M)*x;
y = J(M+1:end,M+1:end) \ y;
z = [x; y];
ztrue = J \ b;
err = norm(ztrue - z) / norm(ztrue);
res = norm(J*z - b) / norm(b);
restrue = norm(J*ztrue - b) / norm(b);

Running this code multiple times the approximate averaged values for $\texttt{err}$, $\texttt{res}$, and $\texttt{restrue}$ were $1\times 10^{-13}$, $1\times 10^{-13}$, and $4\times 10^{-14}$, respectively.