2016-12-05 42 views
3

我正試圖在Halide中實現Cholesky分解。部分常見算法如crout包含三角矩陣上的迭代。以某種方式,分解的對角線元素是通過從輸入矩陣的對角線元素中減去部分列和來計算的。列總和是在輸入矩陣的三角形部分的平方元素上計算的,不包括對角線元素。Halides中的Cholesky分解

使用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);    
    } 
} 

我的問題是,如果這樣的模式是一般表達的鹵化物?如果是這樣,它會是什麼樣子?

+0

您是否最終在Halide中實現了這個功能? – zanbri

回答

1

我認爲安德魯可能有一個更完整的答案,但爲了及時響應,您可以使用RDom謂詞(通過RDom :: where引入)來枚舉三角形區域(或將它們推廣到更多維度)。模式的草圖是:

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

然後在減少使用triangular

0

我曾經在哈立德寫過一本快速的Cholesky。不幸的是我找不到代碼。我把外層循環放在C中,並且寫了一個很好的塊面板更新例程,這個例程一次操作32像素的面板。這是在Halide進行三角迭代之前,所以你現在可以做得更好。