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

なお、上記コードは個人的な勉強のために作成。
（コピペして、機能を適宜追加し、ご自由に使ってください）