fulltosparse.c

这是机械工业张德丰《MATLAB与外部程序接口编程》随书源码

/* fulltosparse.c */
#include <math.h>  //该头文件包含有ceil( )的函数原型
#include "mex.h"
#if defined (NAN_EQUALS_ZERO)
#if defined 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 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.");
	}
	//获取输入数据的指针和大小
	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=(int)ceil((double)m *(double)n *percent_sparse);
	plhs[0]=mxCreateSparse(m,n,nzmax,cmplx);
	sr=mxGetPr(plhs[0]);
	si=mxGetPr(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元素在，若不就将percent_sparse
				//增加10%，并且重新计算nzmax和augment的值
				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]=pr9i0;
				if(cmplx)
				{
					si[k]=pi[i];
				}
				irs[k]=i;
				k++;
			}
		}
		pr+=m;
		pi+=n;
	}
	jcs[n]=k;
}