uBLAS::compressed_matrixに直接データを指定する方法 #C++

以下の様な疎行列

matrix.txt

6, 2, 0, 1, 4,
2, 7, 0, 0, 2,
0, 0, 6, 3, 0,
1, 0, 3, 9, 0,
4, 2, 0, 0, 8,

を、uBLASのCRS形式(boost::numeric::ublas::compressd_matrix)の疎行列

crs.cpp

// 行列の次元数
const int COUNT = 5;

// 非ゼロ要素数
const int NONZERO_COUNT = 15;

// 行列の生成
boost::numeric::ublas::compressed_matrix<double> mat(COUNT, COUNT);

として格納したい場合、普通だと

slow.cpp

mat(0, 0) = 6; mat(0, 1) = 2;                mat(0, 3) = 1; mat(0, 4) = 4;
mat(1, 0) = 2; mat(1, 1) = 7;                               mat(1, 4) = 2;
                              mat(2, 2) = 6; mat(2, 3) = 3;
mat(3, 0) = 1;                mat(3, 2) = 3; mat(3, 3) = 9;
mat(4, 0) = 4; mat(4, 1) = 2;                               mat(4, 4) = 8;

ってな感じで、mat(i, j) = valueを使う方法が多い。
というか、どこをみてもたいていこれ以外の方法が乗ってない。

が、これって実はかなり遅い。特に、元からCRS形式のデータが

crs_data.cpp

// 要素値
double elements[NONZERO_COUNT] ={
6, 2,    1, 4,
2, 7,       2,
      6, 3,   
1,    3, 9,   
4, 2,       8, 
};
// 列番号
unsigned int columnIndeces[NONZERO_COUNT] ={
0, 1,    3, 4,
0, 1,       4,
      2, 3,   
1,    2, 3,   
0, 1,       4, 
};
// 各行先頭の位置
unsigned int rowOffsets[COUNT+1] ={
0,
4,
7,
9,
12,
15 //＝非ゼロ要素数 ※ここ忘れがち
};

という風に用意されていて他と共有する場合は、これらのデータを直接流し込むことができる。

fast.cpp

// 非ゼロ要素数の容量を増やす
mat.reserve(NONZERO_COUNT, false);

// データを複製
std::copy(elements,      elements      + NONZERO_COUNT, mat.value_data().begin());
std::copy(columnIndeces, columnIndeces + NONZERO_COUNT, mat.index2_data().begin());
std::copy(rowOffsets,    rowOffsets    + COUNT + 1,     mat.index1_data().begin());

// 要素数を設定
mat.set_filled(COUNT+1, NONZERO_COUNT);

解説すると、まず、CRS形式データそのものは、data、index1、index2に格納されている。index1とindex2ってなんだって思ったら、1が列番号で、2が各行先頭の位置を保持しているらしい。というわけで、それらを取得してstd::copyで複製してやるというのがこの方法。

ただしそれだけではダメで、まず最初にreserveを使って非ゼロ要素数が格納できる容量を増やしてやる必要がある。reserveの引数#2はデフォルト(true)でも良いが、連結する必要はないのでfalseで。

更に、データを複製をした後（前だとダメ！）に、実際の要素数をset_filledで指定する。引数がfilled1とfilled2になっていて「なんのことだ」と思ったら、これがどうやら、各行先頭位置配列の大きさ（つまり非ゼロ要素数）と列番号配列の大きさ（つまり次元数+1）に対応しているらしい。要するに、先に書いたのと同じく、1が列番号で2が各行先頭位置。

なぜかこれをきちんと解説したサイトがなかった（調べきれてないだけかもしれないけど）ので、ここに書いときます。