c++ - Halide 中的 Cholesky 分解

Question

我正在尝试在 Halide 中实现 Cholesky 分解。部分常见算法(例如 crout)由对三角矩阵的迭代组成。在某种程度上，分解的对角线元素是通过从输入矩阵的对角线元素中减去部分列和来计算的。列和是在输入矩阵的三角形部分的平方元素上计算的，不包括对角线元素。

使用 BLAS，C++ 代码如下所示:

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 < n; ++j) {
    double ajj = a[j + j * n] - ddot(&j, &a[j + n], &n, &a[j + n], &n);
    ajj = sqrt(ajj);
    a[j + j * n] = ajj;
    if (j < n) {
            int i__2 = n - j;
            dgemv("No transpose", &i__2, &j, &c_b10, &a[j + 1 + n], &n, &a[j + n], &b, &c_b12, &a[j + 1 + j * n], &c__1);
            double d__1 = 1. / ajj;
            dscal(&i__2, &d__1, &a[j + 1 + j * n], &c__1);            
    }
}

我的问题是这样的模式是否通常可以用 Halide 表达？如果是这样，它会是什么样子？

Answer 1

我认为 Andrew 可能有更完整的答案，但为了及时响应，您可以使用 RDom 谓词(通过 RDom::where 引入)来枚举三角形区域(或将它们推广到更多维度)。图案的草图是:

Halide::RDom triangular(0, extent, 0, extent);
triangular.where(triangular.x < triangular.y);

然后在缩减中使用triangular。