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:


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

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

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.

...float) Math.sqrt(a [N - 1, j]) ;

Bryan Carpenter 2003-04-15