c++ - Halide 中的胆固醇分解

我正在尝试在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表示吗？如果是这样，它将是什么样子？

最佳答案

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

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

然后在还原中使用triangular。