fulltosparse.c

matlab实用教程

#include <math.h>
#include "mex.h"
 
/* 若使用的编译器中NaN等同于0，那么编译此文件时，必须按照下面的格式
 *     mex -DNAN_EQUALS_ZERO fulltosparse.c
*/
 
#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[]
         )
{
    /* 声明变量 */
    mwSize m,n;
    mwSize nzmax;
    mwIndex *irs,*jcs,j,k;
    int cmplx,isfull;
    double *pr,*pi,*si,*sr;
    double percent_sparse;
    
    /*检查输入/输出参数的个数是否合法 */    
    if (nrhs != 1) {
    mexErrMsgTxt("One input argument required.");
    } 
    if(nlhs > 1){
    mexErrMsgTxt("Too many output arguments.");
    }
    
    /* 检查输入参数的数据类型  */
    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");
    }
 
    /* 获取输入数据的大小和指针 */
    m  = mxGetM(prhs[0]);
    n  = mxGetN(prhs[0]);
    pr = mxGetPr(prhs[0]);
    pi = mxGetPi(prhs[0]);
    cmplx = (pi==NULL ? 0 : 1);
    
    /* 为稀疏矩阵分配内存  */ 
    percent_sparse = 0.2;
    nzmax=(mwSize)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]);
    
    /* 复制非零值 */
    k = 0; 
    isfull=0;
    for (j=0; (j<n); j++) {
      mwSize i;
      jcs[j] = k;
      for (i=0; (i<m ); i++) {
    if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {
 
      /* 检查非零元素是否适合分配的输出数组，否则，percent_sparse增加10%，
* 重新计算nzmax，然后扩大疏矩阵
       */
      if (k>=nzmax){
        mwSize oldnzmax = nzmax;
        percent_sparse += 0.1;
        nzmax = (mwSize)ceil((double)m*(double)n*percent_sparse);
 
        /* 确保nzmax至少增加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(mwIndex)));
       
        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;
}