C++でI×J行列 × J×K行列

  //! @brief IxJ行列の i,j の要素にアクセスする
  //! @param [in] i 行番号
  //! @param [in] j 列番号
  //! @param [in] I 行数
  inline int matijI(const int i, const int j, const int I) {
    return i + I * j;
  }

  //! @brief IxJ行列 × JxK行列 の積
  //! @param [out] C Cik 計算結果格納先
  //! @param [in] A Aijの行列
  //! @param [in] B Bjkの行列
  //! @param [in] I Aの行列の行数
  //! @param [in] J Aの行列の列数
  //! @param [in] K Bの行列の列数
  //! @return C
  template<
    typename real_tC,
    typename real_tA,
    typename real_tB>
  real_tC* mult_matrixIJK(real_tC* C, const real_tA* A, const real_tB* B, const int I, const int J, const int K) {

    for (int i = 0; i < I; i++)
      for (int k = 0; k < K; k++) {
        C[matijI(i, k, I)] = 0.0;
        for (int j = 0; j < J; j++) {
          C[matijI(i, k, I)] += A[matijI(i, j, I)] * B[matijI(j, k, J)];
        }
      }

    return C;
  }

#include <Eigen/Dense>

void mult_matrix_check() {

  std::cout << "---------- Eigen ------------" << std::endl;
  Eigen::MatrixXf mat53(5,3);
  mat53(0, 0) = 00; mat53(0, 1) = 01; mat53(0, 2) = 02;
  mat53(1, 0) = 10; mat53(1, 1) = 11; mat53(1, 2) = 12;
  mat53(2, 0) = 20; mat53(2, 1) = 21; mat53(2, 2) = 22;
  mat53(3, 0) = 30; mat53(3, 1) = 31; mat53(3, 2) = 32;
  mat53(4, 0) = 40; mat53(4, 1) = 41; mat53(4, 2) = 42;

  Eigen::MatrixXf mat34(3, 4);
  mat34(0, 0) = 100; mat34(0, 1) = 101; mat34(0, 2) = 102; mat34(0, 3) = 103;
  mat34(1, 0) = 110; mat34(1, 1) = 111; mat34(1, 2) = 112; mat34(1, 3) = 113;
  mat34(2, 0) = 120; mat34(2, 1) = 121; mat34(2, 2) = 122; mat34(2, 3) = 123;

  std::cout << "---------- Eigen mat 1" << std::endl;
  std::cout << mat53 << std::endl;
  std::cout << "---------- Eigen mat 2" << std::endl;
  std::cout << mat34 << std::endl;

  std::cout << "---------- Eigen mat1*mat2" << std::endl;
  std::cout << mat53 * mat34 << std::endl;

  std::cout << "---------- 自前実装 ------------" << std::endl;
  int I;
  float fmat53[5 * 3];
  I = 5;
  fmat53[matijI(0, 0, I)] = 00; fmat53[matijI(0, 1, I)] = 01; fmat53[matijI(0, 2, I)] = 02;
  fmat53[matijI(1, 0, I)] = 10; fmat53[matijI(1, 1, I)] = 11; fmat53[matijI(1, 2, I)] = 12;
  fmat53[matijI(2, 0, I)] = 20; fmat53[matijI(2, 1, I)] = 21; fmat53[matijI(2, 2, I)] = 22;
  fmat53[matijI(3, 0, I)] = 30; fmat53[matijI(3, 1, I)] = 31; fmat53[matijI(3, 2, I)] = 32;
  fmat53[matijI(4, 0, I)] = 40; fmat53[matijI(4, 1, I)] = 41; fmat53[matijI(4, 2, I)] = 42;

  float fmat34[3 * 4];
  I = 3;
  fmat34[matijI(0, 0, I)] = 100; fmat34[matijI(0, 1, I)] = 101; fmat34[matijI(0, 2, I)] = 102; fmat34[matijI(0, 3, I)] = 103;
  fmat34[matijI(1, 0, I)] = 110; fmat34[matijI(1, 1, I)] = 111; fmat34[matijI(1, 2, I)] = 112; fmat34[matijI(1, 3, I)] = 113;
  fmat34[matijI(2, 0, I)] = 120; fmat34[matijI(2, 1, I)] = 121; fmat34[matijI(2, 2, I)] = 122; fmat34[matijI(2, 3, I)] = 123;

  std::cout << "---------- 自前実装 mat 1" << std::endl;
  std::cout << mat_to_string(fmat53, "%3.0lf", 5, 3);

  std::cout << "---------- 自前実装 mat 2" << std::endl;
  std::cout << mat_to_string(fmat34, "%3.0lf", 3, 4);

  std::cout << "---------- 自前実装 mat1*mat2" << std::endl;
  float fmat54[5 * 4];
  mult_matrixIJK(fmat54, fmat53, fmat34, 5, 3, 4);
  std::cout << mat_to_string(fmat54, "%5.0lf", 5, 4);

}