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

Cholesky decomposition

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

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

where $L$ is a lower triangular matrix. In practice 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 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{displaymath}
\begin{minipage}[t]{\linewidth}\small\begin{verbatim}Proc...
...float) Math.sqrt(a [N - 1, j]) ;
}\end{verbatim}\end{minipage}\end{displaymath}



Bryan Carpenter 2003-04-15