\documentstyle{article}

\begin{document}

\title{Choleski decomposition in HPJava}
%\title{Choleski decomposition in Java${^\#}$}

\author{Guansong Zhang, Bryan Carpenter}

%\institute{NPAC at Syracuse University, Syracuse, NY 13244}

%\date{Sep. 12, 1997}

\maketitle

\section{Introduction}

If $A$ is symmetric positive definite, then an often used algorithm
for direct solving linear equation is Choleski decomposition.

\[A=LL{^T}\]

where $L$ is a lower triangular matrix, followed by the 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 $A\mathbf{x}=b$.

In \cite{Ort89} different kinds of algorithm organization are
discussed.  Among them, column interleaved \emph{outer product} form
and row interleaved \emph{column Choleski} form are considered more
promising ones on a parallel computer system with distributed memory.

\section{Outer product form}

The outer product form of Choleski algorithm is based on the following
mathematic equation:

\begin{displaymath}
A_1=
\left[ \begin{array}{cc}
a_{11}           & \mathbf{a}_1\mathit{^T}   \\
\mathbf{a}_{1}   & A_s                       \\
\end{array} \right]
=
\left[ \begin{array}{cc}
a_{11}                   & a_{11}^{1/2}\mathbf{l}_1\mathit{^T}      \\
a_{11}^{1/2}\mathbf{l}_1 & A_2+ \mathbf{l}_1l_1\mathit{^T}          \\
\end{array} \right]
\end{displaymath}
\begin{displaymath}
=
\left[ \begin{array}{cc}
a_{11}^{1/2}    &  0  \\
\mathbf{l}_{1}  &  I  \\
\end{array} \right]
\left[ \begin{array}{cc}
1 & 0    \\
0 & A_2  \\
\end{array} \right]
\left[ \begin{array}{cc}
a_{11}^{1/2} & \mathbf{l}_1\mathit{^T}   \\
   0         &  I                        \\
\end{array} \right]
\end{displaymath}

the updating of the submatix $A_2$ is accomplished by subtracting the
outer product of $\mathbf{l}_1l_1\mathit{^T}$.

A simple description of the algorithm in a sequential code may be,

\begin{displaymath}
\begin{array}{l}
l_{11}=a_{11}^{1/2} \\
For\; k=1\; to\; n-1 \\
\;\;\;\; For\; s=k+1\; to\; n \\
\;\;\;\;\;\;\;\; l_{sk}=a_{sk}/l_{kk} \\
\;\;\;\;  For\; j=k+1\; to\; n \\
\;\;\;\;\;\;\;\;    For\; i=j\; to\; n \\
\;\;\;\;\;\;\;\;\;\;\;\;      a_{ij}=a_{ij}-l_{ik}l_{jk} \\
\;\;\;\; l_{k+1, k+1}= a_{k+1, k+1}^{1/2} \\
\end{array}
\end{displaymath}

On a parallel computer with a distributed memory, the array may have a
column-interleaved storage, then each processor can update its own
column after geting the first updated column of the matrix or
submatrix.

We write the above algorithm in the following HPJava program. In the
code, a cyclic range is used to allocate the orignal array on
multiprocessors, the result $L$ can overwrite part of its storage.
During the computation, collective communication \texttt{remap} is
used for broadcasting updated columns.

\small
\begin{verbatim}
public class ChKjiC {
  public static void main() {

    final int size=100;
    Range x = new CyclicRange(size, id);

    float a[[*,]] = new float [[size, x]];
    // initialize the array;

    float b[[*]] = new float [[size]]; // used as a buffer

    at (x|0) 
      a[[0,x]]=sqrt(a[[0,x]]);

    for (int k=0; k<size-1; k++) {

      for (int s=k+1; s<size; s++)
        at (x|k)
          a[[s,x]]/=a[[x,x]];

      Adlib.remap(b[k+1:],a[k,x|k+1:]);

      over (x|k+1:) 
        for (int i=x; i<size; i++)
          a[[i,x]]-=b[[i]]*b[[x]];

      at (x|k+1)
        a[[k+1,x]]=sqrt(a[[k+1,x]]);
    }
  }
}
\end{verbatim}
\normalsize

\section{Column Choleski form}

This form is based on the following:

if $A_{i-1}=L_{i-1}L_{i-1}^T$ is the decomposition of the $(i-1)st$
principal submatrix, then the decomposition of $A_i$ is

\begin{displaymath}
A_i=
\left[ \begin{array}{cc}
A_{i-1}          & \mathbf{a}_i   \\
\mathbf{a}_i\mathit{^T}   & a_{ii}                    \\
\end{array} \right]
=
\left[ \begin{array}{cc}
L_{i-1} & 0  \\
\mathbf{l}_i\mathit{^T} & l_{ii}  \\
\end{array} \right]
\left[ \begin{array}{cc}
L_{i-1}^{T} & \mathbf{l}_i  \\
0 & l_{ii}  \\
\end{array} \right]
\end{displaymath}

Thus, the new $ith$ row of $L$ is obtained from
$L_{i-1}\mathbf{l}_i=\mathbf{a}_{i}$,
$l_{ii}=(a_{ii}-\mathbf{l}_i\mathit{^T}\mathbf{l}_i)\mathit{^{1/2}}$.

A description of this form of algorithm in a sequential code can be,

\begin{displaymath}
\begin{array}{l}
l_{11}=a_{11}^{1/2} \\
For\; j=2\; to\; n \\
\;\;\;\; For\; s=j\; to\; n \\
\;\;\;\;\;\;\;\; l_{sj-1}=a_{sj-1}/l_{j-1,j-1} \\
\;\;\;\;  For\; k=1\; to\; j-1 \\
\;\;\;\;\;\;\;\;    For\; i=j\; to\; n \\
\;\;\;\;\;\;\;\;\;\;\;\;      a_{ij}=a_{ij}-l_{ik}l_{jk} \\
\;\;\;\;\;\;\;\; l_{jj}= a_{jj}^{1/2} \\
\end{array}
\end{displaymath}

Written in HPJava with row-interleaved storage, the program should be,

\small
\begin{verbatim}
public class ChJkiR {
  public static void main() {

    final int size=100;

    Range x = new CyclicRange(size, id);

    float a[[,*]] = new float [[x, size]];

    // init array ;
    
    float b[[*]] = new float [[size]];

    at (x|0)
      a[[x,0]]=sqrt(a[[x,0]]);

    for (int j=1; j<size; j++) {
      
      Adlib.remap(b[:j-1], a[x|j-1,:j-1]);

      over(x|j:)
        a[[x, j-1]]/=b[[j-1]];

      for (int k=0; k<j-1; k++) {
        over (x|j:)
          a[[x,j]]-=a[[x,k]]*b[[k]];
        
        at (x|j)
          a[[x,j]]=sqrt(a[[x,j]]);
      }
    }
  }
}
\end{verbatim}
\normalsize

\section{Performance}

There are three major loops in the above algorithms, with each index
as \texttt{i}, \texttt{j} and \texttt{k}. After proper arrangement, we
may get 6 different orders of the combination of these loops. Then in
each permutation, we can either use row interleaved storage or column
interleaved storage. So there will be total 12 different oganizition
of the Choleski decomposition. The above two are considered the best.

testing

(to be added)

%I like the first one, which is easy to understand.

\end{document}
