$$
\newcommand{\RR}{\mathbb{R}}
\newcommand{\GG}{\mathbb{G}}
\newcommand{\PP}{\mathbb{P}}
\newcommand{\PS}{\mathcal{P}}
\newcommand{\SS}{\mathbb{S}}
\newcommand{\NN}{\mathbb{N}}
\newcommand{\ZZ}{\mathbb{Z}}
\newcommand{\CC}{\mathbb{C}}
\newcommand{\HH}{\mathbb{H}}
\newcommand{\ones}{\mathbb{1\hspace{-0.4em}1}}
\newcommand{\alg}[1]{\mathfrak{#1}}
\newcommand{\mat}[1]{ \begin{pmatrix} #1 \end{pmatrix} }
\renewcommand{\bar}{\overline}
\renewcommand{\hat}{\widehat}
\renewcommand{\tilde}{\widetilde}
\newcommand{\inv}[1]{ {#1}^{-1} }
\newcommand{\eqdef}{\overset{\text{def}}=}
\newcommand{\block}[1]{\left(#1\right)}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\trace}[1]{\mathrm{tr}\block{#1}}
\newcommand{\norm}[1]{ \left\| #1 \right\| }
\newcommand{\argmin}[1]{ \underset{#1}{\mathrm{argmin}} }
\newcommand{\argmax}[1]{ \underset{#1}{\mathrm{argmax}} }
\newcommand{\st}{\ \mathrm{s.t.}\ }
\newcommand{\sign}[1]{\mathrm{sign}\block{#1}}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\inner}[1]{\langle #1 \rangle}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\ddd}[2]{\frac{\partial #1}{\partial #2} }
\newcommand{\db}{\dd^b}
\newcommand{\ds}{\dd^s}
\newcommand{\dL}{\dd_L}
\newcommand{\dR}{\dd_R}
\newcommand{\Ad}{\mathrm{Ad}}
\newcommand{\ad}{\mathrm{ad}}
\newcommand{\LL}{\mathcal{L}}
\newcommand{\Krylov}{\mathcal{K}}
\newcommand{\Span}[1]{\mathrm{Span}\block{#1}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\tr}{\mathrm{tr}}
\newcommand{\sinc}{\mathrm{sinc}}
\newcommand{\cat}[1]{\mathcal{#1}}
\newcommand{\Ob}[1]{\mathrm{Ob}\block{\cat{#1}}}
\newcommand{\Hom}[1]{\mathrm{Hom}\block{\cat{#1}}}
\newcommand{\op}[1]{\cat{#1}^{op}}
\newcommand{\hom}[2]{\cat{#1}\block{#2}}
\newcommand{\id}{\mathrm{id}}
\newcommand{\Set}{\mathbb{Set}}
\newcommand{\Cat}{\mathbb{Cat}}
\newcommand{\Hask}{\mathbb{Hask}}
\newcommand{\lim}{\mathrm{lim}\ }
\newcommand{\funcat}[1]{\left[\cat{#1}\right]}
\newcommand{\natsq}[6]{
\begin{matrix}
& #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\
{#1}_{#4} \hspace{-1.5em} &\downarrow & & \downarrow & \hspace{-1.5em} {#1}_{#5}\\
& #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\
\end{matrix}
}
\newcommand{\comtri}[6]{
\begin{matrix}
#1 & \overset{#4}\longrightarrow & #2 & \\
#6 \hspace{-1em} & \searrow & \downarrow & \hspace{-1em} #5 \\
& & #3 &
\end{matrix}
}
\newcommand{\natism}[6]{
\begin{matrix}
& #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\
{#1}_{#4} \hspace{-1.5em} &\downarrow \uparrow & & \downarrow \uparrow & \hspace{-1.5em} {#1}_{#5}\\
& #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\
\end{matrix}
}
\newcommand{\cone}[1]{\mathcal{#1}}
$$
Lower-Diagonal Solves
Let \(M = L + D + L^T\) be a symmetric matrix with positive diagonal
\(D\). We look at the various in-place lower-diagonal system solves
found e.g. in Gauss-Seidel method.
Solving \(\block{L + D}x = y\)
This gives:
\[\sum_{j < i} M_{ij} x_j + M_{ii} x_i = y_i\]
Thus, solving incrementally for \(x_i\) gives:
\[x_i = \frac{1}{M_{ii}} \block{y_i - \sum_{j < i} M_{ij} x_j}\]
which can also be performed in-place on vector \(y\):
\[y_i := \frac{1}{M_{ii}} \block{y_i - \sum_{j < i} M_{ij} y_j}\]
Solving \(\block{L + D}x = -L^T y\)
This gives:
\[\sum_{j < i} M_{ij} x_j + M_{ii} x_i = -\sum_{j > i} M_{ij} y_j\]
Thus, solving incrementally for \(x_i\) gives:
\[x_i = \frac{1}{M_{ii}} \block{- \sum_{j > i} M_{ij} y_j - \sum_{j < i} M_{ij} x_j}\]
which can also be performed in-place on vector \(y\):
\[y_i := \frac{1}{M_{ii}} \block{- \sum_{j > i} M_{ij} y_j - \sum_{j < i} M_{ij} y_j}\]
Alternatively:
\[y_i := y_i - \frac{1}{M_{ii}} \sum_j M_{ij} y_j\]
Solving \(\block{L + D}x = -L^T y - z\)
Mixing the above gives, for all \(i\):
\[\sum_{j < i} M_{ij} x_j + M_{ii} x_i = -\sum_{j > i} M_{ij} y_j - z_i\]
Thus, solving incrementally for \(x_i\) gives:
\[x_i = \frac{1}{M_{ii}} \block{- \sum_{j > i} M_{ij} y_j - \sum_{j < i} M_{ij} x_j - z_i}\]
which can also be performed in-place on vector \(y\):
\[y_i := \frac{1}{M_{ii}} \block{- \sum_{j > i} M_{ij} y_j - \sum_{j < i} M_{ij} y_j - z_i}\]
Alternatively:
\[y_i := y_i - \frac{1}{M_{ii}} \block{z_i + \sum_j M_{ij} y_j}\]
This one is used in Gauss-Seidel method:
\[\begin{align}
x_{k+1} &= x_k - \block{L + D}^{-1}\block{Mx_k + q} \\
&= - \block{L + D}^{-1} \block{L^T x_k + q} \\
\end{align}\]
e-mail comments/suggestions