はじめに

　「数値計算に興味があるし，データのやり取りが簡単だと嬉しいけど，全ての動作が高速じゃないとイヤ！！！」で， matlab になかなか手を出せていない全国 $3\uparrow \uparrow \uparrow 3$ ¹ 人の皆さんこんにちは．このエントリは，高速計算のために取れる最終手段の 1 つ，mex function を紹介し，皆さんのハッピーエンジョイストレスフリーな数値計算ライフに寄り添うものです．また，mex function は octave でも利用できますが，ここではケアしません²．

mex function に，ひいては高速化に手を出す前に

　改めて，mex function は高速化における 最後の砦 です．求める計算が組み込み関数で提供されていないか，置き換える対象は十分にボトルネックであるかを，それぞれ精査する必要があります．

前者を満たさず，組み込み関数が存在した場合，開発に掛けた労力が無駄になりますし³，自分の実装より高速で性能が良かった，という結末もザラです．特に for 文，if 文が遅くなるのはインタプリタ型言語の宿命です．まずは組み込み演算子や組み込み関数で回避を試みるべきです．
後者を満たさない場合，掛けた労力に対しプログラム全体の実行時間への効果は薄く，やはり徒労に終わるはずです．並列計算における議論ですが，アムダールの法則は効果的な高速化について論じています．安易な高速化はただ可読性と保守性を引き下げるだけの 危険な行為 です．

これらの条件をパスした部分は高速化のため， mex function を用いる価値があります．matlab の高級で柔軟な仕組みを利用しつつ，ボトルネックを mex function で解消するような使い方が安全と思われます．また mex では，matlab の 強力なメモリの庇護から抜けます．メモリリークには十分注意が必要です⁴．

　ところで，matlab は強力である一方，大変高価なソフトウェアです．欲しい関数が，未購入 Toolbox に含まれていること，ありますよね．現場からは以上です． ~~（弊学はフルライセンスを全学にバラ撒いてくれ～～たのむ～～～～～）~~

mex function

　c や c++ による高速な関数を matlab から呼べる関数として提供する仕組みです⁵．matlab では，ソースファイルから mex function をビルドする関数 mex が用意されており，裏でコンパイル⁶を行います．c++ において，クラスを用いた記述方法もありますが，本エントリでは c 言語における記述を紹介します．

ゲートウェイルーチン

一般的な c プログラムが main 関数を持つように，mex でもエントリポイントとなる関数が存在します．

mexExample.c

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    /* write code */
}

この，mexFunction が，mex に必須な関数です．データのやり取りに関わる引数については次の通りです．

int nlhs: 出力（左辺）の数．後述の*plhs[]の，配列サイズと等価．
mxArray *plhs[]: 出力（左辺）のデータへのポインタの配列．複数出力のため，ポインタの配列．
int nrhs: 入力（右辺）引数の数．後述の*prhs[]の，配列サイズと等価．
mxArray *prhs[]: 入力（右辺）のデータへのポインタの配列．複数入力のため，ポインタの配列．

実例と plhs について

matlab から mex fuction を呼ぶ場合，ファイル名が関数の識別子として用いられます．

example.m

A = [1 2 3; 4 5 6; 7 8 9;];    % 3×3 の行列
b = 1.0;                       % スカラー

[x, y, z] = mexExample(A, b);  % mex function の呼び出し

このように matlab から呼び出された場合，引数にはそれぞれ次の値，参照（正しくはアドレス値）が入ります．

nlhs : 3
nrhs : 2
prhs[0] : 行列 A への参照
prhs[1] : スカラー b への参照（$1 \times 1$ 行列として取得できる⁷）

plhs[] に関しては，出力したい形状の mxArray を作成⁸（create）します．この例の場合，plhs[0]～plhs[2] のそれぞれに mxArray を作成します．横着して prhs[] を再利用することは不可能です⁹．

mxArray について

　ここで mxArray について掘り下げ，データを c の配列としてアクセスする方法を紹介します．前述の通り， mxArray は入出力のデータに対応する型です．matlab におけるデータとは，実数行列だけでなく，様々な表現を内包するものです．列挙すれば，文字列行列，論理行列，複素数行列，疎行列表現，構造体，cell 行列，関数ハンドル……，実に様々です¹⁰．また，次元数も次元の要素数も任意です．この状況に対して，mxArray のフィールドには，これらの情報を保持する仕組みがあり，また mxArray 作成時に設定できます．この情報は例えば，入力されたデータに関するエラーハンドルに役立ちます．

