外部インタフェース/API

スパース配列の操作

MATLAB APIは、ユーザのMEX-ファイル内部からスパース配列を作成し操作できる関数を提供しています。これらのAPIルーチンは、スパース配列に関連する2つのパラメータirとjcにアクセスしたり操作します。MATLABのスパース配列の保存法は、「MATLAB配列」を参照してください。

つぎの例は、スパース配列の作成方法を説明しています。

/
*==============================================================
* fulltosparse.c
* This example demonstrates how to populate a sparse
* matrix.  For the purpose of this example, you must pass in a
* non-sparse 2-dimensional argument of type double.

* Comment: You might want to modify this MEX-file so that you can 
* use it to read large sparse data sets into MATLAB.
*
* This is a MEX-file for MATLAB.  
* Copyright (c) 1984-2000 The MathWorks, Inc.
* 
*=============================================================*
/

/* $ Revision: 1.5 $ */

#include <math.h> /* Needed for the ceil() prototype. */
#include "mex.h"

/* If you are using a compiler that equates NaN to be zero, you 
 * must compile this example using the flag  -DNAN_EQUALS_ZERO.
 * For example:
 *
 *     mex -DNAN_EQUALS_ZERO fulltosparse.c
 *
 * This will correctly define the IsNonZero macro for your C
 * compiler.
 */

#if defined(NAN_EQUALS_ZERO)
#define IsNonZero(d) ((d)!=0.0 || mxIsNaN(d))
#else
#define IsNonZero(d) ((d)!=0.0)
#endif

void mexFunction(
         int nlhs,       mxArray *plhs[],
         int nrhs, const mxArray *prhs[]
         )
{
    /* Declare variable. */
    int j,k,m,n,nzmax,*irs,*jcs,cmplx,isfull;
    double *pr,*pi,*si,*sr;
    double percent_sparse;
    
    /* Check for proper number of input and output arguments. */    
    if (nrhs != 1) {
    mexErrMsgTxt("One input argument required.");
    } 
    if(nlhs > 1){
    mexErrMsgTxt("Too many output arguments.");
    }
    
    /* Check data type of input argument. */
    if (!(mxIsDouble(prhs[0]))){
    mexErrMsgTxt("Input argument must be of type double.");
    }
    
    if (mxGetNumberOfDimensions(prhs[0]) != 2){
    mexErrMsgTxt("Input argument must be two dimensional\n");
    }

    /* Get the size and pointers to input data. */
    m  = mxGetM(prhs[0]);
    n  = mxGetN(prhs[0]);
    pr = mxGetPr(prhs[0]);
    pi = mxGetPi(prhs[0]);
    cmplx = (pi==NULL ? 0 : 1);
    
    /* Allocate space for sparse matrix. 
     * NOTE:  Assume at most 20% of the data is sparse.  Use ceil
     * to cause it to round up. 
     */

    percent_sparse = 0.2;
    nzmax=(int)ceil((double)m*(double)n*percent_sparse);

    plhs[0] = mxCreateSparse(m,n,nzmax,cmplx);
    sr  = mxGetPr(plhs[0]);
    si  = mxGetPi(plhs[0]);
    irs = mxGetIr(plhs[0]);
    jcs = mxGetJc(plhs[0]);
    
    /* Copy nonzeros. */
    k = 0; 
    isfull=0;
    for (j=0; (j<n ); j++) {
      int i;
      jcs[j] = k;
      for (i=0; (i<m ); i++) {
    if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {

      /* Check to see if non-zero element will fit in 
       * allocated output array.  If not, increase percent_sparse 
       * by 10%, recalculate nzmax, and augment the sparse array.
       */
      if (k>=nzmax){
        int oldnzmax = nzmax;
        percent_sparse += 0.1;
        nzmax = (int)ceil((double)m*(double)n*percent_sparse);

        /* Make sure nzmax increases atleast by 1. */
        if (oldnzmax == nzmax) 
          nzmax++;

        mxSetNzmax(plhs[0], nzmax); 
        mxSetPr(plhs[0], mxRealloc(sr, nzmax*sizeof(double)));
        if(si != NULL)
          mxSetPi(plhs[0], mxRealloc(si, nzmax*sizeof(double)));
        mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));
       
        sr  = mxGetPr(plhs[0]);
        si  = mxGetPi(plhs[0]);
        irs = mxGetIr(plhs[0]);
      }
      sr[k] = pr[i];
      if (cmplx){
        si[k]=pi[i];
      }
      irs[k] = i;
      k++;
        }    
      }
      pr += m;
      pi += m;
    }
    jcs[n] = k;
}

MATLABプロンプトで

full = eye(5)
full =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1

と入力すると、5行5列のフル単位行列を作成します。フル行列に対してfulltosparseを使うと、対応するスパース行列を生成します。

spar = fulltosparse(full)
spar =
   (1,1)        1
   (2,2)        1
   (3,3)        1
   (4,4)        1
   (5,5)        1

多次元数値配列の総さ

C MEX-ファイルからの関数の呼び出し