Cholesky decomposition

If $A$ is a symmetric positive definite matrix, associated linear equations are often solved using Choleski decomposition:

$$A=LL{^T}$$

where $L$ is a lower triangular matrix. In practice this is followed by forward and back substitutions:

$$\mathit{L}\mathbf{y}=\mathbf{b}, \; \mathit{L{^T}}\mathbf{x}=\mathbf{y}$$

to complete the solution of $A\mathbf{x}=b$. A pseudocode algorithm for Cholesky decomposition is

$$\begin{array}{l} For\; k=1\; to\; n-1 \\ \;\;\;\; l_{kk}=a_... ...ij}=a_{ij}-l_{ik}l_{jk} \\ l_{nn}=a_{nn}^{1/2} \\ \end{array}$$

A parallel version, assuming the main array is stored by columns with the rows cyclically distributed, is given in figure 4.4. The $l$ array is accumulated in the lower part of the input array a. Note that the array m has a replicated distribution, so the remap operation is a broadcast of the relevant part of column $k$.

Figure 4.4: Cholesky decomposition.

$$\begin{minipage}[t]{\linewidth}\small\begin{verbatim}Proc... ...float) Math.sqrt(a [N - 1, j]) ; }\end{verbatim}\end{minipage}$$