checkError.c

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    int mrows = mxGetM(prhs[0]);
    int ncols = mxGetN(prhs[0]);
    if( !mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) ||
    !(mrows==1 && ncols==1) ) {
        mexErrMsgIdAndTxt( "MATLAB:timestwo:inputNotRealScalarDouble",
                           "Input must be a noncomplex scalar double.");
    }

    /* write code */
}

このコードでは，入力された 1 つ目のデータが，実数で double のスカラーでない場合に，matlab へエラーを出します．
　次に肝心のデータへのアクセスは，mxGet* 関数を使います．* には，データの型に応じてキーワードが入ります¹¹．matlab で一般的な倍精度実数行列のデータへアクセスするには，mxGetDoubles を使います¹²．

accessCall.m

x = [1 2 3 4]';     % 転置して 要素 4 の列ベクトルとする． 
y = accessData(x);

accessData.c

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    double *vec = mxGetDoubles(prhs[0]);
    // i 番目の要素に vec[i] として参照できる．
    double sum = 0.0;
    for (int i = 0; i < 4; i++) {
        sum += vec[i];
    }
    /* write code */
}

多次元配列となかよく

　本エントリの本題に入ります．mex function における入出力について紹介しましたが，多次元行列のアクセスはどうなるでしょうか．まず試しに，二次元行列のアクセスを行います．なお，c や c++ で二次元配列を操作するとき，多くの場合行優先ですが，mex においては 列優先 である点に注意が必要です．

access2DCall.m

A = [1 2 3; 4 5 6];     % 2×3 の行列
B = access2DData(A);

access2DData.c

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    double *mat = mxGetDoubles(prhs[0]);
    // i 行目 j 列目の要素は mat[i+2*j] で参照できる．
    double sum = 0.0;
    for (int j = 0; j < 3; j++) {
        for (int i = 0; i < 2; i++) {
            sum += mat[i+2*j];
        }
    }
    /* write code */
}

mxGet*() は，任意形状の行列に対応するため，単純に行列の先頭ポインタのみ返します．そのため，プログラム書く人間が上手いこと元形状に即した変換を考える必要があります．
　……これ，多次元になると本当に厄介で，例えば 6 次元行列の表現を考えます¹³．適当な 6 次元行列として，X = rand(2,3,4,5,6,7) を与えたとき，X(i,j,k,l,m,n) の要素を mex で参照するには，次のように記述します（簡単のため，mex での識別子と共通に書きます）．

X[i+2*(j+3*(k+4*(l+5*(m+6*n))))]

これは何？長い，見づらい，バグの素って感じですね……

配列へのポインタでなんとかする

　そこで，~~無理くり~~多次元行列となかよしになりましょう．具体的には，X[n][m][l][k][j][i] という書き方を手に入れる方法です．matlab の次元の並びと逆順になることに注意です．

access6DCall.m

X = rand(2,3,4,5,6,7);     % 先程の 2×3×4×5×6×7 の行列
Y = access6DData(X);

access6DData.c

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    double *ten = mxGetDoubles(prhs[0]);
    double (*multen)[6][5][4][3][2] = (double(*)[6][5][4][3][2])ten; // これが配列へのポインタ
    double sum = 0.0;
    for (int n = 0; n < 7; n++) {
        for (int m = 0; m < 6; m++) {
            for (int l = 0; l < 5; l++) {
                for (int k = 0; k < 4; k++) {
                    for (int j = 0; j < 3; j++) {
                        for (int i = 0; i < 2; i++) {
                            sum += multen[n][m][l][k][j][i];         // バッチリ！
                        }
                    }
                }
            }
        }
    }
    /* write code */
}

以上， 配列へのポインタ¹⁴ でなんとかする方法です．先に取り上げた二次元行列にも使えます．
ex)：double (*mulmat)[2] = (double(*)[2])mat; → mulmat[j][i]

