3

我正在尝试在 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 是否可以普遍表达这样的模式？如果是这样，它会是什么样子？

4

2 回答 2

1

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

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

然后triangular在减少使用。

于 2016-12-06T07:01:31.827 回答
0

我曾经有一个用 Halide 写的快速 Cholesky。不幸的是我找不到代码。我将外部循环放在 C 中，并编写了一个很好的块面板更新例程，该例程一次在类似于 32 宽的面板上运行。这是在 Halide 进行三角迭代之前，所以也许你现在可以做得更好。

于 2016-12-07T17:34:44.923 回答