Go to Qiita Advent Calendar 2024 Top
LoginSignup
17
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

信州大学 kstmAdvent Calendar 2019

Day 1
@iyselee07

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

Last updated at Posted at 2019-11-30

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

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

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

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

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

mex function

　c や c++ による高速な関数を matlab から呼べる関数として提供する仕組みです5．matlab では，ソースファイルから mex function をビルドする関数 mex が用意されており，裏でコンパイル6を行います．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$ 行列として取得できる7）

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

mxArray について

　ここで mxArray について掘り下げ，データを c の配列としてアクセスする方法を紹介します．前述の通り， mxArray は入出力のデータに対応する型です．matlab におけるデータとは，実数行列だけでなく，様々な表現を内包するものです．列挙すれば，文字列行列，論理行列，複素数行列，疎行列表現，構造体，cell 行列，関数ハンドル……， 実に様々です10．また，次元数も次元の要素数も任意です．この状況に対して，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* 関数を使います．* には，データの型に応じてキーワードが入ります11．matlab で一般的な倍精度実数行列のデータへアクセスするには，mxGetDoubles を使います12．

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 次元行列の表現を考えます13．適当な 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 */
}

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

　また，配列へのポインタを使うには，コンパイル時に要素数が静的に定まっている必要があります．そのため，定数や const 修飾子付き変数などで設定します．しかし，mex コンパイラに gcc を選択した場合には，要素数を普通の変数で設定できました15．どうして？？？ ただし，macOS では mex のビルドに gcc が使えません16． （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 の構成員は自らのアンテナに従い思い思いに活動しているため，それぞれが何かしら尖ったモノを持っています．今年も実用的な技術ネタから，技術一辺倒ではない伸び伸びとした話が出揃うはずです．明日以降も引き続きお楽しみくださいませ．
　

  1. 2019年現在，国連の予測では，世界人口は 2100 年に 110 億人へ達し，その後は減少すると見込まれている． ↩

  2. octave では，mex function と似た oct-file という仕組みも提供されている． ↩

  3. 教育などを目的とした再開発では，この限りではない（そんなことわざわざ mex でやるだろうか:thinking:）． ↩

  4. メモリリークは， 最悪ブルスクを引き起こす（n 敗）． ↩

  5. 紹介した言語の他に Fortran も利用できる． ↩

  6. このとき使用されるコンパイラは，任意のものに変更できる． ↩

  7. mxGetScalar 関数を使い，明示的にスカラーとして 値 を読み込むことも可能． ↩

  8. mxArray の create 周りの詳細情報は ここ を参照されたい． ↩

  9. 出力と全く同じ形状の行列を出力するとき，ディープコピー関数，mxDuplicateArray を使う手がある． ↩

  10. これらを網羅するには別エントリができてしまうので，ここ や，このあたりを参照されたい． ↩

  11. 詳しくは ここ ↩

  12. matlab R2018a 以降を想定している．matlab R2017b 以前は，(double*)mxGetData(prhs[0]) と呼ぶ（型に関してはキャストでなんとかする）． ↩

  13. 6 次元の行列なんて使うのか？って感じだが，独立なデータの列挙に次元を振ったり，matlab で for を回避するテクを使ったりすると，案外使う． ↩

  14. 間違っても配列のポインタではない（1 敗）．ポインタにはたくさんの種類がある．（このサイトが今回一番の収穫まである） ↩

  15. Stack Overflow の回答者も，強調して，コンパイル時に明確になってるべきって言ってるんだけどなぁ．なんでだろ． ↩

  16. Windows では MinGW を用いて，gcc を利用できる． ただし，OpenMP に対応しないのが痛い． （lcc-win64 コンパイラに関しての制限でした．MinGW の gcc で OpenMP を使用できます．）Linux であれば，ネイティブで gcc が使える．詳しい表はここ． ↩

17
8
5

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
17
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?