9

MATLAB で組み込みの LAPACK/BLAS ルーチンを呼び出す方法を知りたいです。私は MATLAB と mex ファイルの経験がありますが、実際には LAPACK または BLAS ライブラリを呼び出す方法がわかりません。このような関数の mex ファイルを作成する必要がないため、呼び出しを簡素化するファイル交換のゲートウェイ ルーチンを見つけました。MATLAB とこれらの組み込みライブラリとの間の基本的なメッセージングを学習するには、おもちゃの例が必要です。行列の乗算や LU 分解などのおもちゃの例は大歓迎です。

4

2 に答える 2

10

上記のFEX送信からのlapack.mファイルの内部を見ると、関数の使用方法に関するいくつかの例が表示されます。

例：DGESVDを使用したSVD分解：

X = rand(4,3);
[m,n] = size(X);
C = lapack('dgesvd', ...
     'A', 'A', ...           % compute ALL left/right singular vectors
      m, n, X, m, ...        % input MxN matrix
      zeros(n,1), ...        % output S array
      zeros(m), m, ...       % output U matrix
      zeros(n), n, ....      % output VT matrix
      zeros(5*m,1), 5*m, ... % workspace array
      0 ...                  % return value
);
[s,U,VT] = C{[7,8,10]};      % extract outputs
V = VT';

（注：出力変数にこれらのダミー変数を使用した理由は、Fortran関数はすべての引数が参照によって渡されることを期待しているためですが、MATLABのMEX関数では入力を変更できないため、すべての入力のコピーを返すように記述されています。変更を加えたセル配列）

我々が得る：

U =
     -0.44459      -0.6264     -0.54243       0.3402
     -0.61505     0.035348      0.69537      0.37004
     -0.41561     -0.26532      0.10543     -0.86357
     -0.50132      0.73211     -0.45948    -0.039753
s =
       2.1354
      0.88509
      0.27922
V =
     -0.58777      0.20822     -0.78178
      -0.6026     -0.75743      0.25133
     -0.53981      0.61882      0.57067

これは、MATLAB独自のSVD関数と同等です。

[U,S,V] = svd(X);
s = diag(S);

それは与える：

U =
     -0.44459      -0.6264     -0.54243       0.3402
     -0.61505     0.035348      0.69537      0.37004
     -0.41561     -0.26532      0.10543     -0.86357
     -0.50132      0.73211     -0.45948    -0.039753
s =
       2.1354
      0.88509
      0.27922
V =
     -0.58777      0.20822     -0.78178
      -0.6026     -0.75743      0.25133
     -0.53981      0.61882      0.57067

編集：

完全を期すために、DGESVDルーチンのFortranインターフェースを直接呼び出すMEX関数の例を以下に示します。

良いニュースは、MATLABがライブラリlibmwlapackとlibmwblas2つの対応するヘッダーファイルblas.hを提供し、lapack.h使用できることです。実際、ドキュメントには、MEXファイルからBLAS/LAPACK関数を呼び出すプロセスを説明するページがあります。

この場合、lapack.h次のプロトタイプを定義します。

extern void dgesvd(char *jobu, char *jobvt, 
  ptrdiff_t *m, ptrdiff_t *n, double *a, ptrdiff_t *lda,
  double *s, double *u, ptrdiff_t *ldu, double *vt, ptrdiff_t *ldvt,
  double *work, ptrdiff_t *lwork, ptrdiff_t *info);

svd_lapack.c

#include "mex.h"
#include "lapack.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex m, n, lwork, info=0;
    double *A, *U, *S, *VT, *work;
    double workopt = 0;
    mxArray *in;

    /* verify input/output arguments */
    if (nrhs != 1) {
        mexErrMsgTxt("One input argument required.");
    }
    if (nlhs > 3) {
        mexErrMsgTxt("Too many output arguments.");
    }
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) {
        mexErrMsgTxt("Input matrix must be real double matrix.");
    }

    /* duplicate input matrix (since its contents will be overwritten) */
    in = mxDuplicateArray(prhs[0]);

    /* dimensions of input matrix */
    m = mxGetM(in);
    n = mxGetN(in);

    /* create output matrices */
    plhs[0] = mxCreateDoubleMatrix(m, m, mxREAL);
    plhs[1] = mxCreateDoubleMatrix((m<n)?m:n, 1, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(n, n, mxREAL);

    /* get pointers to data */
    A = mxGetPr(in);
    U = mxGetPr(plhs[0]);
    S = mxGetPr(plhs[1]);
    VT = mxGetPr(plhs[2]);

    /* query and allocate the optimal workspace size */
    lwork = -1;
    dgesvd("A", "A", &m, &n, A, &m, S, U, &m, VT, &n, &workopt, &lwork, &info);
    lwork = (mwSignedIndex) workopt;
    work = (double *) mxMalloc(lwork * sizeof(double));

    /* perform SVD decomposition */
    dgesvd("A", "A", &m, &n, A, &m, S, U, &m, VT, &n, work, &lwork, &info);

    /* cleanup */
    mxFree(work);
    mxDestroyArray(in);

    /* check if call was successful */
    if (info < 0) {
        mexErrMsgTxt("Illegal values in arguments.");
    } else if (info > 0) {
        mexErrMsgTxt("Failed to converge.");
    }
}

64ビットWindowsでは、MEXファイルを次のようにコンパイルします。mex -largeArrayDims svd_lapack.c "C:\Program Files\MATLAB\R2013a\extern\lib\win64\microsoft\libmwlapack.lib"

これがテストです：

>> X = rand(4,3);
>> [U,S,VT] = svd_lapack(X)
U =
   -0.5964    0.4049    0.6870   -0.0916
   -0.3635    0.3157   -0.3975    0.7811
   -0.3514    0.3645   -0.6022   -0.6173
   -0.6234   -0.7769   -0.0861   -0.0199
S =
    1.0337
    0.5136
    0.0811
VT =
   -0.6065   -0.5151   -0.6057
    0.0192    0.7521   -0.6588
   -0.7949    0.4112    0.4462

対。

>> [U,S,V] = svd(X);
>> U, diag(S), V'
U =
   -0.5964    0.4049    0.6870    0.0916
   -0.3635    0.3157   -0.3975   -0.7811
   -0.3514    0.3645   -0.6022    0.6173
   -0.6234   -0.7769   -0.0861    0.0199
ans =
    1.0337
    0.5136
    0.0811
ans =
   -0.6065   -0.5151   -0.6057
    0.0192    0.7521   -0.6588
   -0.7949    0.4112    0.4462

U（との固有ベクトルの符号Vは任意であるため、2つを比較すると符号が反転する可能性があることに注意してください）

于 2011-06-22T13:42:11.233 に答える
-2

matlab で SVD を使用する方法の良い/実用的な例は、ここで説明されています: キャプチャされた座標を画面座標に変換する

lapackを使用してobjective-cでsvdを計算する方法の詳細は、ここに書かれています xcodeのLAPACKからSVDを使用して、objective-CでA = USVtからVを計算します

于 2012-08-09T21:23:07.140 に答える