実装

４階テンソルと２階テンソルの積を関数として以下のように実装する。

/* ４階テンソルと２階テンソルの積 */
void WH_Tensor3D__whT4_colon_whT_OUT_whT
(double whT4_C[3][3][3][3],
 double whT_X[3][3],
 double OUT__whT_Y[3][3])
/*
  入力引数:
    whT4_Cは４階テンソル
    whT_X はテンソル
  出力引数:
    OUT__whT_Yはテンソルで、whT4_C と whT_X の積
*/
{
  int i, j, k, l;

  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      double tmp = 0;

      for (k = 0; k < 3; k++) {
        for (l = 0; l < 3; l++) {
          tmp += whT4_C[i][j][k][l] * whT_X[k][l];
        }
      }
      OUT__whT_Y[i][j] = tmp;
    }
  }
}

/* ２階テンソルと４階テンソルの積 */
void WH_Tensor3D__whT_colon_whT4_OUT_whT
(double whT_X[3][3],
 double whT4_C[3][3][3][3],
 double OUT__whT_Y[3][3])
/*
  入力引数:
    whT_X はテンソル
    whT4_Cは４階テンソル
  出力引数:
    OUT__whT_Yはテンソルで、whT_X と whT4_C の積
*/
{
  int i, j, k, l;

  for (k = 0; k < 3; k++) {
    for (l = 0; l < 3; l++) {
      double tmp = 0;

      for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
          tmp += whT_X[i][j] * whT4_C[i][j][k][l];
        }
      }
      OUT__whT_Y[k][l] = tmp;
    }
  }
}