lasso经典matlab源码

/* mex -c mexutil.c mex repmat.c mexutil.obj to check for warnings: gcc -Wall -I/cygdrive/c/MATLAB6p1/extern/include -c repmat.c */ #include "mexutil.h" #include <string.h> /* repeat a block of memory rep times */ void memrep(char *dest, size_t chunk, int rep) { #if 0 /* slow way */ int i; char *p = dest; for(i=1;i<rep;i++) { p += chunk; memcpy(p, dest, chunk); } #else /* fast way */ if(rep == 1) return; memcpy(dest + chunk, dest, chunk); if(rep & 1) { dest += chunk; memcpy(dest + chunk, dest, chunk); } /* now repeat using a block twice as big */ memrep(dest, chunk<<1, rep>>1); #endif } void repmat(char *dest, const char *src, int ndim, int *destdimsize, int *dimsize, const int *dims, int *rep) { int d = ndim-1; int i, chunk; /* copy the first repetition into dest */ if(d == 0) { chunk = dimsize[0]; memcpy(dest,src,chunk); } else { /* recursively repeat each slice of src */ for(i=0;i<dims[d];i++) { repmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1], ndim-1, destdimsize, dimsize, dims, rep); } chunk = destdimsize[d-1]*dims[d]; } /* copy the result rep-1 times */ memrep(dest,chunk,rep[d]); } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { const mxArray *srcmat; int ndim, *dimsize, eltsize; const int *dims; int ndimdest, *destdims, *destdimsize; char *src, *dest; int *rep; int i,nrep; int extra_rep = 1; int empty; if(nrhs < 2) mexErrMsgTxt("Usage: xrepmat(A, [N M ...])"); srcmat = prhs[0]; if(mxIsSparse(srcmat)) { mexErrMsgTxt("Sorry, can't handle sparse matrices yet."); } if(mxIsCell(srcmat)) { mexErrMsgTxt("Sorry, can't handle cell arrays yet."); } ndim = mxGetNumberOfDimensions(srcmat); dims = mxGetDimensions(srcmat); eltsize = mxGetElementSize(srcmat); /* compute dimension sizes */ dimsize = mxCalloc(ndim, sizeof(int)); dimsize[0] = eltsize*dims[0]; for(i=1;i<ndim;i++) dimsize[i] = dimsize[i-1]*dims[i]; /* determine repetition vector */ ndimdest = ndim; if(nrhs == 2) { nrep = mxGetN(prhs[1]); if(nrep > ndimdest) ndimdest = nrep; rep = mxCalloc(ndimdest, sizeof(int)); for(i=0;i<nrep;i++) { double repv = mxGetPr(prhs[1])[i]; rep[i] = (int)repv; } if(nrep == 1) { /* special behavior */ nrep = 2; rep[1] = rep[0]; } } else { nrep = nrhs-1; if(nrep > ndimdest) ndimdest = nrep; rep = mxCalloc(ndimdest, sizeof(int)); for(i=0;i<nrep;i++) { rep[i] = (int)*mxGetPr(prhs[i+1]); } } for(i=nrep;i<ndimdest;i++) rep[i] = 1; /* compute output size */ destdims = mxCalloc(ndimdest, sizeof(int)); for(i=0;i<ndim;i++) destdims[i] = dims[i]*rep[i]; for(;i<ndimdest;i++) { destdims[i] = rep[i]; extra_rep *= rep[i]; } destdimsize = mxCalloc(ndim, sizeof(int)); destdimsize[0] = eltsize*destdims[0]; for(i=1;i<ndim;i++) destdimsize[i] = destdimsize[i-1]*destdims[i]; /* for speed, array should be uninitialized */ plhs[0] = mxCreateNumericArray(ndimdest, destdims, mxGetClassID(srcmat), mxIsComplex(srcmat)?mxCOMPLEX:mxREAL); /* if any rep[i] == 0, output should be empty array. Added by KPM 11/13/02. */ empty = 0; for (i=0; i < nrep; i++) { if (rep[i]==0) empty = 1; } if (empty) return; src = (char*)mxGetData(srcmat); dest = (char*)mxGetData(plhs[0]); repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); if(mxIsComplex(srcmat)) { src = (char*)mxGetPi(srcmat); dest = (char*)mxGetPi(plhs[0]); repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); } }