例11-8.m

这是一个MATLAB7.0基础与提高例题的所有源码

/* fulltosparse.c */
#include <math.h> 
#include "mex.h"

/* 如果读者使用的编译器使用符号NaN表示0，那么在编译此例子
* 代码的时候读者必须使用标志-DNAN_EQUALS_ZERO ，例如：
* mex -DNAN_EQUALS_ZERO findnz.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[]
        )
{
  /* 声明变量 */
  int j,k,m,n,nzmax,*irs,*jcs,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);

  /* 为稀疏矩阵分配内存空间 
   * 注意: 当最多有20% 的数据为空的时候，使用ceil
   *  函数将导致舍入。 
   */

  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]);
    
  /* 复制非0值 */
  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]))) {

        /* 检查非0元素是否适合分配输出数组，
         * 如果不适合，将稀疏矩阵比例增加， 
         * 重新计算nzmax，并且扩大稀疏矩阵。
         * 此过程一直持续到10％为止
*/
        if (k >= nzmax) {
          int oldnzmax = nzmax;
          percent_sparse += 0.1;
          nzmax = (int)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(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;
}