BLASとかで高速化が難しい線形計算

SNL の計算でいまボトルネックとなっている部分が BLAS などでは高速化できないようなので手こずっている。

計算は数学的には簡単で、3 x n 行列 X に対して n 次元ベクトル a があって、このときに 3 x n 行列 Y を求めるものである。ただし、Y(i,j)  = X(i,j) * a(j) である。

Matlab の場合は、
Y = X*spdiags(a,0,length(a),length(a));
でできるにはできるが、これがかなり遅い。

計算式がY(i,j) = X(i,j) *a(j)であるため、各 j について BLAS の daxpy を使うこともできるのだが、iのループが3しかなく、BLAS を呼ぶオーバーヘッドが大きくなり BLAS も効率的に使うことができない。

現状は
 for j=1:n
  val = a(j);
  Y(1,j) = X(1,j) * val;
  Y(2,j) = X(2,j) * val;
  Y(3,j) = X(3,j) * val;
end
でループ展開をして回していて、上記の BLAS, Matlab よりも高速になっているが、それでもボトルネックのままなので、高速化をしたいところである。

なお、n=200万ぐらい必要で、この計算が 3000 回程度必要となっている。

今日の作業内容：Matlab 高速化 3h + 論文読み 3h
今日のランチ：らく 鶏照り焼き冷麦
明日の予測作業時間：2h