CUDAカーネル内で行列やベクトルの演算が必要になった際、コードが複雑になりそうだったため、
以下の部品を作ってみた。動的配列を使いたかったが、メモリ確保に数十ミリ秒を要したため、
静的配列を用いた。このとき、サイズ変更に柔軟に対応するために、サイズをテンプレート引数
として与えるようにした。
いまのところ、メンバ関数は、転置と、掛け算と、3x3の逆行列で、いずれもCUDAカーネル内から
呼ぶことができる。
template<typename T, unsigned int w, unsigned int h> struct K_mat {
T val[h][w];
__device__ void transpose(const K_mat<T, h, w>& in) {
for (unsigned int i = 0; i < h; i++) {
for (unsigned int j = 0; j < w; j++) {
val[i][j] = in.val[j][i];
}
}
}
template<unsigned int d> __device__ void mul(const K_mat<T, d, h>& in0, const K_mat<T, w, d>& in1) {
for (unsigned int i = 0; i < h; i++) {
for (unsigned int j = 0; j < w; j++) {
val[i][j] = 0.0;
for (unsigned int k = 0; k < d; k++) {
val[i][j] = val[i][j] + in0.val[i][k] * in1.val[k][j];
}
}
}
}
__device__ void inverse_3x3(const K_mat<T, h, w>& in) {
float det = abs(in.val[0][0] * in.val[1][1] * in.val[2][2]
+ in.val[0][1] * in.val[1][2] * in.val[2][0]
+ in.val[0][2] * in.val[2][1] * in.val[1][0]
- in.val[0][2] * in.val[1][1] * in.val[2][0]
- in.val[0][0] * in.val[2][1] * in.val[1][2]
- in.val[0][1] * in.val[1][0] * in.val[2][2]);
if (1e-7 < det) { // machine epsilon of 32bit is 1.192e-7
for (unsigned int i = 0; i < h; i++) {
for (unsigned int j = 0; j < w; j++) {
val[i][j] = in.val[(i + 1) % 3][(j + 1) % 3] * in.val[(i + 2) % 3][(j + 2) % 3]
- in.val[(i + 2) % 3][(j + 1) % 3] * in.val[(i + 1) % 3][(j + 2) % 3];
val[i][j] = __fdividef(val[i][j], det);
}
}
}
else {
for (unsigned int i = 0; i < h; i++) {
for (unsigned int j = 0; j < w; j++) {
val[i][j] = 0.0;
}
}
}
}
};
なお、上記コードは個人的な勉強のために作成。
（コピペして、機能を適宜追加し、ご自由に使ってください）