行列同士の掛け算_clapack

行列同士の掛け算

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

ここでは例として、行列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行列の管理とあわせて関数化しておくとよいと思います。