ハライドでのコレスキー分解

Question

私はハライドでコレスキー分解を実装しようとしています。 croutのような一般的なアルゴリズムの一部は、三角マトリックス上の反復で構成されています。分解の対角要素は、入力行列の対角要素から部分列和を減算することによって計算されます。列合計は、対角要素を除いて、入力行列の三角形部分の二乗要素にわたって計算されます。このようなパターンはハライドで表現できる一般的である場合

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

私の質問は：？次のようになり++コードは、Cの場合とBLASを使用して

もしそうなら、どのように見えるでしょうか？

Answer 1

私はAndrewがより完全な答えを持っているかもしれないと思っていますが、適時の応答のために、RDom述語（RDom :: whereから導入）を使用して三角領域を列挙します。パターンのスケッチである。

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

その後減少triangularを使用します。

Answer 2

私はかつてHalideで書かれた速いコレスキーを持っていました。残念ながら私はコードを見つけることができません。私は外側のループをCに入れて、一度に32ワイドパネルのようなもので動作する良いブロックパネル更新ルーチンを書いた。これは、Halideが三角反復をする前であったので、おそらくあなたは今より良くできます。