next up previous contents
Next: Matrix multiplication with reduced Up: Array Sections Previous: Two-dimensional Fourier transform   Contents

Cholesky decomposition

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

\begin{displaymath}
A=LL{^T}
\end{displaymath}

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

\begin{displaymath}
\mathit{L}\mathbf{y}=\mathbf{b}, \; \mathit{L{^T}}\mathbf{x}=\mathbf{y}
\end{displaymath}

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

\begin{displaymath}
\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}\end{displaymath}

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 b has a replicated distribution, so the remap operation is a broadcast of the relevant part of column $k$.

Figure 4.4: Choleksy decomposition.
\begin{figure}\small\begin{verbatim}Procs1 p = new Procs1(P) ;
on(p) {
Rang...
... [N - 1])
a [N - 1, j] = Math.sqrt(a [N - 1, j]) ;
}\end{verbatim}\end{figure}



Bryan Carpenter 2002-07-12