　また，配列へのポインタを使うには，コンパイル時に要素数が静的に定まっている必要があります．そのため，定数や const 修飾子付き変数などで設定します．しかし，mex コンパイラに gcc を選択した場合には，要素数を普通の変数で設定できました¹⁵．どうして？？？ ただし，macOS では mex のビルドに gcc が使えません¹⁶．（12/1 19:11 修正）
　また，例では要素数がコンパイル時に定まっている場合を取り上げましたが，要素数の指定に変数も扱えます．C++ としてビルドする場合ではケアされないので，注意です（n敗）．また，windows で msvc を用いてビルドする場合には，コンパイル時に数値が定まっている必要があります．適宜，minGW による gcc を用いてこの問題を回避しましょう．

最後に

　以上，mex function 入門として，簡単な入出力と，よく知る多次元配列表記のために配列へのポインタを用いる方法を紹介しました．matlab は非商用限定の個人利用ライセンスが15,000 円ほどで買えるようになりましたし，学生であれば 5000 円ほどで買えます．さらに，学生は，MATLAB and Simulink Student Suite （matlab と Simlink に加え，主要な Toolbox が 10 個も付いてくるお買い得パック）が購入できます．個人利用ライセンスの枠で同じ内容を揃えると 65,000 円ほどかかるところ，僅か 10,000 円 程度です．Web を通して購入すれば，同時購入する Toolbox が 3,000 円から 1,000 円程度に割り引かれるキャンペーンもあります．買いましょう． ~~（ダイマ楽しい）~~

　個人的には，matlab と mex function の布教ができたことはもちろん，後半部分の話ができて満足です．多次元配列へのアクセスで，何か上手い方法はあるだろうと思っていましたが，配列へのポインタを発掘したときはポインタのこと何も分からんってなりました．c は今や古典ですが，古典として楽しめるところが結構いい感じですね．
　
　また，このエントリは，信州大学のモノづくりサークル kstm（カスタム）のアドベントカレンダー一発目でした．kstm の構成員は自らのアンテナに従い思い思いに活動しているため，それぞれが何かしら尖ったモノを持っています．今年も実用的な技術ネタから，技術一辺倒ではない伸び伸びとした話が出揃うはずです．明日以降も引き続きお楽しみくださいませ．
　

2019年現在，国連の予測では，世界人口は 2100 年に 110 億人へ達し，その後は減少すると見込まれている． ↩
octave では，mex function と似た oct-file という仕組みも提供されている． ↩
教育などを目的とした再開発では，この限りではない（そんなことわざわざ mex でやるだろうか）． ↩
メモリリークは，最悪ブルスクを引き起こす（n 敗）． ↩
紹介した言語の他に Fortran も利用できる． ↩
このとき使用されるコンパイラは，任意のものに変更できる． ↩
mxGetScalar 関数を使い，明示的にスカラーとして値を読み込むことも可能． ↩
mxArray の create 周りの詳細情報はここを参照されたい． ↩
出力と全く同じ形状の行列を出力するとき，ディープコピー関数，mxDuplicateArray を使う手がある． ↩
これらを網羅するには別エントリができてしまうので，ここや，このあたりを参照されたい． ↩
詳しくはここ ↩
matlab R2018a 以降を想定している．matlab R2017b 以前は，(double*)mxGetData(prhs[0]) と呼ぶ（型に関してはキャストでなんとかする）． ↩
6 次元の行列なんて使うのか？って感じだが，独立なデータの列挙に次元を振ったり，matlab で for を回避するテクを使ったりすると，案外使う． ↩
間違っても配列のポインタではない（1 敗）．ポインタにはたくさんの種類がある．（このサイトが今回一番の収穫まである） ↩
Stack Overflow の回答者も，強調して，コンパイル時に明確になってるべきって言ってるんだけどなぁ．なんでだろ． ↩
Windows では MinGW を用いて，gcc を利用できる． ~~ただし，OpenMP に対応しないのが痛い．~~ （lcc-win64 コンパイラに関しての制限でした．MinGW の gcc で OpenMP を使用できます．）Linux であれば，ネイティブで gcc が使える．詳しい表はここ． ↩

mex function ことはじめ ～N次元配列となかよし編～