!!!行列同士の掛け算

ここでは、行列AがMa x Na、行列BがMb x Nbで構成されるときの乗算について
プログラム主体で書きます。

{{math A = \begin{bmatrix} a00&a01&a02\\ a10&a11&a12 \end{bmatrix} }}

{{math B = \begin{bmatrix} b00&b01&b02&b03\\ b10&b11&b12&b13\\ b20&b21&b22&b23\\ \end{bmatrix} }}

ここでは例として、行列Aを2 x 3、行列Bを3 x 4とします。
「行列」はその名の通り行と列で構成されるもので、「M x N」とすると行はM個、列はN個、となります。

二次元配列で書くとすると、
 A[0][0] = a00;
 A[0][1] = a01;
 A[0][2] = a02;
 A[1][0] = a10;
 A[1][1] = a11;
 A[1][2] = a12;
のようになります（配列添え字の前が行のインデックス、後ろが列のインデックス）。

さて、CLAPACKでは行と列を逆転して扱います。上の行列Aを例にとると、
 double matA[3 * 2];
 matA[0 * 2 + 0] = a00;
 matA[0 * 2 + 1] = a10;
 matA[1 * 2 + 0] = a01;
 matA[1 * 2 + 1] = a11;
 matA[2 * 2 + 0] = a02;
 matA[2 * 2 + 1] = a12;
と格納されます。行列Bも同じくです。

CLAPACKを使った計算は、すべてにおいて列行の順の一次配列を行列として使うので、入り口と出口でコンバート（転置）するようにして、内部では一貫して列行で扱うように心がけるようにします（でないとこんがらがるので）。

!!MN行列の掛け算の法則

行列AがMa x Na、行列BがMb x Nbで構成される場合、計算後の行列Cは
Ma x Nbの行と列になります。これは
 Ma x Na   Mb x Nb
の左端と右端を採用した、と覚えるといいかもしれません。
また、「Na = Mb」である必要があります。上の真ん中2つはイコール、と覚えることができます。

たとえば、行列Aが2 x 3、行列Bが3 x 4の場合は、乗算後の行列は（両端をとって）2 x 4となることになります。

行列乗算くらいならCLAPACKを使う必要もありませんね。
ということでプログラムです。

!!行列同士の掛け算プログラム

 /**
  * CLAPACKの行列同士の乗算（列と行が逆になっている点に注意）。
  * @param[in]  pM1    行列M1
  * @param[in]  m1M    行列M1での縦方向の要素数
  * @param[in]  m1N    行列M1での横方向の要素数
  * @param[in]  pM2    行列M2
  * @param[in]  m2M    行列M2での縦方向の要素数
  * @param[in]  m2N    行列M2での横方向の要素数
  * @param[out] pRetM  m1M x m2N 行列が返る
  */
 bool MatrixMultiply(double *pM1, const int m1M, const int m1N,
     double *pM2, const int m2M, const int m2N,
     double *pRetM)
 {
     int i, j, k;
     int iPos, iPos2, iPos3, iRPos;
     double dVal;
 
     if(m1N != m2M) return false;
 
     for(i = 0; i < m1M; i++) {
         iRPos = i;
         iPos2 = 0;
         for(j = 0; j < m2N; j++) {
             iPos  = i;
             iPos3 = iPos2;
             dVal = 0.0;
             for(k = 0; k < m1N; k++) {
                 dVal += (*(pM1 + iPos)) * (*(pM2 + iPos3));
                 iPos += m1M;
                 iPos3++;
             }
             *(pRetM + iRPos) = dVal;
 
             iPos2 += m2M;
             iRPos += m1M;
         }
     }
 
     return true;
 }

pM1とpM2にて、行と列が入れ替わっているCLAPACK用の行列要素を入れます。
結果がpRetMの配列に返ります。

代数計算系のプログラムでは、この行列乗算が結構出てきますので、MN行列の管理とあわせて関数化しておくとよいと思います。