テンプレートによる多重ループの静的展開

#include <cstddef>
#include <iostream>

template<
  std::size_t I, std::size_t L, std::size_t K, std::size_t M
, std::size_t J, std::size_t N
>
struct j_loop
{
  static void invoke(
    double const (&a)[L][M], double const (&b)[M][N], double (&c)[L][N]
  )
  {
    c[I][J] += a[I][K] * b[K][J];
    j_loop<I, L, K, M, J + 1, N>::invoke(a, b, c);
  }
};

template<
  std::size_t I, std::size_t L, std::size_t K, std::size_t M, std::size_t N
>
struct j_loop<I, L, K, M, N, N>
{
  static void invoke(
    double const (&)[L][M], double const (&)[M][N], double (&)[L][N]
  )
  {}
};

template<
  std::size_t I, std::size_t L, std::size_t K, std::size_t M, std::size_t N
>
struct k_loop
{
  static void invoke(
    double const (&a)[L][M], double const (&b)[M][N], double (&c)[L][N]
  )
  {
    j_loop<I, L, K, M, 0, N>::invoke(a, b, c);
    k_loop<I, L, K + 1, M, N>::invoke(a, b, c);
  }
};

template<std::size_t I, std::size_t L, std::size_t M, std::size_t N>
struct k_loop<I, L, M, M, N>
{
  static void invoke(
    double const (&)[L][M], double const (&)[M][N], double (&)[L][N]
  )
  {}
};

template<std::size_t I, std::size_t M, std::size_t J, std::size_t N>
struct j_init_loop
{
  static void invoke(double (&a)[M][N])
  {
    a[I][J] = 0.0;
    j_init_loop<I, M, J + 1, N>::invoke(a);
  }
};

template<std::size_t I, std::size_t M, std::size_t N>
struct j_init_loop<I, M, N, N>
{
  static void invoke(double (&)[M][N])
  {}
};

template<std::size_t I, std::size_t L, std::size_t M, std::size_t N>
struct i_loop
{
  static void invoke(
    double const (&a)[L][M], double const (&b)[M][N], double (&c)[L][N]
  )
  {
    j_init_loop<I, L, 0, N>::invoke(c);
    k_loop<I, L, 0, M, N>::invoke(a, b, c);
    i_loop<I + 1, L, M, N>::invoke(a, b, c);
  }
};

template<std::size_t L, std::size_t M, std::size_t N>
struct i_loop<L, L, M, N>
{
  static void invoke(
    double const (&)[L][M], double const (&)[M][N], double (&)[L][N]
  )
  {}
};

template<std::size_t L, std::size_t M, std::size_t N>
void static_matrix_multiplication(
  double const (&a)[L][M], double const (&b)[M][N], double (&c)[L][N]
)
{
  i_loop<0, L, M, N>::invoke(a, b, c);
}

using namespace std;

int main()
{
  double a[3][3] =
  {
    6.0, 1.0, 8.0,
    7.0, 5.0, 3.0,
    2.0, 9.0, 4.0
  };

  double b[3][3] =
  {
    2.0, 7.0, 6.0,
    9.0, 5.0, 1.0,
    4.0, 3.0, 8.0
  };

  double c[3][3];

  static_matrix_multiplication(a, b, c);
  // 以下のコードに相当するコードを静的に展開する．
  // for(std::size_t i = 0; i < 3; ++i){
  //   for(std::size_t j = 0; j < 3; ++j){
  //     c[i][j] = 0.0;
  //   }
  //   for(std::size_t k = 0; k < 3; ++k){
  //     for(std::size_t j = 0; j < 3; ++j){
  //       c[i][j] += a[i][k] * b[k][j];
  //     }
  //   }
  // }

  cout << c[0][0] << " " << c[0][1] << " " << c[0][2] << endl;
  cout << c[1][0] << " " << c[1][1] << " " << c[1][2] << endl;
  cout << c[2][0] << " " << c[2][1] << " " << c[2][2] << endl;
  
  return 0;
}

再帰関数テンプレートを用いてベクトルの内積などの静的展開をやって見せているコードはどこにでもあるけれど，上でやっているような多重ループの静的展開は見たことがない．何でだろう？
テンプレートによる多重ループの静的展開で問題になるのは関数テンプレートを部分特殊化できない，ということ．上ではクラステンプレートの静的メンバ関数によって関数テンプレートの部分特殊化をエミュレートするテクニックを使っている．
元々はテンプレート引数で大きさを指定するような小型の行列での計算の静的展開や，中規模以上の行列計算のアンロールに使うことを想定して考えたやつだったり．
ところで上のループはikjだけれど静的展開でもikjが速いのかな？直感的にはikjが一番速そうなんだけれど・・・．
上のコードVC++7.1の最適化だとcの値が即値になりますな．恐ろしや．