Cholesky decomposition in Halide
<p>I'm trying to implement a Cholesky decomposition in Halide. Part of common algorithm such as crout consists of an iteration over a triangular matrix. In a way that, the diagonal elements of the decomposition are computed by subtracting a partial column sum from the diagonal element of the input matrix. Column sum is calculated over squared elements of a triangular part of the input matrix, excluding the diagonal element.</p>

<p>Using BLAS the code would in C++ look as follows:</p>

<pre><code>double* a; /* input matrix */
int n; /* dimension */
const int c__1 = 1;
const double c_b12 = 1.;
const double c_b10 = -1.;

for (int j = 0; j &lt; n; ++j) {
    double ajj = a[j + j * n] - ddot(&amp;j, &amp;a[j + n], &amp;n, &amp;a[j + n], &amp;n);
    ajj = sqrt(ajj);
    a[j + j * n] = ajj;
    if (j &lt; n) {
            int i__2 = n - j;
            dgemv("No transpose", &amp;i__2, &amp;j, &amp;c_b10, &amp;a[j + 1 + n], &amp;n, &amp;a[j + n], &amp;b, &amp;c_b12, &amp;a[j + 1 + j * n], &amp;c__1);
            double d__1 = 1. / ajj;
            dscal(&amp;i__2, &amp;d__1, &amp;a[j + 1 + j * n], &amp;c__1);            
    }
}
</code></pre>

<p>My question is if a pattern like this is in general expressible by Halide? And if so, how would it look like?</p>